Abstract
The paper examines the theoretical issues of using borehole temperature survey data to control a frozen wall formed around the sinking mine shafts of the Nezhinsk mining and processing plant potash mine. We consider adjusting the parameters of the mathematical model of the frozen soil based on temperature measurements in boreholes. Adjustment of the parameters of the mathematical model (thermophysical properties of the soil) is usually carried out by minimizing the discrepancy functional between the experimentally measured and model temperatures in the temperature control boreholes. An important question about the form of this functional and the existence of minima remained after the previous studies. The study aimed at this question included analysis of heat transfer in two horizontal layers (sand and chalk) for two shafts under construction using artificial ground freezing. It was shown that the discrepancy functional minimum under certain conditions moves over time or is nonunique. This phenomenon results in ambiguity in adjusting the mathematical model parameters in the frozen soil to fit the borehole temperature survey data. At the stage of the frozen wall growth, the effective thermal conductivity in the frozen zone can be determined ambiguously from the temperature measurements in the boreholes—its value can change over time. At the stage of maintaining the frozen wall, the solution turns out to be dependent on the ratio of effective thermal conductivities in the frozen and unfrozen zones.
1. Introduction
Construction of mine shafts and underground tunnels in complex hydrogeological conditions requires special techniques, specifically the artificial ground freezing (AGF) technique [1–4]. The purpose of the AGF is the formation of a frozen wall (FW) of a given thickness, sufficient to resist the pressure of the surrounding unfrozen rocks and the pore water contained therein to prevent groundwater ingress into the mine shaft.
The current safety regulations in many countries enjoin organizing systematic monitoring of the frozen soil state. According to Construction Standard VSN 189-78 in Russia, monitoring the freezing process and the state of the FW should be carried out using hydrogeological and temperature control boreholes. Conducting temperature surveys in temperature control (TC) boreholes is currently the main and most informative method of analyzing the actual FW properties. This experimental method is used for the construction of mine shafts and underground tunnels, both in Russia [5–7] and in other countries [8–11].
As a rule, experimental control of the FW state is accompanied by theoretical calculations of the formation time of an FW of a given thickness [4, 12, 13]. In practice, the experimental and theoretical techniques are often used jointly: in this case, experimentally measured temperatures of the soil in TC boreholes are used to adjust the mathematical model of the frozen soil, and the latter is used to determine the temperature field in the whole soil volume subject to the thermal effect of the freezing system [14–16]. This approach is called back analysis [14] or the solution of the coefficient inverse Stefan problem [16]. Adjusting the soil’s thermal conductivity values is carried out in such a way as to find the best match between the calculated and measured temperatures in TC boreholes [17]. The thermal conductivities are usually chosen as adjustable parameters because the accuracy of their laboratory determination in core samples is usually lower than that for other parameters of the problem.
As shown in [16], there may exist more than one solution to the inverse Stefan problem in practice in some situations. Physically, this means that the soil’s thermal properties may be ambiguously determined from temperature measurements in TC boreholes. In this situation, there is a high risk of incorrect selection of the soil’s thermal properties from the obtained set of the inverse problem solutions. This result is unfavorable because, in this situation, there may appear an uncontrollable error in determining the soil temperature field at a spatial distance from the TC boreholes, as well as in the theoretical forecasting of the state of the frozen soil in the future, without relying on the TC boreholes’ measurement data.
For this reason, it is important to determine the conditions under which there appears ambiguity when defining the thermal soil properties as a result of an adjustment of the mathematical model using the measured temperatures in TC boreholes. Thus far, in the existing scientific literature, this particular issue remains unaddressed.
This paper includes the analysis of the conditions when the thermal soil properties defined from the data measured in TC boreholes become ambiguous. The analysis was made on actual temperature survey data obtained in monitoring the FW formation around the sinking vertical mine shafts of the Nezhinsk mining and processing plant in the Republic of Belarus.
2. Mathematical Model of the Frozen Soil
A brine (indirect) system of AGF was considered. A contour of vertical freezing pipes was mounted around the projected mine shaft. The frozen soil volume was divided into several horizontal layers, each of which was considered separately (see Figure 1).
Theoretical analysis of heat transfer in each soil layer during its artificial freezing was implemented by solving the generalized Stefan problem [18]. The mathematical formulation of the Stefan problem included the following assumptions for each soil layer: (1)Thermal soil properties in frozen and unfrozen zones are homogeneous and isotropic(2)Vertical heat flow is absent; the temperature distribution is uniform over the entire layer(3)Pore water migration is not considered(4)The soil is completely saturated with moisture, and its pore space contains no gas(5)There is a local heat balance between different phases (dry matrix, pore water, and pore ice)(6)Complete phase transition of pore moisture occurs at a certain finite temperature range
When specifying the locations of the freezing pipes and TC boreholes, their actual horizontal deviations from the design locations were considered. Figure 1 schematically illustrates the freezing contour without showing the actual deviations for each freezing pipe.
Given the assumptions made, the generalized 2D Stefan problem in enthalpy formulation is as follows [16, 19]: in which is the soil layer number; is the specific enthalpy of the soil (J/m3); are physical coordinates (m); is the physical time (s); are thermal conductivities of the soil in the unfrozen and frozen zones, respectively (W/(m·°C)); are the specific heat capacities of the soil in the unfrozen and frozen zones, respectively (J/(kg·°C)); are the soil densities in the unfrozen and frozen zones, respectively (kg/m3); is the liquidus temperature for pore water (°С); is the solidus temperature for pore ice (°С); is the volume fraction of ice in the soil pores (iciness) (m3/m3); is the specific heat of water crystallization (J/kg); is the soil moisture content [20] (kg/kg); is the coolant temperature (°С); is the temperature of virgin rock (°С); is the heat transfer coefficient at the freezing pipe wall (W/(m2 °С)); is the Nusselt number; is the Peclet number; is the coolant’s heat conductivity (W/(m·°C)); is the specific heat of the coolant (J/(kg·°C)); is the coolant density (kg/m3); is the radius of the freezing pipes (m); is the dimensionless factor to account for the upward coolant flow occupying only part of the cross-section of the freezing pipe; is the average coolant velocity in the freezing pipe (m/s); is the height of the freezing pipes (m); represents the surface of freezing pipes; represents the surface of the outer boundary; and is the coordinate along the normal to the surface (m).
Assumption no. 2 about the absence of vertical heat flows requires additional clarification. In fact, vertical heat flows can have a strong effect on the temperature field near the contact of soil layers with different thermophysical properties. However, if we consider a soil layer with homogeneous thermophysical properties and a sufficiently large thickness, then the temperature field in its middle horizontal section will not be affected by vertical heat transfer with adjacent soil layers for a sufficiently long time. In this sense, problems (1), (2), (3), (4), (5), (6), (7), (8), and (9) should be considered the problem of heat transfer in the middle horizontal section of a homogeneous soil layer.
The mathematical model is implemented numerically using the finite difference method in polar coordinates. The computational domain was divided into several concentric rings (see Figure 2). The boundary conditions of equal temperatures and heat fluxes were set between the adjacent rings. The densest grid was set in a ring with freezing pipes. The grid thickening near freezing pipes can be seen from Figure 2, where the freezing pipes are schematically shown using black circles. A preliminary simulation defined the size of the grid cells to ensure the independence of the solution from the gridding technique. The minimum cell size near the freezing pipes was taken to be 0.5 cm.
A Forward Time Centered Space method was used for the diffusion equation. Algorithmization was carried out in Microsoft Visual Studio 2017. A simulation was conducted for two soil layers in the freezing interval for the conditions of the sinking shafts of the Nezhinsk mining and processing plant: a layer of chalk at the depth interval of and a layer of sand at the depth interval of . The thermal properties of these layers, taken from the site survey reports of the considered area [21], are summarized in Table 1.
It can be seen from Table 1 that some thermophysical properties (especially density) of the sand layer are sufficiently different for the first and second shafts. For the chalk layer, these differences are minimal. Such differences are due to the fact that the thermophysical properties of soils for each shaft were determined using core samples taken from different exploratory boreholes.
Table 2 shows the problem geometry. Figure 3 presents the time charts of the temperature and coolant flow rate in the freezing pipes. The coolant was a 25% CaCl2 solution with thermal properties taken from [22]: , , and . The maximum coolant flow rate in the freezing system (263 m3/s) corresponds to the average coolant velocity in the freezing pipe’s outer section, equal to approximately 0.22 m/s, and the Reynolds number equal to approximately 580. In this case, the flow is laminar (). Equations (6) and (7) are used to calculate the heat transfer coefficient corresponding to the laminar flow regime.
(a)
(b)
(c)
(d)
The frozen soil state in each of the two sinking shafts was observed using four TC boreholes located as shown in Figure 1. Temperature profiles measured in TC boreholes and shown in Figure 4 were used to refine the soil’s thermal properties (thermal conductivity of frozen and unfrozen zones, moisture content) and adjust the corresponding mathematical model of heat transfer in frozen soil. Thermal properties of the soil were adjusted by minimizing the discrepancy functional of the temperature in TC boreholes measured experimentally () and calculated theoretically () by solving the direct Stefan problems (1), (2), (3), (4), (5), (6), (7), (8), and (9): in which is the number of TC boreholes; is the temperature difference, typical for the considered problem (°С); and is the modeling time (days).
(a)
(b)
The thermal soil properties obtained as a result of the functional minimizing (10) are effective thermal properties, which implicitly include other physical processes and factors that are not considered in the model: soil heterogeneity, error in the determination of actual deviations of borehole positions, heat transfer due to possible moisture migration, and others.
2.1. Simulation Results
In this work, we investigated the form of functional in the phase space of the adjustable properties of the problem ( and ). As a result of the multiparametric numerical simulation of soil freezing performed at different values of thermal conductivities and , the isolines of the functional and its minimum location in the phase space of thermal conductivities at different time points for two mine shafts and two studied soil layers (see Figures 5–7) were obtained.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Generally, in all four considered cases (2 shafts and 2 soil layers), we obtained a similar qualitative result for the functional minimum position and its form. At the initial time interval from beginning the freezing system (0–50 days), the minimum was unique, and its position gradually changed over time (shifting from right to left). In the time interval of 50–100 days, the minimum changed shape from a point to a line. The position of the line at subsequent times did not essentially change (times up to 250 days from the freezing start were studied).
On some plots (Figures 5(d) and 7(d)), instead of the line of minima, several local minima are on the same line—the result of the specifics of interpolation of the two-dimensional field of the functional based on a limited number of calculated discrete points ( grid).
Quantitative analysis of the functional isolines made it possible to acquire the following interesting features. At the initial time interval of 0–50 days, the functional minimum moved almost parallel to the -axis for both soil layers; i.e., during this time interval, thermal conductivity in the frozen zone strongly depended on time. Consequently, the determination of the thermal conductivity value in the frozen zone from the borehole temperature survey data may appear unreliable in this time interval. Furthermore, the functional minimum positions at the initial time interval for both soil layers are different in different shafts, both along the - and -axes. This is especially noticeable for short times (25 days; see Figures 5(a) and 7(a)) and is probably due to the strong influence of soil heterogeneity, deviation survey data errors, and other factors not considered in the model.
For simulation times of more than 100 days, the positions of the lines of minima of the functional for different shafts were approximately the same for both chalk and sand. The inclination of the line of minima of the functional for the chalk ranged from 75 to 84°, while the inclination of the line of minima of the functional for the sand varied from 73 to 81°. These angle values were close enough to 90° indicating that at such times, the solution of the Stefan problem is much more sensitive to thermal conductivity in the frozen zone than in the unfrozen one. This result also indicates that the solution of the inverse Stefan problem at greater simulation times does not provide for the correct definition of thermal conductivity in the frozen and unfrozen zones simultaneously.
3. Theoretical Interpretation
The obtained dependencies of the change in the functional (10) position and form were also analyzed qualitatively using a simplified one-dimensional heat balance model. For simplicity, it was assumed that the soil had a sharp phase boundary separating the frozen and unfrozen zones. In this case, the heat balance equation for some small section of such a boundary is described as follows: in which is a coordinate on the normal to the phase boundary (m), is the coordinate describing the position of the phase boundary (m), «» index indicates that the formula is written for the region to the right of the boundary (unfrozen zone), and «» index is for the region to the left of the boundary (frozen zone).
Formula (11) is the boundary condition at the phase transition front in the classic one-dimensional formulation of the Stefan problem [23]. The first summand on the left (11) characterizes the heat outflow to the freezing pipe, and the second one characterizes the heat inflow from the unfrozen soil. The difference between them forces displacement of the position of the phase transition front and FW growth (see Figure 8).
The temperature distributions to the left and the right of the phase transition front can be simplified using the one-dimensional model proposed in [12]: in which is the phase transition point (°С), is the radius of the freezing pipe (m), and is the thermal influence radius of the freezing system (m). Index «» corresponds to the frozen zone, and index « » corresponds to the unfrozen zone.
Equations (12) and (13) correspond to the single freezing pipe case. Coordinates of the phase transition front and the thermal influence radius are functions of time and other thermal properties of the soil. According to [12, 23, 24], they increase with time proportionally to .
First, we analyzed the contribution of each of the summands in equation (11). Thus, at the initial times (ice growing or active freezing stage), the FW thickness was small, and the temperature gradient was maximum, due to which the first summand in (11) was much greater than the second one. This follows both from the form of the derivative expressions (12) and (13) on the radial coordinate and from the results of numerical simulation on two-dimensional model (1) – (9) (see Figure 9(a)). In this figure, is the heat flow from the phase transition boundary to the freezing pipes and is the heat flow from the surrounding unfrozen soil to the phase transition boundary. The growth of FW occurs mainly due to the first summand on the left in (11). If we neglect the second summand on the left, then the solution to the Stefan problem will depend on the ratio. However, the solution remains indirectly dependent on the value through the thermal influence radius .
(a)
(b)
It is known that during the AGF, both a decrease and an increase in the moisture content of the soils can occur [25, 26]. This effect is not taken into account in the model we are considering. However, this effect can affect the measured temperatures in TC boreholes on a real object. And this, in turn, can affect the procedure for adjusting the parameters of the model. In our case, only two thermal conductivities are adjusted, and the moisture content does not change. For this reason, the thermal conductivity in the frozen zone will change in such a way as to best reflect the change in the moisture content of the frozen rocks on a real object.
Then, during the ice holding (passive freezing) stage ( days), when the FW growth slowed and/or the designated FW thickness was maintained, the heat flows on the left side in (11) equalized and almost compensated each other. In this case, the value slightly changes or remains unchanged. So, if one neglects the right-hand term (11) as the smallest summand in terms of the absolute value, it will be easy to see that the solution will depend on the ratio of thermal conductivities . This ratio, according to (12) and (13), is
The temperature ratio on the right side of (14) during the ice holding stage is greater than 1, and the logarithm ratio is also sufficiently greater than 1 since the frozen zone is usually much smaller than the unfrozen zone of thermally disturbed soil. The thermal conductivity ratio is significantly higher than unity, which corresponds to a large inclination of the line of minima in Figures 5–7 at times over 100 days.
A qualitative estimate of the value (14) for the considered parameters of the freezing system gives the value 3.4, wherein the functional line of minima for chalk at a time of 150 days (Figure 5(d)) can be set approximately by the following approximating function: . Thus, the coefficient near is fairly well predicted using an example on a simplified one-dimensional heat balance model. However, the approximating dependence for the line of minima also has a rather large intercept term of 6.81. This is due to the fact that in reality, the FW thicknesses for sand and chalk in Figure 9(b) grow with time during both the ice growth and ice holding stages and do not remain constant, as we assumed in this simplified model. The increase in FW thickness during the ice holding stage is due to the fact that a nonoptimal operating mode of the freezing station was selected in practice. However, during the ice holding stage, the time dependence is close to linear. For this reason, the right-hand side of (11) is constant during this period. Consequently, the conclusion drawn about the dependence of thermal conductivities remains valid for the case under consideration, but the very form of the dependence is linear. If we quantify the intercept term in the equation then we get the value which also correlates well with the intercept in the approximating function for the chalk. When carrying out a numerical estimate, it was assumed that corresponds to the half-thickness of the FW .
The temporary change of the position of the functional minimum (10) during the ice growing stage can have another possible explanation associated with the error in the data on the position of the boreholes. We conducted an additional sensitivity analysis of the solution to the temperature measurement error introduced through the error of the borehole deviation survey. Based on our experience, the borehole deviation survey error can be large and exceed 0.5 m for depths greater than 100 m. At the initial time interval of freezing, the temperature gradient of the soil near the freezing pipes is high, which can lead to significant temperature measurement errors in TC boreholes [27].
Figure 10 shows an example where the borehole deviation survey from two different contractors gave significantly different borehole trajectories. This example does not apply to the potash mine under consideration, but to another potash mine under construction. Near each borehole, the figure shows its number and its design depth.
If one assumes that the TC borehole displacement error from well deviation survey data is (m), then the temperature error defined from (12) and (13) will be as follows:
The choice of the specific function used depends on where the TC borehole resides—the frozen or unfrozen zone. In the initial period before the phase boundary reaches the considered TC borehole, function (18) should be used, followed by function (17). Function (17) decreases with time, whereas function (18) may be deemed almost constant since, according to [12], the ratio of and does not depend on time.
The phase transition front quite quickly reaches all TC boreholes, except for TC-3, located more than 2 m from the freezing contour. This results from the relationship of FW thickness vs. time for shaft no. 1 (see Figure 9(b)). Therefore, after 25–50 days from the start of freezing, formula (17) should be used. Under this formula, under conditions of the sinking shafts of the Nezhinsk mining and processing plant, the temperature error in TC boreholes induced by the inaccurate deviation survey can be approximately calculated as follows:
Here, is the distance from the TC borehole to the closest freezing pipe (1 m), and the value of the phase front radius is also 1 m. If one assumes that such a temperature error is typical for all TC boreholes, then the functional (10) error will be as follows:
This error is comparable with the values of the discrepancy functional (10), which follows from Figures 5–7, and therefore can lead to a significant change of the functional minimum position. If the actual position of the TC borehole appears closer to the freezing contour than in the mathematical model with incorrectly set displacements of boreholes, then the actual thermal effect of freezing pipes on a TC borehole will occur faster than in the model and will lead to overestimation of thermal conductivity of the soil, which happened at the initial times (as shown in Figures 5–7).
Over time, error (20) will reduce logarithmically and starting from a certain moment will stop significantly influencing the position of the functional minimum (10). Thus, the deviation survey error hypothesis also explains the temporal dynamics of the position of the minimum of functional (10) at the ice growing stage.
4. Conclusions
The paper describes the study of the discrepancy functional of experimentally measured and model temperatures in the temperature control boreholes. The functional was constructed in the phase space of the adjustable parameters of the model-thermal conductivities in frozen and unfrozen zones. It shows that under certain conditions, the discrepancy functional minimum changes with time or is nonunique. This result leads to the ambiguity of the adjustment of the parameters in the mathematical model of thermal processes in the frozen soil with temperature measurement data control.
The analysis of the ambiguity triggering events showed that at the frozen wall growth stage, the unreliable determination of thermal conductivity in the frozen zone is possible due to the substantial temperature gradient of the soil near the temperature control boreholes and deviation survey error. Another explanation for the unreliable determination of thermal conductivity in the frozen zone is associated with a possible change in the moisture content of soils during freezing. At the ice holding stage, both heat conductivities (in the frozen and the unfrozen zones) may be defined incorrectly; as in the given time interval, the solution of the direct Stefan problem depends mostly on their ratio.
The general conclusion is that the adjustment of soil thermal properties at AGF should be implemented with consideration for the time performance of the position and shape of the discrepancy functional minimum in the temperature control boreholes. Determining the soil’s thermal properties from the borehole temperature survey data for certain time points will likely be unreliable.
The article did not investigate other possible ways of unambiguous determination of soil thermal properties, for example, other methods of the inverse Stefan problem regularization (by Tikhonov [27–29]) or other methods of calculating the discrepancy functional. Also, the change in soil moisture content during freezing and its effect on the discrepancy functional are not studied. These issues are the subject of future research by the authors.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The results of Sections 1 and 2 were obtained with financial support from the Russian Science Foundation (RSF) (project no. 17-11-01204). The results of Section 3 were obtained with financial support from the Ministry of Science and Higher Education of the Perm Krai (Russian Federation) (project no. “C-26/563”).