Modeling Seepage Flow and Spatial Variability of Soil Thermal Conductivity during Artiﬁcial Ground Freezing for Tunnel Excavation

: Artiﬁcial ground freezing (AGF) technology has been commonly applied in tunnel construction. Its primary goal is to create a frozen wall around the tunnel proﬁle as a hydraulic barrier and temporary support, but it is inevitably affected by two natural factors. Firstly, seepage ﬂows provide large and continuous heat energy to prevent the soil from freezing. Secondly, as a key soil parameter in heat transfer, the soil thermal conductivity shows inherent spatial variability, binging uncertainties in freezing effects and efﬁciency. However, few studies have explored the inﬂuence of spatially varied soil thermal conductivity on AGF. In this study, a coupled hydro-thermal numerical model was developed to examine the effects of seepage on the formation of frozen wall. The soil thermal conductivity is simulated as a lognormal random ﬁeld and analyzed by groups of Monte-Carlo simulations. The results conﬁrmed the adverse effect of groundwater ﬂow on the formation of frozen wall, including the uneven development of frozen body towards the downstream side and the higher risk of water leakage on the upstream face of the tunnel. Based on random ﬁnite element analysis, this study quantitively tabulated the required additional freezing time above the deterministic scenario. Two levels of the additional freezing time are provided, namely the average level and conservative level, which aim to facilitate practitioners in making a rule-of-thumb estimation in the design of comparable situations. The ﬁndings can offer practitioners a rule of thumb for estimating the additional freezing times needed in artiﬁcial ground freezing, accounting for the seepage ﬂow and spatial variation in soil thermal conductivity. place seepage The temperature development and formation time of frozen wall Monte-Carlo simulations are compared to those of deterministic scenarios. The results indicate that the spatial variation of thermal conductivity


Introduction
As a temporary ground improvement and waterproof technique, artificial ground freezing (AGF) has been widely applied in tunneling or shaft engineering [1,2]. As the fluid refrigerant flows in the freeze pipes, it extracts heat from the surrounding ground and freezes pore water gradually into ice, hence increasing the strength and stability of soil [3]. The formed frozen wall can not only stabilize the ground structure, but also isolate the connection between groundwater and underground engineering, preventing potential water leakage accidents [4]. For AGF originating from natural ground freezing, many researches on natural frozen soil have been reported [5][6][7][8]. Compared with the natural soil freezing, the freezing process and effects of artificial soil freezing are more complicated.
The AGF technique is vulnerable to many factors, and previous studies have explored this issue. Vitel et al. [3] investigated the freezing efficiency of different refrigerants. Seven different refrigerants were compared and the chloride was found to be the most efficient coolant under the conditions of their study. Liu et al. [9,10] found that uncertain drilling inaccuracy or pipe inclination of freeze pipes had considerable influences on freezing effects and an additional freezing time was required to be longer with the increase of an inaccuracy level. Alzoubi et al. [11] designed a controlled laboratory scale experimental rig to mimic the AGF process and the results indicated that the flow rate of the refrigerant had a high effect on the closure time of frozen wall. Alzoubi et al. [12] reported that the spacing between freeze pipes could significantly impact the freezing time and shape of the frozen body. These factors are expected to be more carefully considered or even solved with corresponding measures and with a greater maturity of technology. However, two natural and inevitable factors, namely, groundwater seepage and uncertain soil thermal properties, are likely to affect the freezing evolution and final freezing effects of the AGF technique.
Firstly, groundwater flow with a high velocity can provide steady and continuous heat to delay or even prevent the closure of a frozen wall. Zhou and Meschkeet [13] investigated the impact of groundwater flow on the formation of a frozen wall and found that the connection between adjacent frozen piles would be delayed when the seepage flow is high. Vitel et al. [14] carried out three-dimensional ground freezing simulations under high seepage flow velocity conditions, and the progression of the freezing front was found to be asymmetrical, which extended more in the flow direction. Both Marwan et al. [15] and Huang et al. [16] suggested that seepage would affect the position and the shape of the formed frozen wall and two different optimization methods were presented to get a better arrangement of the freeze pipes.
Secondly, the mechanical properties of the frozen soil show a strong dependence of temperature and its strength increases with the decrease in temperature [7,8]. Therefore, the thermal properties of soil should be paid more attention to in the theoretical research and construction design of the AGF technique. In this regard, Papakonstantinou et al. [17] performed a numerical study for the AGF technique applied in a Naples underground project, and they found that the soil thermal conductivity was a key parameter in ground freezing. The frozen zone developed more slowly in areas with a low thermal conductivity. Semin et al. [18] pointed out that the earth heat infiltration in the AGF process showed a remarkable dependence on the thermal properties of rock mass, such as the thermal conductivity and heat capacity. However, it is known that the physical and thermal parameters of soil show inherent and great spatial variations [19][20][21][22][23], which directly influence the engineering performance, probabilistic design and reliability analysis of geotechnical structures [24][25][26][27][28]. Therefore, it is necessary to investigate whether and how the spatial variations of soil thermal properties affect the AGF technique. Wang et al. [29] investigated the impact of uncertain thermal parameters on the heat transfer around a single freeze pipe, and the soil thermal conductivity was found to be the main factor affecting thermal stability. Liu et al. [10] explored the effect of uncertain soil thermal conductivity on two freeze pipes and the results indicated that the spatially varied soil thermal conductivity had some influences on the freezing performance. However, few similar researches discussed the effects of uncertain soil thermal conductivity on a specific AGF project subject to seepage flow. No quantitative compensation remedy has been proposed for the spatial variability of soil thermal conductivity.
For a general tunnel excavation project using AGF technology, Figure 1 depicts the adverse effects of seepage flow and spatial variability of soil thermal conductivity. Distinct to the soil freezing process in homogeneous soil without seepage, groundwater flow affects the shape of the frozen body and its growth on the upstream surface of the tunnel, leading to potential leakage accidents during excavation. In addition, the spatial variability of soil thermal conductivity brings the uncertainty of freezing efficiency around each freeze pipe, increasing the adverse effects caused by seepage flow. To ensure the anti-seepage and supporting capacity of the frozen wall, more additional freezing time is required. Different to previous work, this study combines the effects of seepage flow and uncertain soil thermal conductivity on the soil freezing evolution of a metro tunnel project. In this work, the finite element software COMSOL Multiphysics ® (Version 5.3) is employed to develop a couple hydro-thermal model and simulate the temperature evolution of the AGF system. The numerical simulation is based on a full coupling of temperature and seepage field by involving the heat transfer in porous media and the physical interfaces of Darcy's law. The spectral representation method (SRM) is used to generate lognormal random fields of soil thermal conductivity, and groups of Monte-Carlo simulations are carried out for statistical analysis.

Basic Assumptions
To investigate the influences of groundwater seepage and uncertain soil thermal properties in an AGF project, this study developed a coupled hydro-thermal numerical model to simulate the AGF process, which considers the heat transfer and seepage flow. For simplicity, this study is based on the following assumptions: The AGF is simulated within a fully saturated porous medium, which comprises a soil matrix containing pores filled with water or crystal ice phases; 2.
Only heat conduction is considered, and convection and radiation are negligible; 3.
Only the phase change of water is considered instead of the water heave and water migration in soils [30,31]; 4.
The mechanical aspect of soil freezing (i.e., stress and strain) is not considered; 5.
The cooling plan is directly loaded on to the outer face of the freeze pipes, and the heat exchange between the freeze pipes and the refrigerant is ignored.

Energy Conservation
Tan et al. [32] proposed the governing equations for a fully coupled hydro-thermal model deduced from energy conservation, which considered the effects of latent heat, thermal conductivity and seepage velocity on temperature distribution. The temperature governing equation in freezing porous media is expressed as: where C eq is equivalent to the volumetric thermal capacity; T is temperature; t is freezing time; v w is the velocity field of fluid water. ρ and C are density and specific heat capacity, respectively. The solid matrix, liquid water and crystal ice are subscripted by s, w and i, respectively. q is conductive heat flux that can be written based on the Fourier law: where λ e is the effective thermal conductivity of porous ground structure, which mainly depends on the porosity of the ground, thermal conductivities and spatial distribution of each phase (soil, ice and unfrozen water). The volumetric average approximation method is usually used to calculate the value of effective thermal conductivity in the porous media [11,12,[33][34][35]. It computes the effective thermal conductivity as the weighted arithmetic mean of fluid and porous matrix conductivities: where n is porosity; λ represents the thermal conductivity; γ is the unfrozen water content within the pore fluid, its value decreases from 1 to 0 as the freezing development, regarded as a function of temperature as [6]: where w and m are material constants related to pore size distribution in soil; T 0 is the freezing point. Compared with the unfrozen zone; water/ice phase transition and release of latent heat take place in the frozen zone, so the equivalent volumetric thermal capacity; C eq , is recorded by Tan et al. [32] as below: where L f is the heat latent. C f and C u are the volumetric thermal capacities in the frozen zone and unfrozen zone, respectively. They can be expressed as the sum of volumetric heat capacities of three phases divided by their volume fractions in the frozen or unfrozen zone. To simplify Equation (5), C eq can be expressed as [16]:

Continuity Equations
According to the law of mass conservation, the differential governing equation for water flow has been derived by Tan et al. [32] as below: In this study, the flow is driven by gradients in the hydraulic field, and we can assume that the groundwater velocity field is governed by Darcy's law: where K r is relative permeability; κ is intrinsic permeability; g is gravity acceleration; P is fluid pressure; µ w is dynamic viscosity of fluid, which is a function of freezing temperature [36]: The relative permeability, K r , varies between 0 and 1 depending on the degree of freezing and unfrozen water content, γ. The widely used formula is defined by: where z is a material constant. In this work, the relations of γ-T, and K r -γ during freezing are depicted in Figure 2.

Model Description
Vitel et al. [3,37] proposed that compared with the axial heat flux between ground strata, the horizontal heat flux between freeze pipes and ground structure is more dominant. Therefore, the AGF model in this study simplified the shield tunnel of the Nanjing Metro Line 10 from a general three-dimensional model to a two-dimensional geometry [38]. As shown in Figure 3a, this tunnel, with a radius of 5.82 m, rests in a soil layer with an initial temperature of 18 • C. A stable seepage flows from the left to the right inside saturated soil. Fifty-two freeze pipes with a diameter of 0.108 m and a final cooling temperature of −30 • C are distributed on the circumference of the pipe circle, the radius of which is 6.15 m. In order to monitor the freezing process and the development of the frozen body, 17 monitoring points are set up to analyze the heat transfer process, illustrated in Figure 3b. Their specific coordinates are listed in Table 1. It can be noted that the monitoring points (except origin) form a ring with a thickness of 1.5 m, which is the minimal thickness of frozen wall to meet the demand of frozen soil strength [39]. The material parameters involved in this study are listed in Table 2.

Coordinate (m)
x (Horizontal) y (Vertical)  Table 2. Material parameters involved in the hydrothermal coupling model.

Symbol Value Unit
Density of soil ρ s 1900

Boundary Conditions
As illustrated in Figure 3a, the upper and lower boundaries are impermeable, the seepage water is driven from left to right by an added head difference between the upstream and downstream surfaces of the model, and the right boundary is set as a liquid-free surface.
The initial temperature of soil and groundwater is 18 • C, four sides of the model are treated as thermal insulation boundaries. The temperature of the outer face of each freeze pipe decreases with the time. Table 3 provides the cooling plan of freeze pipes during active freezing. The final temperature of freeze pipes is −30 • C.

Random Field of Soil Thermal Conductivity
Wang et al. [29] analyzed the stochastic thermal regime of frozen soil by taking the uncertain thermal parameters of frozen soil as random variables, stochastic processes and random fields, respectively. The results suggested that the situation of random fields was more rational. Therefore, the spatial variation of soil thermal conductivity is taken into consideration by the concept of random field in this work.
In various methods of simulating random fields, the spectral representation method (SRM) is one of the most often used methods to generate second-order stationary random fields of soil property discretization [40,41]. This research generates lognormal random fields reflecting the spatial variability of soil thermal conductivity by SRM, and groups of Monte-Carlo simulations are carried out to analyze the AGF process in such spatially variable media. The SRM usually discretizes the random field into a series of trigonometric functions with random phase angles; for this, two independent N-order random number matrices R 1 and R 2 satisfying uniform distribution are first generated. The frequencies are set to: where ω i and ω j are frequency coordinate values; ∆ω is the discrete interval on the frequency coordinate axis. The power spectral density functions, G i and G j , are real nonnegative functions of ω i and ω j , defined as below: where l x and l y are correlation distances in the xand y-directions, respectively. The standard deviations, σ i and σ j , are expressed by: Therefore, the random field f (x, y) is obtained as: Referring to Hu et al. [38]'s work, the soil thermal conductivity ranges from 1.87 to 3.40 W/(m·K). For this reason, the mean value of λ s is selected as 2.9 W/(m·K) and the coefficient of variation (COV) as 0.25 in this study; the correlation distances along the horizontal and vertical directions are set to 20 m and 5 m, respectively. Figure 4a illustrates a typical realization of soil thermal conductivity random field. Figure 4b shows the mesh size of established finite element model and temperature field at 10 days when the seepage velocity is set as 0.5 m/d.

Model Validation
Prior to numerical simulations, the developed coupled hydro-thermal model  In order to monitor the temperature development, 15 thermistors are distributed along the measuring lines (ML1 and ML2) at the middle plane (z = 0.6 m), which is the selected plane for validation study in this work. The specific locations of thermistors and material parameters involved can refer to [42]. The seepage velocity of 1.5 m/d is selected, and the numerical analysis is conducted using the developed coupled hydro-thermal model. Figure 5b shows the mesh size of the validation model and temperature field at 40 h. The red dots denote where the thermistors are placed.
In Figure 6, the simulation results of the validation model are compared with the experimental data (sourced from [42]) at five typical freezing durations (i.e., 0 h, 1 h, 5 h, 20 h, 40 h). It can be observed that the simulation results are basically in good agreement with the experimental data, except at the boundary of the laboratory model (see Figure 6c, when y = 1.0 m). This is because it is difficult to achieve completely heat insulated boundary conditions in the lab, leading to overestimated results of the freezing effect in numerical calculations. Table 4 summarizes the computational errors of numerical analysis to experimental tests at different freezing durations. The total error measure E tot considers the errors of all temperature monitoring points along the ML1 and ML2 at each moment. It is defined as: where T s and T e are the temperature values of simulated and experimental results, respectively. The comparison results show that the calculation error along the ML1 is within 19.85%, which serves to validate the developed numerical model.  The validation study is mainly used to verify the effectiveness of developed coupled hydro-thermal model in simulating the AGF process. To investigate the effects of seepage flow and spatial variability of soil thermal conductivity, a series of simulation groups are established. Firstly, four reference cases when the soil thermal conductivity keeps constant (2.9 W/(m·K)) with different seepage velocities (i.e., 0 m/d, 0.5 m/d, 1 m/d, and 1.5 m/d) are set up. Then, four groups of Monte-Carlo simulations (µ λs = 2.9 W/(m·K) and COV = 0.25) with the same seepage velocities (i.e., 0 m/d, 0.5 m/d, 1 m/d, and 1.5 m/d) are respectively calculated for 100 times, in which each simulation represents one possible working field where λ s spatially varies.

Effects of Seepage
The development of temperature field of four reference cases can be visually seen in Figure 7, which depicts the −10 • C and 10 • C temperature isotherm curves at four typical freezing durations (i.e., 10 days, 30 days, 50 days and 70 days). It can be seen that the frozen zone develops obviously towards the downstream side affected by seepage flow.  Figure 9a indicates the symmetrical development of the temperature field as well. After about 24 days without seepage, the curves of two inner points drop slightly faster than that of two outer points. This is because a stable and closed frozen wall is formed at this time, and the freezing effect of freeze pipes on the internal zone is more concentrated.
As shown in Figure 8b, when the seepage velocity is 0.5 m/d, two points at the upstream direction of the corresponding freeze pipes (A and D) cool slower than points B and E. It can be explained that the water flow from upstream releases heat to hinder the freezing effect of pipes, but the passed cold water even helps to cool the downstream zone, which is clearly illustrated by Figure 9a,b. After about 15 days, the curves of the two outer points (A and E) drop slower than that of the two inner points (D and B), respectively. This is because the frozen wall is initially formed, and the groundwater flow is blocked into the tunnel, but the flow along the outface of the frozen wall, illustrated in Figure 10 (larger arrows represent higher velocities), provides continuous heat energy for the two outer points and makes them freeze slower.
When the seepage velocity increases to 1 m/d and 1.5 m/d, even if it is far away from the freeze pipes, the origin (point C) cools down in a shorter time due to the rapid flow of passed cold water. The overall trend of curves is unchanged but the discrepancy of the freezing speed between the upstream and downstream points becomes more pronounced.    A circular closed frozen wall with a thickness of 1.5 m is often desired in geotechnical practices [39], so the closed ring formed by monitoring points (except origin) is treated as the desired region for the frozen wall in this research (see Figure 3b). Considering the high energy costs of soil freezing, minimizing the time required to achieve a fully frozen body within the desired zone is of great benefit to the economy. Yamamoto and Springman [43] suggested that freezing the frozen wall to −3 • C can ensure its strength and stability. Therefore, in this study, when all the sixteen temperature monitoring points drop below −3 • C, it is considered that a desired frozen wall with enough thickness and strength has been completely formed. However, the presence of ground water flow has considerable influence on the formation of a closed and stable frozen wall. Firstly, the water flow provides continuous heat energy to the soil within the gaps between adjacent freeze pipes, which delay the closure of the frozen wall, as shown in Figure 11. With the increase in seepage velocity from 0 m/d to 1.5 m/d, the adverse effects of water flow are more considerable. It takes more time to achieve a frozen ring around the tunnel profile. Secondly, the water flow releases more heat to the upstream zone of pipes and inhibits the expansion of the frozen body in the upstream direction, which is why the temperature curves of points A and D drop slower in Figure 8. This means that the frozen body cannot form concentrically around the freeze pipes, but expands to downstream faster and larger. Similarly, in the simulation results of Marwan et al. [15], an oval-shaped zone toward the downstream side is obtained. In this research, the frozen zone with a freezing temperature lower than −3 • C is defined as the effective frozen zone. As depicted in Figure 12, the final effective frozen zones after 80 days under different velocities of seepage are made up of two red −3 • C isotherm curves, and the desired frozen region is marked with a ring for comparison. With the increase in seepage velocity, the effective frozen zone expands to downstream more clearly, and the upstream soils of freeze pipes get more difficult to cool down, and so the lowest temperatures which points A and D drop to at the end of the simulation gradually rise when the seepage velocity is faster (see Figure 8). When the seepage velocity increases to 1.5 m/d, a completed desired frozen wall cannot even be achieved, which may lead to serious water leakage accidents in construction. Therefore, when there is a groundwater flow with high seepage velocities, it is suggested that the arrangement of freeze pipes should be optimized in practical engineering and more freeze pipes should be arranged on the upstream face of the tunnel.  Figure 12 indicates that the formation of a completed frozen wall is directly related to the soil freezing at points A and D. The seepage flow delays the expansion of an effective frozen zone in the upstream direction. Another factor affecting the heat conduction is the soil thermal conductivity, λ s . Papakonstantinou et al. [17] found that a lower soil thermal conductivity could result in a slower freezing speed during the AGF process. To make sense of the importance of soil thermal conductivity on soil freezing, the deterministic analyses of soil thermal conductivity without spatial variation are established before the Monte-Carlo simulations. In this part, the value of λ s ranges from 0.9 to 4.9 (unit: W/(m·K)) and keeps constant in each single simulation. The corresponding temperature curves for points A and D are plotted in Figure 13. A similar conclusion can be observed in that the freezing efficiency shows a positive correlation with soil thermal conductivity. Especially when 0.5 m/d seepage flow exists, the cost time cooling point A to −3 • C even increases from 20 days (λ s = 4.9 W/(m·K)) to 60 days (λ s = 0.9 W/(m·K)) (see Figure 13c). This huge discrepancy in freezing efficiency caused by varying soil thermal conductivity should not be ignored. As for the Monte-Carlo simulation groups (µ λs = 2.9 W/(m·K) and COV = 0.25), 100 times are respectively calculated with each velocity (i.e., 0 m/d, 0.5 m/d, 1 m/d, and 1.5 m/d). Figure 14 summarizes the temperature curves of five horizontal monitoring points, the black curves represent the deterministic results in reference cases where λ s is constant at 2.9 W/(m·K), and the red curves represent 100 repetitive Monte-Carlo calculations. The average, standard deviation, and range of each group's temperature at the end of the simulation (80 days) are extracted in Table 5. It is obvious that the temperature curves of points A and D fluctuate in a greater way than other temperature monitoring points, and the fluctuation range gradually increases with the acceleration of seepage velocity.   For points A and D, Figure 15 depicts the relationship between temperature and seepage velocity at three typical freezing periods (i.e., 20 days, 40 days, and 60 days). The effects of spatially varied soil thermal conductivity can be reflected by comparing the deterministic result with the mean value of the random field result. In addition, considering the high operating cost of artificial freezing, the unfavorable scenario with the worst freezing efficiency should be given priority by engineers. Therefore, the maximal value of temperature among 100 repetitive simulations is also plotted in Figure 15. It can be observed that the adverse effect of seepage flow on point A is more evident, when the seepage velocity exceeds 1 m/d, soil at point A cannot descend to the target temperature (−3 • C) at 40 days. If the seepage velocity is controlled within 1 m/d, the desired frozen wall should be formed within 60 days in most cases. One of the most concerned points for engineers is the formation time of the desired frozen wall, which is a key factor relating to the project schedule and cost in practical engineering. In reference groups, when the seepage velocity is 0 m/d, 0.5 m/d and 1 m/d, the formation time of the desired frozen wall is computed as 22 days, 29 days and 49 days, respectively. However, when the seepage velocity increases to 1.5 m/d, the desired frozen wall has not been formed, even at the end of the simulation. Similar observations were reported in the literature (e.g., [14,15]). As for the Monte-Carlo simulation groups, the formation times of the desired frozen wall are summarized in the frequency histograms and scatter plots in Figures 16 and 17. The mean formation time is calculated as 23.3 days, 32.9 days and 55.8 days respectively, marked by black dashed lines in Figure 16. Compared with the deterministic results represented by red dashed lines, the average time of forming the desired frozen wall in Monte-Carlo groups is obviously longer. With the increase in seepage velocity, the difference between them becomes greater. When the seepage velocity reaches 1 m/d, 7% of the simulations fail to achieve the desired frozen wall even after 80 days, so 1 m/d is roughly considered as a critical reference velocity by writers, and additional freezing time and anti-seepage measures are needed when applying AGF technique in seepage fields higher than 1 m/d.  In addition, Figure 17 clearly demonstrates that when no seepage flow occurs, the formation time of the frozen wall is relatively stable around the mean value, which indicates that the spatial variation of soil thermal conductivity shows little effect. However, the increase in seepage velocity leads to obvious fluctuations of the formation time. It is reasonable to believe that the spatial variation of soil thermal conductivity has a remarkable influence on soil freezing when high-speed seepage exists. To ensure the formation of the desired frozen wall, the additional freezing time is required to be longer, with an increase in seepage velocity.

Effects of Soil Thermal Conductivity
Based on the data from simulation results, Figure 18 provides the additional freezing time above the deterministic value against the seepage velocity. The solid curve demonstrates the adverse effect caused purely by seepage flow, while the joint effect of seepage flow and spatially varied soil thermal conductivity is represented by the dashed curves. Figure 18 can be used as a rule of thumb for practitioners to estimate the additional freezing time above the deterministic value in order to account for the seepage flow and spatial variation in soil thermal conductivity. Locatelli et al. (2019) [44] pointed out that the mean velocity of groundwater flow is about 0.345 m/d, which means around 32% additional freezing time is required under an average level, and 55% under the conservative level according to Figure 18.

Conclusions
This study developed a coupled hydro-thermal numerical model to investigate the performance of an AGF system in tunnel excavations. Groundwater flow and spatial variation of soil thermal conductivity were considered as influencing factors. The concept of random field is adopted to characterize the spatial variation of soil thermal conductivity. The coupled hydro-thermal process was simulated by a finite element method incorporated with spatially random thermal conductivity. Groups of Monte-Carlo simulations have been carried out to assess the effects of soil heterogeneity on the AGF under different velocities of seepage flow. The model is validated by laboratory test results.
The simulation results shown that the AGF in tunnel excavation is vulnerable to seepage and spatial variability of soil thermal conductivity, since a small portion of leakage can undermine the closure of a frozen wall. The place most prone to seepage accidents is at the upstream face of the tunnel. The temperature development and formation time of a desired frozen wall from Monte-Carlo simulations are compared to those of deterministic scenarios. The results indicate that the spatial variation of thermal conductivity adversely affects the freezing performance, especially for scenarios with a velocity of seepage flow higher than 1 m/d, where the adverse effect of the spatially varied soil thermal conductivity on freezing efficiency is even more pronounced.
As a compensation remedy to account for the coupled hydro-thermal effects caused by seepage flow and soil heterogeneity, additional freezing time is required to ensure the formation of a desired frozen wall. Based on random finite element analysis, this study quantitively tabulated the required additional freezing time above the deterministic scenario. Two levels of the additional freezing time are provided, namely the average level and conservative level, aiming to facilitate practitioners in making a rule of thumb estimation in the design for comparable situations.
The developed model is limited to the study of the hydro-thermal process of the AGF technique; however, the mechanical behaviors of soil freezing are ignored in this model. To obtain more practical recommendations, it is necessary to consider the fundamental mechanisms of the AGF system, as well as the spatial randomness of permeability in future work. In addition, only typical random field parameters were considered in the current study. These limitations will be considered in future work.

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