Numerical simulation of frozen wall formation in water-saturated rock mass by solving the Darcy-Stefan problem

Artificial ground freezing (AGF) is a common technology of shaft sinking through water-bearing strata. The AGF technique is used to create a frozen wall preventing shaft flooding. An additional factor that makes shaft sinking more complicated is associated with external groundwater flows occurred due to hydrostatic pressure gradients. In this paper, we study the influence of groundwater seepage on the frozen wall formation in fluidsaturated rock mass in the framework of the two-dimensional two-phase Darcy-Stefan problem. The results of numerical simulation of the thermal and hydraulic properties of the sandstone layer at the site of Petrikov Mining and Processing Plant are presented. It has been found that the external groundwater flow has a significant effect on the growth of a frozen wall in the case when the groundwater velocity magnitude is greater than or equal to 50 mm/day. This critical seepage velocity strongly depends on how quickly the water content and rock mass permeability decrease with decreasing temperature, or on the parameters of the rock mass freezing characteristic curve and permeability versus temperature curve. The proper setting of these parameters is a sine qua non for creating adequate mathematical models of heat and mass processes in the artificially frozen water-saturated rock mass.


INTRODUCTION
haft sinking through water-bearing strata requires special methods to prevent the flow of water in the excavation. Artificial ground freezing (AGF) is the most widely used and highly efficient method of ground water control. The ground freezing process involves drilling and installing a series of relatively closely spaced pipes around the future shaft and circulating a coolant through these pipes (see Fig. 1). This method is used when other conventional methods such as dewatering, shoring and grouting are not feasible [1]. S Freezing continues until a frozen cylinder of sufficient strength to withstand the hydrostatic pressure is formed. The frozen cylinder is also known as the frozen wall or the ice wall [2,3]. When the required thickness of the frozen wall is reached, the procedure of shaft sinking begins. The time at which the frozen wall reaches the required thickness can be determined by both experimental monitoring [4,5] and mathematical modeling [6,7]. Due to the high velocity of groundwater in a porous rock mass, the correct mathematical prediction of the temperature field in the water-saturated rock mass is possible only in the framework of a coupled thermohydraulic problem, also known as a Darcy-Stefan problem [8 -10]. Literature review shows that there are many works devoted to the Darcy-Stefan problem, in which the freezing and thawing processes arising in soils and rocks in various natural and technical systems are studied. In [11,12], the heat and mass transfer processes in partially frozen porous media are considered in relation to permafrost thaw. Numerical analysis of heat and mass transfer in artificially frozen rocks during the construction of tunnels and mine shafts was carried out in [10, 13 -15]. In [13], the authors investigated the influence of coolant temperature, freeze pipe spacing and seepage temperature on the closure time and shape of the frozen wall. In [10], a thermo-hydraulic model of fluid-saturated rock mass was used to analyze how the heterogeneity of rock mass properties affects the formation of a frozen wall. In (15], a coupled thermo-hydromechanical mathematical model is proposed for studying the frozen wall formation during tunnel excavation. In [14], the unsteady salinity distribution near the phase transition zone was considered. The seepage of groundwater is often neglected in the calculation of temperature fields in artificially frozen rocks and soils [6, 16 -18]. The reason is that the effect of groundwater seepage is negligible in most practical cases when the content of groundwater is low, or the groundwater is stagnant. There is nevertheless a need to determine the right conditions for simulation of the frozen wall formation without consideration of groundwater seepage, and it is the purpose of this paper to find such conditions. To this end, we have formulated a two-dimensional two-phase Darcy-Stefan problem for the horizontal layer of water-saturated rock mass and developed an algorithm for a numerical solution of the stated problem of artificial ground freezing using the finite difference method.

MATHEMATICAL MODEL
he problem of artificial freezing of the water-saturated rock mass is considered. It is assumed that the physical processes occurred in the rock mass and having a strong impact on temperature distribution are as follows:  groundwater seepage;  diffusion and convection heat transport; T  phase transition of groundwater from liquid to solid state;  heat exchange between the rock mass and the coolant circulating through freezing pipes. We make the following assumptions in order to develop the model: 1. The rock mass is isotropic and fully water saturated. 2. The initial temperature field in the rock mass is homogeneous. 3. Heat and mass transfer in the vertical direction is negligible compared to that in the horizontal direction. 4. Filtration of groundwater occurs under the action of a horizontal pressure gradient. 5. Ground water and ice are incompressible. 6. Freezing pipes have deviations from the vertical position due to the error in determining the direction of drilling. It should be noted that the validity of the second assumption was verified by analyzing the conditions at the site of Petrikov deposit of potash salt in Belarus. The 3D numerical simulations of the frozen wall formation have shown that this assumption is valid when the simulation time of artificial freezing is less than 200 days and the thickness of the considered rock mass layer is more than 5 m [19]. The assumptions listed above make it possible to reduce the dimensions to 2D for the horizontal layer of the rock mass. The computational domain of this problem is shown in Fig. 2. The equations of mass and energy balance in the horizontal layer of the rock mass can be written as where l  is the groundwater density, kg/m 3 ; w is the unfrozen water content; n is the porosity; i  is the ice density, kg/m 3 ; l v is the vector of groundwater velocity, m/s; tot H is the total specific enthalpy of the water-saturated rock mass (including water and ice in pores), J/(kg°C); l H is the specific enthalpy of the groundwater, J/(kg°C);  is the thermal conductivity of water-saturated rock mass, W/(m°С); T is the temperature of water-saturated rock mass, °С; t is time, s. The thermal conductivity of the water-saturated rock mass depends on the unfrozen water content according to the following law: where 1  is the thermal conductivity of the frozen rock mass, W/(m°С); 2  is the thermal conductivity of the unfrozen rock mass, W/(m°С). The groundwater velocity is calculated from the Darcy's law: where r k is the water relative permeability; k is the absolute permeability of rock mass, m 2 ; l  is the dynamic viscosity of groundwater, Pas; l p is the hydrostatic pressure in the pore space, Pa. The dependence of the total specific enthalpy tot H on the temperature T of the water-saturated rock mass is expressed by the formula: where 1  is the density of the frozen rock mass, kg/m 3 ; 2  is the density of the unfrozen rock mass, kg/m 3 ; 1 c is the specific heat capacity of the frozen rock mass, J/(kg°C); 2 c is the specific heat capacity of the unfrozen rock mass, J/(kg°C); L is the specific heat of groundwater phase transition, J/kg; sc T is the temperature when the groundwater freezing begins, °C. The specific enthalpy l H can be represented as the function of temperature T : As follows from (6), when the unfrozen water content w is equal to zero, the specific enthalpy l H of water in the pores is also assumed to be zero, since there is no unfrozen water in the pores at such temperatures. The dependence of the unfrozen groundwater content w on the temperature T (or the soil freezing characteristic curve) can be written in the form: where B is the empirical parameter, which characterizes the reduction of the water content with decreasing temperature. It is assumed that the relative permeability of the water in pores is a temperature-dependent parameter. Hence we can write where M is the empirical parameter, which characterizes the reduction of the water relative permeability with decreasing temperature. The problem (1) -(8) is supplemented with boundary and initial conditions: where 0 sc T T  is the initial temperature of water-saturated rock mass, °С; 0 v is the initial value of groundwater velocity vector, m/s; is the heat exchange coefficient between the rock mass and the coolant in the freezing pipe, W/(m 2 °С); ( ) fr T t is the temperature of coolant, °С; out  is the outer boundary of the domain; fr  is the boundary of the freezing pipes.

NUMERICAL SIMULATION
numerical simulation of the problem (1) -(13) was undertaken using the finite difference method. An explicit scheme with first-order accuracy in time and second-order accuracy in space was applied. For the discretization of the system of Eqns. (1) -(13), we used a non-uniform rectangular grid. The grid sizes at the outer boundary out  and at the boundary fr  of the freezing pipes were selected by performing preliminary calculations to determine the temperature distribution independent of the grid parameters. At the boundary fr  , the heat flow was calculated using the effective heat exchange area eff S , which is formed by the faces of the cells intersecting with the freezing pipe circle (see Fig.   3). The inequality of the physical and effective heat exchange areas is cancelled out by means of a reduction factor R , which is equal to the ratio of the physical and effective heat exchange areas. The finite difference scheme is implemented in C # programming language. The calculation was carried out for the conditions at the site of Petrikov Mining and Processing Plant. The sandstone layer at the depth interval 136-146 m was considered. The necessary input data for numerical calculation is given in Tab. 1. The deviation of the position of a freezing pipe with increasing depth was obtained from the inclinometry data at the site. The maximum deviation of the freezing pipes at a depth of 146 m from the design position is 1.85 m, whereas the minimum deviation is 0.7 m (see Fig. 2). The diagram given in Fig.4   The results of the numerical solution of the Darcy-Stefan problem are presented in Fig. 5 -7. Temperature distributions are shown in Fig. 5. Fig. 6 demonstrates the distribution of the relative velocity after 60 days of freezing. Different values of the external seepage velocity are considered. The relative velocity is the ratio of the seepage velocity magnitude to the constant external seepage velocity magnitude: where r is the radius vector of the point, m; t is time, s. It is seen that the relative seepage velocity is equal to 1 at the outer boundary.
The groundwater flow has a significant impact on the formation of the frozen wall when the external seepage velocity is greater than or equal to 50 mm/day. In this sense, 50 mm / day is the critical seepage velocity for the considered thermal properties of the rock mass. The significant impact means that the time of closed-loop frozen wall formation (or frozen wall closure time) changes by more than 5% due to the presence of groundwater seepage. Qualitatively, this fact corresponds to the results obtained by the authors of the study [19] for the rock mass layer with similar thermal properties. The quantitative analysis shows that the results obtained in [19] give the critical seepage velocity 1.5 -2 times lower than that obtained in our calculations. This is due to the fact that another function ( ) r k T was used in this study (Heaviside step function). The performed numerical simulation showed that the critical seepage velocity strongly depends on the empirical parameters M and B in formulas (7) - (8). Fig. 7 shows the frozen wall closure time as a function of the parameters M and B. These results were determined for the external seepage velocity 150 mm/day. The closure of the frozen wall is assumed to occur when the water relative permeability value becomes equal or less than 0.02. As shown in Fig. 7, the dependence of the frozen wall closure time on the parameter M has a pronounced hyperbolic character. In the interval 0.2 0.3 M   , all three curves in the graph have asymptotes. This means that, for sufficiently low values of M, the closure of the frozen wall does not occur. When the parameter M increases from 0.3 to 2.0, the closure time decreases more than fourfold. For large values of M, when the water relative permeability falls to zero almost immediately after the temperature reaches the value sc T , the effect of parameter B on the frozen wall closure time is minimal.

CONCLUSIONS
n this paper, the two-dimensional two-phase Darcy-Stefan problem is formulated and applied to the problem of frozen wall growth simulation in the presence of external groundwater flow. A numerical algorithm for solving this problem is described. The results of a multiparameter numerical simulation of the frozen wall formation obtained in the study of the thermal and hydraulic properties of the sandstone layer at the site of Petrikov Mining and Processing Plant are presented. It was found that the external groundwater flow has a significant effect on the frozen wall growth when its velocity magnitude is greater than or equal to 50 mm/day. This critical seepage velocity depends on the thermal and hydraulic properties of the layer. Specifically, it strongly depends on how quickly the water content and rock mass permeability decrease with decreasing temperature, or on the parameters of the rock mass freezing characteristic curve and permeability versus temperature curve. Therefore, the development of an adequate mathematical model of heat and mass transfer in the artificially frozen water-saturated rock mass requires effective identification of these parameters on the basis of the laboratory test results.