*Article* **Case Study of Heat Transfer during Artificial Ground Freezing with Groundwater Flow**

**Rui Hu 1,†, Quan Liu 2,\*,† and Yixuan Xing <sup>2</sup>**


Received: 13 August 2018; Accepted: 21 September 2018; Published: 25 September 2018

**Abstract:** For the artificial ground freezing (AGF) projects in highly permeable formations, the effect of groundwater flow cannot be neglected. Based on the heat transfer and seepage theory in porous media with the finite element method, a fully coupled numerical model was established to simulate the changes of temperature field and groundwater flow field. Firstly, based on the classic analytical solution for the frozen temperature field, the model's ability to solve phase change problems has been validated. In order to analyze the influences of different parameters on the closure time of the freezing wall, we performed the sensitivity analysis for three parameters of this numerical model. The analysis showed that, besides the head difference, the thermal conductivity of soil grain and pipe spacing are also the key factors that control the closure time of the frozen wall. Finally, a strengthening project of a metro tunnel with AGF method in South China was chosen as a field example. With the finite element software COMSOL Multiphysics® (Stockholm, Sweden), a three-dimensional (3D) numerical model was set up to simulate the change of frozen temperature field and groundwater flow field in the project area as well as the freezing process within 50 days. The simulation results show that the freezing wall appears in an asymmetrical shape with horizontal groundwater flow normal to the axial of the tunnel. Along the groundwater flow direction, freezing wall forms slowly and on the upstream side the thickness of the frozen wall is thinner than that on the downstream side. The actual pipe spacing has an important influence on the temperature field and closure time of the frozen wall. The larger the actual pipe spacing is, the slower the closing process will be. Besides this, the calculation for the average temperature of freezing body (not yet in the form of a wall) shows that the average temperature change of the freezing body coincides with that of the main frozen pipes with the same trend.

**Keywords:** artificial ground freezing method; groundwater flow; temperature field; freezing wall; effective hydraulic conductivity

#### **1. Introduction**

With its good water barrier and non-pollution characteristics, the artificial ground freezing (AGF) method is widespread in environmental engineering (e.g., excavation of pollutants) [1,2], mining engineering (e.g., mine shafts construction), and especially in underground engineering (e.g., subsea tunnel and urban subway construction) [3–6]. The basic principle of this method is to circulate a fluid refrigerant (ca. −30 ◦C) or liquid nitrogen (ca. −196 ◦C) through a pre-buried pipe network in the subsurface in order to form a freezing wall for strengthening construction or blocking groundwater flow.

Since its first application in Germany in 1887, the AGF method has been considerably developed in theory and practice. Especially in the conventional groundwater environment, for example, in freshwater and low flow rate conditions, this technology has been successfully applied in a large

number of projects and it tends to be mature [7–10]. However, when this method is implemented under some special groundwater conditions, such as high groundwater flow velocity area, seawater invasion area, or groundwater heavily polluted area, the freezing effect will be affected by groundwater flow velocity, salt content, and different types of pollutants [11–15]. For instance, in the frozen wall project of the Fukushima Power Plant, the "Ice wall" failed to prevent groundwater infiltration. The most likely cause was the complex hydrological environment and the high concentration of radioactive toxic water [1]. This will bring new challenges to the typical application of the AGF method. In this paper, we mainly focus on the effect of flow velocity on the freezing effect from the perspective of groundwater flow.

During the development of freezing walls, groundwater flow and formation temperature distribution will interact with each other [11,16–22]. On the one hand, groundwater flow will take away the heat (coldness) around the freezing columns and slow down the cooling process of the formation. On the other hand, with the lowering of the formation temperature, the groundwater flow process will also be impacted. After the formation temperature reaches freezing point, an ice-water phase transition will occur in the pores, which will result in the formation of ice. Thereby, this will cause changes in the physical properties of the soil, such as permeability, density, thermal conductivity, and strength. The higher the flow velocity is, the stronger that the interaction will be. As the flow velocity increases, the heat transferred by the groundwater flow gradually replaces the heat extracted by the freezing pipes, slowing down the growth of the freezing columns around individual freezing pipes. This is a potential challenge for the AGF method.

In the past, many researchers have carried out research on the effect of groundwater flow on AGF method. First, some scholars have put forward corresponding numerical models from the viewpoint of heat-fluid coupling. Harlan (1973) [23] presented a mathematical model that coupled heat and fluid flow based on analogy and applied to the analysis of heat-fluid phenomena in partially frozen soil. Due to the lack of reliable data, quantitative comparisons of observed and simulated profile responses were not made. Using a finite element scheme, Guymon and Luthin (1974) [24] developed a one-dimensional a coupled model for heat and moisture based on the Richards equation and the heat conduction equation, including convective components, and applied this to analyze moisture movement and storage in Arctic soils. He pointed out that the feasibility of the model depends primarily on the assumption of equations and the computation costs. Hansson et al. (2004) [25] presented a new method that accounts for phase change in a fully implicit numerical model. This model was evaluated by comparing predictions with data from laboratory columns freezing experiments and applied this to simulate hypothetical road problems involving simultaneous heat transport and water flow. Based on some parametric assumptions, McKenzie et al. (2007) [26] simulated groundwater flow with energy transport numerically by modifying the SUTRA computer code and applied this to freezing processes in peat bogs. The model was validated by the analytical solution from Lunardini (1985) without an advection term [27]. Kurylyk et al. (2014) [28] reanalyzed this analytical solution and compared the results between the numerical solution and the Lunardini analytical solution with advection term. In the study of Scheidegger et al. (2017) [29], the development of permafrost was simulated with a two-dimensional (2D) thermo-hydro (TH) model using COMSOL software and this model was validated by this analytical solution. Tan et al. (2011) [30] established a TH coupling model based on the theories of continuum mechanics, thermodynamics, and segregation potential. The model considered the Soret effect, as well as the influence of segregation potential on groundwater flow velocity and water pressure distribution. This model has been verified by comparison with the well-known TH coupling laboratory test that was conducted by Mizoguchi (1990) [31].

Secondly, some scholars have focused on experimental studies of the AGF method combined with groundwater flow. Anton Sers (2009) [32] studied the growth of the freezing body with various groundwater flow velocities through a small-scale sandbox. Pimentel et al. (2012) [11] developed a large-scale physical experimental device for the AGF with high groundwater flow velocities and compared the experimental results of different experimental setups [33–35]. The results showed that, for the conditions of no groundwater flow or single freezing pipe, there were no substantial differences with respect to temperature change. But, with multiple pipes and high groundwater flow velocity, the discrepancies become larger as advection becomes more dominant in relation to heat conduction. Huang et al. (2012) [36] conducted a series of physical experiments and numerical simulations with controlled groundwater flow and boundary conditions. These results have shown that the growth of a freezing body is affected by groundwater flow.

It has been proven that the groundwater flow will affect the temperature distribution around the freezing pipes and further prevent the formation of the freezing wall. However, the quantification of the effect of groundwater flow velocity on the evolution of the freezing wall is still in dispute. Jessberger and Jagow (1996) [37] studied the effect of groundwater flow velocity on the AGF method with brine circulation and pointed to a flow velocity of 2 m/day as the critical limit. That is to say, the frozen body between the pipes will not close when the flow rate is more than 2 m/day. Besides this, Andersland and Ladanyi (2004) [38] claimed that, even with flow velocities between 1 and 2 m/day, no closure of the freezing body will occur. Other scholars have also given some different critical flow velocities [39,40].

In this work, we introduced the parameter of effective hydraulic conductivity and the saturation function of unfrozen water. Based on this parameter and function, we set up a fully coupled three-dimensional (3D) numerical model (COMSOL Multiphysics® model with finite element) with groundwater flow field and frozen temperature field. This model is validated through the comparison of the results from numerical calculation and the analytical solution of the classic frozen problem from Lunardini [27] Its ability to solve the phase change problem is evaluated. With different numerical case studies, sensitivity analysis was performed for three parameters: the pipe spacing, thermal conductivity of soil grain and groundwater flow velocity (derived from hydraulic gradient or head difference). The corresponding influences of each parameter on the closure time of the freezing wall was analyzed. Finally, a strengthening project of a metro tunnel with AGF method in South China was chosen as an example for the field study. A numerical model was set up to simulate the artificial ground freezing process of this project. The regular pattern of frozen wall development was analyzed and the model was validated through the comparison of calculated and in-situ measured temperatures.

#### **2. Methods**

During the AGF process, thermal-hydraulic (TH) coupling is a key issue to be solved. As soil temperatures decrease below the freezing point, free water in the pores starts to freeze and the resulting ice crystals clog the original pore channels of groundwater flow. Simultaneously, due to the existence of convection, the change of the groundwater flow field will affect heat transfer in soil in return. In particular, the water-ice phase transition will lead to small changes in the physical properties of soil materials. The typical phase diagram of the freezing process is shown in Figure 1. In a saturated representative volume element (RVE), the pore space *ε* consists of ice, free water, and unfrozen film water, which are denoted by subscripts *i*, w, and *rs*, respectively. The saturation *S*, which describes the proportion of the volume occupied by the components in the pores, is introduced.

**Figure 1.** Schematic of a typical phase diagram for the freezing process.

Based on the above coupling mechanism, we construct the mathematical model of this problem by establishing the governing equations of the groundwater flow field and the freezing temperature field, respectively. Through the multi-physics finite element software COMSOL Multiphysics®, this problem can be numerically solved. Assumptions on this issue are required, as followed:


#### *2.1. Water Mass Conservation Equation*

Due to the barrier effect of ice, water mass transfer is inevitably affected by the crystallization, which is controlled by temperature. Moreover, water-ice phase transition will be accompanied by density change between liquid and solid, resulting in a mass loss. Thus, the total water mass conservation could be governed by Darcy's law, as follows:

$$(1 - \varepsilon S\_i) S\_{OP} \frac{\partial p}{\partial t} + \nabla \left[ -k\_r K \cdot (\nabla p + \rho\_w g \nabla D) \right] = Q\_S + Q\_T \tag{1}$$

where *SOP* is the specific pressure storativity; *p* is the pressure (Pa); *ε* is the porosity; *t* is the time (s); *kr* is the effective hydraulic conductivity that used to represent the decrease of permeability in the water-ice phase transition zone; *K* is the hydraulic conductivity (m/s); *ρ<sup>w</sup>* denotes the fluid density (kg/m3); ∇*<sup>D</sup>* is the gravitational potential gradient tensor; and, *QS* denotes sources and sinks. *QT* represents temperature-induced mass increments, and in this case, this term can be defined as

$$Q\_T = \varepsilon (\rho\_w - \rho\_i) \frac{\partial S\_i}{\partial t} \tag{2}$$

where *Si* is the saturation of ice in the pores; and, *ρ<sup>i</sup>* denotes the ice density (kg/m3).

#### *2.2. Energy Conversation Equation*

Due to the groundwater flow, heat transfer in the form of convection becomes non-negligible. The higher the groundwater flow velocity is, the more intense its impact on heat transfer will be. Moreover, the heat transfer, in this case, contains the water-ice phase transition process. It causes a dramatic energy shift in a narrow range of phase transition. Therefore, in addition to heat conduction, the energy conversation equation should also include heat convection and phase transition terms,

$$\mathbb{C}\_{\epsilon q} \frac{\partial T}{\partial t} - \nabla \left( \lambda\_{\epsilon q} \nabla T \right) + \mathbb{C}\_{w} \vec{u}^{\dagger} \nabla T - \rho\_{w} L \frac{\partial \mathcal{S}\_{w}}{\partial T} = Q\_{G} \tag{3}$$

where *T* is the temperature (K); <sup>→</sup> *u* is seepage velocity vector, as calculated by Darcy's law ( → *u* = −*krK* ∇*p <sup>ρ</sup>wg* + ∇*D* ); *L* is the latent heat released during phase transition (J/kg); *Cw* is the fluid heat capacity ((J/(kg·K)); *Sw* is the saturation of water in the pores; and, *∂Sw*/*∂T* represents the water volume fraction per volume of pore medium that is produced by freezing as *T* decreases. *C* is the bulk heat capacity and *λ* is the bulk thermal conductivity. The subscript *eq* denotes that the physical variables are calculated by using the volumetrically weighted arithmetical mean method, which is,

$$\mathcal{L}\_{\varepsilon\eta} = \varepsilon (\mathcal{S}\_{\overline{w}} \rho\_{\overline{w}} c\_{\overline{w}} + \mathcal{S}\_{\overline{i}} \rho\_{\overline{i}} c\_{\overline{i}}) + (1 - \varepsilon) \rho\_{\overline{s}} c\_{\overline{s}} \tag{4}$$

$$
\lambda\_{\epsilon\eta} = \varepsilon (S\_w \eta\_w + S\_i \eta\_i) + (1 - \varepsilon)\eta\_s \tag{5}
$$

where *c* is the specific heat capacity (J/(kg·K)) and *η* is the thermal conductivity (W/(m·K)). The subscript *i*, w, and *s* denote ice, water, and grain particles, respectively.

#### *2.3. TH Coupling Parameters Description*

Based on the theory of groundwater flow and heat transfer in porous media, the key coupling factors of TH governing equations are effective hydraulic conductivity and water saturation. They will have a major impact on the distribution of temperature and groundwater flow fields

#### 2.3.1. Effective Hydraulic Conductivity

With the decreasing of temperature, the ice that forms in the pore will block the original groundwater flow channel and it results in lower soil permeability. The effective hydraulic conductivity, defined as a function of temperature, is introduced to describe this permeability reduction. Related research shows that the hydraulic conductivity during freezing depends on temperature and soil type [41–43]. The experimental results from Burt and Williams (1976) [42] show that the hydraulic conductivity will drop down to 10−<sup>4</sup> times below the original value within a narrow temperature range below freezing point. Some other scholars describe this process by introducing the parameter of relative permeability [26,44]. The common functions that describe the relationship between the physical parameters and changing temperatures with phase change are piecewise linear function [28] and exponential expression [26,29]. In addition to the well-known geotechnical software GeoStudio [45], numerical software using the smoothed step function to describe parameter changes during phase transitions is less common. Due to the release of latent heat during the ice-water phase change, the nonlinearity of the convection-diffusion equation is strengthened. In this paper, we adopt smoothed Heaviside step function to characterize the relationship between permeability and temperature. The effective hydraulic conductivity is defined as

$$k\_{ef} = (1 - k\_r) f l c 2 \text{hs} \left( T - T\_{p\prime} \delta T \right) + k\_r \tag{6}$$

where *kr* is residual hydraulic conductivity, i.e. the hydraulic conductivity in the frozen body, and it is usually set very small but cannot be equal to zero; and, *flc2hs* represents the smoothed Heaviside function with a continuous second derivative without overshoot, which is a built-in step-function in COMSOL. The function *flc*2*hs*(*x*,*d*) is equal to 0 for *x* < −*d*, and to 1 for *x* > *d*. In between (−*d* < *x* < *d*), the function goes from 0 to 1, and at *x* = 0, *flc2hs* = 0.5 (more details refer to [29]); *Tp* is the freezing point; and, *δT* is the radius of the phase transition interval.

#### 2.3.2. The Function of Water Saturation and Freezing

The function of water saturation describes the changing water content in the pores during the phase change. Obviously, it is affected by temperature and therefore a function of temperature. The saturated soils in the unfrozen zone are filled with water in their pores. In this case, the water saturation is 1. As the soil temperature drops below the freezing point of pure water, phase change of the pore water will take place, resulting in a gradual decrease in water saturation. Due to the presence

of a certain amount of dissolved salt in the groundwater, the purification effect of ice causes a gradual increase of pore water concentration in the local area. In addition, some of the water exists in the form of a thin film. Hence, when more water turns into ice, the capillary and the absorption force will be enlarged, which makes it harder for the water to freeze [46]. As a result, some super-cooled unfrozen water remains in the frozen zone (*Srw*). Besides the descriptions with linear function, some scholars use exponential functions to describe the residual saturation [15,34]. In order to facilitate the numerical model, the water saturation function is described by the following Heaviside function,

$$S\_{w} = (1 - S\_{rs}) f l c 2 \text{hs} \left( T - T\_{p\_{\prime}} \delta T \right) + S\_{rw} \tag{7}$$

The freezing function is used to determine the amount of released latent heat in Equation (3) during the phase change, as −*∂Sw*/*∂T*.

#### *2.4. Numerical Implementation*

The numerical implementation of the mathematical model is performed through the finite element numerical software COMSOL Multiphysics®. This software can provide an interface between the interconnected variables of different physical fields. In the study of Scheidegger et al. (2017) [29], the TH coupling permafrost model was also established with this software. Compared with the permafrost model in their study, our AGF model faces a more complicated TH coupling process (mainly reflected in the interactional influence of freezing front shape and seepage velocity). In order to characterize this process, the numerical calculation requires much more detailed grids, especially at the vicinity of the freezing front. In our numerical study, we therefore utilized the self-adapting gridding function to initialize the detailed finite elements in order to improve the convergence of the model. Two smoothed step functions are used to describe the change of effective thermal conductivity and unfrozen water saturation with changing temperature, respectively, with the expectation that the non-linearity at the freezing front could be decreased. The feasibility of the new functions will be verified through the subsequent analytical solution and the case study of the engineering project. Previous experience has shown the obvious advantages of this software in solving complex models of coupled multi-physics. Therefore, the numerical solution of the artificial ground freezing method combined with groundwater flow in this paper is implemented in a COMSOL model, which is based on a full coupling of temperature and flow fields by combining the physical interfaces of Darcy's Law and Heat Transfer in Porous Media. The effective hydraulic conductivity (*kr*) and the function of water saturation (*Sw*) and its derivative with respect to *T* (*∂Sw*/*∂T*) in this model are defined with the smoothed step function. The relationships between them and various temperatures are shown in Figure 2.

**Figure 2.** The temperature dependent parameters employed in the numerical implementation: (**a**) Effective hydraulic conductivity curve with various temperatures; and, (**b**) The function of water saturation and freezing with various temperatures.

#### **3. Validation and Simulation Studies**

#### *3.1. Comparison with Analytical Solution Based Results*

The established numerical model for the heat transfer with phase change is verified by comparing the numerically calculated temperature with the results that are based on the classical analytical solution of freezing temperature field proposed by Lunardini (1985) [27] (hereafter referred to as the "Lunardini solution"). The Lunardini solution is firstly developed to solve a one-dimensional, semi-infinite soil thawing problem. It accommodates heat advection via water flow. In the study of Kurylyk et al. (2014) [28], the authors compared their numerical results with this analytical solution. In this study, we also use this analytical solution to validate our numerical model. The employed geometry, initial conditions, and boundary conditions in this study are shown in Figure 3. The detailed physical parameters of materials are given by Kurylyk et al. (2014) [28].

**Figure 3.** Comparison of the results of analytical and numerical solutions. (**a**) Model geometry, initial conditions and boundary conditions employed in COMSOL; and, (**b**) position of the thawing front versus time at three Darcy velocities: 0.001 m/year (light gray, negligible advection), 10 m/year (medium gray), and 100 m/year (black) calculated by Lunardini's analytical solution (lines) and numerical solution (squares).

To verify the performance of this model at different groundwater flow velocities, we specified three different boundary conditions at the top: 0.001 m/year, 10 m/year, 100 m/year. The comparison of the results of analytical and numerical solutions shows a good consistency in Figure 3, which indicates that this numerical model is reliable for solving the heat transfer problem with phase change and existing groundwater flow. Beyond this, this model will be further verified by the in-situ case study in Section 4.

#### *3.2. Simulation Studies*

Groundwater flow will affect the soil temperature distribution in the freezing process. During this process, the groundwater flow rate is one of the key factors. A large amount of research work has been done on the influence of the groundwater flow rate on the closure time of the freezing wall [27,37,47–49]. However, apart from the groundwater flow rate, other parameters, such as pipe spacing, thermodynamic parameters, and so on are also important factors that affect the freezing temperature field [50]. In our study through a series of numerical experiments, we simulate and analyze the effects of pipe spacing, thermal conductivity of soil particles, and groundwater flow velocity on the temperature distribution and the closure time of the freezing wall. The geometric model used in the numerical experiment is shown in Figure 4. The 16 freezing pipes are arranged in a square with the same pipe radius of r0 = 0.054 m. The ambient temperature (*T*0) is 15 ◦C and the initial hydraulic head has a difference of =0.8 m (Δ*H*) between left and right boundary. The temperature of the freezing pipes (*Tp*) is −30 ◦C. For detailed parameter setup, please refer to Table 1.

Based on the numerical model introduced above, three sets of parametric sweep experiments are designed for three factors, which are the head difference (Δ*H*), thermal conductivity of soil grain (*λs*), and pipe spacing (*l*): (1) when *l* = 1.1 and *λ<sup>s</sup>* = 2.0, Δ*H* values are 0.5, 1.0, 1.5; (2) when *l* = 1.1 and Δ*H* = 1.0, the *λ<sup>s</sup>* values are 1.5, 2.0, 2.5; and, (3) when Δ*H* = 1.5 and *λ<sup>s</sup>* = 2.0, the values of *l* are 1.0, 1.1, 1.2 (for the individual unit please refer to Table 1). The experimental results of temperature distribution are shown in Figure 5. Parametric sweep results show that as the head difference and pipe spacing increase, the formation of freezing walls becomes more difficult. On the contrary, as the thermal conductivity of soil grain increases, the thickness of the frozen wall gradually increases.

**Figure 4.** Model geometry, initial and boundary conditions for the simulation of the artificial ground freezing (AGF) project.


**Table 1.** Parameters used in numerical experiments.

Note: <sup>1</sup> Value for parametric sweep.

**Figure 5.** The temperature distribution results of parametric sweep experiments.

Except for the shape of the frozen wall, the closure time of the frozen wall is the key factor in evaluating the ground freezing project. The case studies introduced above show that the parameters of *λs*, *l*, and Δ*H* have an important influence on the closure time, among which the roles of *λ<sup>s</sup>* and *l* are explicit. The third parameter, Δ*H*, is the controlling factor of the groundwater velocity. The influence of groundwater velocity on the ground freezing process has been discussed for years. In the ground freezing project, if the groundwater flow exceeds a certain speed, which is defined as the critical velocity, the freezing wall will not be able to close. However, this definition is difficult to quantify, as it is essentially an economic problem. As shown in Figure 6, the closure time increases exponentially with increasing groundwater flow velocity, especially after the velocity has exceeded 1.5 m/day.

**Figure 6.** Closure times of freezing wall with various groundwater flow velocities and definition of the critical velocity.

In order to study further the determination and influence of critical velocity, we firstly consider the fact in our case that the freezing wall is planned to be closed within 40 days and define the allowed maximum groundwater flow velocity as the critical velocity. As shown in Figure 6, when *l* = 1.1 and *λ<sup>s</sup>* = 2.0, the critical velocity is 1.5 m/day. Subsequently, we set these two parameters within certain ranges (*l* ∈ [0.5, 1.5], *λ<sup>s</sup>* ∈ [1.0, 3.0]) and randomly sample 500 pairs of data to calculate numerically the corresponding critical velocity values. With respect to this critical value, linear regression analyses of pipe spacing and thermal conductivity are performed, and the results are shown in Figure 7. The slope of the regression line indicates that the pipe spacing is negatively correlated

to the critical velocity and has more influence on the critical velocity than the thermal conductivity. As shown in Figure 7a, the degree of dispersion with a lower pipe spacing value (smaller than 0.8 m) has an increasing trend, which indicates that the pipe spacing value is no more linearly correlated to the critical velocity. On the contrary, as shown in Figure 7b, the sampling points become more scattered with increasing thermal conductivity. This indicates that, when the thermal conductivity value becomes larger, the critical velocity becomes less sensitive to it. Based on both figures, the critical velocity can be mainly found within the range of 0.25 m/day~1 m/day. If the groundwater velocity does not exceed 1 m/day, in most (about 90.6%) cases the freezing wall should be able to close within 40 days. Once exceeded, extra measures (e.g., lowering the freezing temperature, decreasing the pipe spacing, etc.) should be taken to ensure the closure of the freezing wall.

**Figure 7.** Scatter plot of changes in critical groundwater flow velocity with respect to (**a**) Pipe spacing; and, (**b**) Soil grain thermal conductivity based on 500 random sampling pairs. The red line represents the linear regression analysis.

#### **4. Engineering Application**

#### *4.1. Engineering Background and the 3D Model*

The application case of this study is an AGF strengthening project for the jointing section of a tunnel and a subway station located in Southern China. According to the geological survey, a saturated aquifer is at the depth of the frozen zone with silty fine sand. According to the in-situ groundwater flow test, the horizontal groundwater flow (the arrows shown in the Figure 8b) is from the left to the right side almost perpendicular to the crosscut of the tunnel. In this AGF project, 35 frozen pipes are arranged with the diameter of 108 mm each, the designed pipe spacing 1.1 m, and the length of ca. 7 m each. According to the in-situ measurements, the arrangement of freezing pipes and monitoring the hole is shown in Figure 8. The 3D model has a geometry of 20 m × 20 m × 15 m, in which the excavated tunnel is 8 m × 7.5 m × 5 m.

**Figure 8.** The arrangement of freezing pipes and monitoring points in three-dimensional (3D) model. (**a**) Side view with blue lines representing the freezing pipes; (**b**) Front view (the arrows represent groundwater flow direction).

#### *4.2. Boundary and Intial Conditions*

According to the prior in-situ temperature monitoring data at the desired freezing depth, the initial ground temperature is about 15 ◦C. The box chart of initial temperatures in different thermo-observation holes is shown in Figure 9.

**Figure 9.** The box chart of initial ground temperatures in different thermo-observation holes.

The lateral wall of the freezing pipe is the cooling source of the freezing system. The change in lateral wall temperature has a great effect on the temperature distribution in the whole system. The temperature monitoring values of the main-pipe, through which the refrigerant can flow to the freezing wall, can be used to represent the lateral wall temperature. According to the measured temperature of the loop main-pipe (*Tout*) during the freezing period (50 days), the lateral wall temperature change curve is shown in Figure 10.

**Figure 10.** The measured temperature of the loop main-pipe during freezing period.

The groundwater velocity obtained through the in-situ field test is 0.5 m/day. Based on Darcy's law, the head difference between up and downstream boundary is calculated as 0.8 m. The parameters used for the engineering example in this model is shown in Table 2.


**Table 2.** Parameters used for engineering example.

#### *4.3. Simulation Results, Verification and Discussion*

The simulation results of the temperature distribution and freezing body development at various times (10 days, 20 days, 30 days, 40 days, and 50 days) are shown in Figure 11. The freezing body and streamlines development show a clear coupling effect of groundwater flow field and the freezing temperature field. After 50 days, the freezing wall has been completely closed. But, the thickness of the frozen wall at the upstream side (left side) is relatively thin, especially at the junction of this side to the top of the frozen wall (top left corner). The simulated temperature of the freezing wall at the section of *Y* = 3.5 m shows obvious non-uniform distribution. Due to the arrangement of freezing pipes and the groundwater flow, the temperatures at the top and bottom of the frozen wall are lower than those at both sides of the up- and downstream. The actual pipe spacing has an important influence on the temperature field and the closure time of the freezing wall. For example, at the upper left corner of the freezing circle, the decrease of temperature is slower, resulting in a slower closing process.

**Figure 11.** The development of freezing body with streamlines (black lines) (**left**) and temperature distribution with velocity vector (black arrows) at *Y* = 3.5 m (**right**) at various times.

According to the study in Section 3.2, the freezing wall should be able to close within 40 days if the groundwater velocity is 0.5 m/day, the soil thermal conductivity is 2.14 W/(m·K), and the pipe spacing is around 1 m. However, Figure 11 shows a closure time of 50 days. On the 40th day the freezing wall was not yet completely closed. This is due to the fact that the freezing pipes are not uniformly distributed and the pipe spacing values vary. For example, at the upper left corner the pipe spacing is 1.5 m, while at the bottom, the smallest distance between the pipes is only 0.5 m. As the frozen wall was built on the other side, and therefore the groundwater flow was limited for the left side, the unfrozen area went frozen in the next 10 days after the 40th day.

The comparison between the in-situ measured and the numerically calculated temperatures at the five monitoring points (positions shown in Figure 8) is shown in Figure 12. After the 50-day freezing period, the temperature changes of the top three monitoring points (M1, M2, M3) are relatively gentle and among which the middle measuring point M2 has a lower temperature than M1 and M3. Both monitoring points at the upstream side (M4) and the downstream side (M5) have reached below 0 ◦C after 50 days of freezing. The temperature of the downstream monitoring point decreased more rapidly and reached 0 ◦C earlier than the upstream monitoring point. Overall, the calculated curves of the temperatures at the monitored points coincide with the ones that are measured in-situ.

**Figure 12.** The temperature comparison curve of monitoring points between simulated and measured values.

The average temperature of the frozen body is one of the most important indices to evaluate the strength of the frozen wall. According to the calculated results, the curve of the average temperature of the frozen wall is plotted against time in Figure 13. The results show that the average temperature change of the frozen wall coincides with that of the main frozen pipes (Figure 10) with the same trend. After a 50-day freezing period, the average temperature of the frozen body reached below −10 ◦C.

**Figure 13.** Average temperature curve of the frozen body versus freezing time.

#### **5. Conclusions**

In this study, based on the theories of heat transport and groundwater flow in a porous medium, we set up a fully coupled numerical model with finite element software COMSOL Multiphysics®. The coupled process of heat transport and groundwater flow with phase change was simulated. Through the model validation with an analytical solution, the sensitivity analysis, as well as the simulation of a field example with in-situ measurement data, we have drawn the following conclusions.

During the phase change, the effective hydraulic conductivity and the saturation of unfrozen water are functions of the temperature that controls the heat transport process. This function is generally described with a piecewise linear function. In this paper, we decided to use the smoothed step function. The comparison of numerically calculated and analytically solved results, as well as the validation of the field model have shown that it is feasible to use the smoothed step function to describe the change of effective hydraulic conductivity *kr* and the saturation of unfrozen water *Sw* during the phase change.

With numerical case studies, sensitivity analysis was performed for three parameters: head difference (·*H*), thermal conductivity of soil grain (*λs*), and pipe spacing (*l*). The corresponding influence of these parameters on the critical groundwater flow velocity was analyzed. The results show that when the groundwater flow velocity is within a range of 0.25 m/day~1.0 m/day, in about 90% of cases the freezing wall can be closed within 40 days. Extra measures are recommended to ensure the closure of the freezing wall, if the groundwater flow velocity exceeds this range.

The 3D model for the field example with further validation has shown that, along the groundwater flow direction, the freezing wall forms slowly and on the upstream side the thickness of the freezing wall is thinner than on the downstream side. The actual pipe spacing has an important influence on the temperature field and closure time of the freezing wall. The larger the actual pipe spacing is, the slower the closing process will be. Besides this, the average temperature change of the freezing body coincides with that of the main frozen pipes with the same trend.

**Author Contributions:** R.H. and Q.L. performed the studies and wrote the paper. Y.X. has contributed to the numerical modelling.

**Acknowledgments:** The work is supported by the Ministry of Education of the People's Republic of China through the program "Research on Mechanism of Groundwater Exploitation and Seawater Intrusion in Coastal Areas" (Project Code 20165037412) and by "the Fundamental Research Funds for the Central Universities" ("Research on the hydraulic tomographical method for aquifer characterization ", Project Code 2015B29314). It is also supported by Jiangsu Provincial Department of Education, Project Code 2016B1203503.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
