
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12703-X
10.1016/j.heliyon.2024.e36672
e36672
Research Article
Analytical solution for surrounding rock temperature of cold-region tunnel considering phase change
Wu Wentao a
Guo Jiaqi gjq519@163.com
ab⁎
Wang Xiaochuan c
Hu Huanmeng c
Zhao Pengyu b
a School of Civil Engineering, Henan Polytechnic University, Jiaozuo, Henan, China
b School of Highway, Chang'an University, Xi'an, Shaanxi, China
c CCCC-SHEC Fourth Highway Engineering Co., Ltd, Luoyang, China
⁎ Corresponding author. School of Civil Engineering, Henan Polytechnic University, Jiaozuo, Henan, China. gjq519@163.com
22 8 2024
15 9 2024
22 8 2024
10 17 e3667229 5 2024
16 8 2024
20 8 2024
© 2024 Published by Elsevier Ltd.
2024

https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
The temperature of the surrounding rock in cold-region tunnels is crucial for antifreeze design, and the water-ice phase transition is essential to addressing the temperature field. This paper proposes a refined method that equates the latent heat of the ice-water phase transition to heat capacity and establishes a one-dimensional radial heat transfer model considering phase change. By defining an average thermal diffusivity coefficient through the concept of equal accumulated temperature, this method overcomes the limitations of classical heat transfer theory in directly solving the temperature field of three zone (unfrozen zone, freezing zone and frozen zone). Additionally, by employing the variable separation method and Fourier integral transformation method, the analytical formula for the transient temperature field considering phase change is derived. Then, the analytical solution was verified based on the field data. The results calculated using this method exhibit greater consistency with field temperature data and outperform the modified Stephan formula in determining the maximum frozen depth of the surrounding rock. Finally, the simplified form of the established analytical solution was further discussed. The research results can provide a theoretical basis for the analysis of the temperature field of the surrounding rock of the tunnel in cold regions and its antifreeze design.

Keywords

Cold-region tunnels
Surrounding rock with phase change
Porous media
Temperature field
Analytical solution
Latent heat
==== Body
pmc1 Introduction

To further improve the national highway network, promote economic development, and ensure the stability of the border areas, a large number of traffic tunnels will be built in the high-latitude areas in the northeast and high-altitude cold areas in the west of China [1],[2]. H o wever, due to the periodic occurrence of negative temperature with seasonal changes, these tunnels after completion are prone to frost damage such as frost heaving and cracking of tunnel lining, icing of tunnel pavement, which seriously affects the normal operation of these tunnels [[3], [4], [5]]. Frost damage has become the bottleneck for the construction and maintenance of tunnels and the sustainable development of traffic in high-latitude, high-altitude and cold regions. Study shows that frost damage is affected by two factors, low temperature and pore water [6],[7]. The study of temperature field in tunnels in cold regions can be a breakthrough to solve the problem of frost damage in tunnels. Meanwhile, the maximum frozen depth of surrounding rock determines the depth of drainage facilities and is also the main parameter for the design of thermal insulation layer in the antifreeze design of tunnels in cold regions. Therefore, it is of great significance to study the temperature field of tunnel surrounding rock in cold regions.

In view of the important role of tunnel surrounding rock temperature field in the prevention and control of frost damage, scholars at home and abroad have carried out numerous related research by using the methods of field test [8],[9], analytical calculation [[10], [11], [12]] and numerical simulation [[13], [14], [15]]. Field test is the most direct method to obtain the distribution of temperature field and explore the space-time evolution of temperature field. Zhang et al. [16] monitored five cross-sections of the internal environmental temperature field and two cross-sections of the surrounding rock temperature field, thereby establishing a three-dimensional temperature field mathematical model for cold region tunnels. However, long-term field test of tunnel surrounding rock temperature field is often limited by cost and time, and it is difficult to carry out systematic research under multiple working conditions. In addition, due to the limitation of the number of monitoring elements, it is impossible to obtain the full field information of the surrounding rock temperature, and also incapable to obtain the temperature information in advance, thus fails to provide a basis for the tunnel antifreeze design.

The numerical simulation method and analytical calculation are not subject to the above restrictions. In the aspect of numerical simulation analysis of the temperature field of the surrounding rock of the tunnel in cold regions, Zhang et al. [17] gave a finite element calculation equation based on the coupling of water and heat, and used the computer program to predict the freezing condition of the Kunlun Mountain Tunnel on the Qinghai-Tibet Railway. Zhan et al. [18] analyzed the frost heave deformation of the surrounding rock of the tunnel using COMSOL Multiphysics. Ma et al. [19] established a water-heat coupling model of surrounding rock considering water-ice phase change and thermal convection, and used COMSOL numerical simulation software to analyze the influence of insulation layer on the temperature field of tunnels in cold regions. Based on the principles of energy conservation and mass conservation, Li et al. [20] established a water-heat coupling model for tunnels in cold regions, simulated the water-heat changes around the tunnel through FORTRAN program and determined the optimal thickness of insulation layer. Zhu et al. [21] established a 3D numerical model for the unsteady temperature transfer between airflow and surrounding rock, and analyzed the impact of cold and hot air on the temperature of the tunnel's surrounding rock. To sum up, the numerical simulation method can better reproduce the actual scene, obtain the temperature information of the surrounding rock of the whole field, and visually evaluate the antifreeze effect of the insulation layer. However, the analysis process is complex, the calculation is time-consuming, and the requirements for the numerical analysis and software operation ability of engineering technicians are demanding, so it is difficult to be widely applied to engineering projects.

In terms of the analytical calculation of the temperature field of the surrounding rock of the tunnel in cold regions, Lai et al. [22] used the dimensionless and perturbation techniques to give the analytical solution of the temperature field of the circular cross-section tunnel in cold regions, but this method is only applicable to the case where the initial ground temperature is 0 °C. Zhang et al. [23] used the orthogonal and expansion theorem of Bessel function to derive the analytical solution of the temperature field of the surrounding rock of the tunnel in the cold region with a circular section, which takes into account the influence of the insulation layer. Xia et al. [24] obtained the explicit analytical solution of the transient temperature field of the surrounding rock of the tunnel in the cold region with thermal insulation layer by using the method of separating variables and Laplace transform. Lu et al. [25] established a new method of heat conduction analysis suitable for the Dirichlet transient boundary conditions. By expanding the boundary temperature with Fourier series, and taking advantage of the periodicity of boundary changes, an explicit analytical solution of heat conduction of multilayer composite plates was obtained. Liu et al. [26] established the analytical calculation model of the three-dimensional temperature field of the tunnel in cold regions, and gave the solution method through Laplace integral transformation and Fourier integral transformation. Gao et al. [27] built a convection-conduction model under the effect of train-induced wind and proposed the analytical solution of the temperature field of the tunnel in the cold area under the action of train wind. The above studies on analytical solution provide theoretical basis for the design of tunnel engineering in cold regions. However, as an important factor of the surrounding rock temperature field of cold-region tunnel, the water-ice phase change is often ignored in the analytical calculation of the surrounding rock temperature field due to the difficulty of nonlinear solution of the phase change latent heat.

In order to accurately describe the water-ice process of tunnels in cold regions and calculate its temperature field variation, it is necessary to consider the water-ice phase change in the analytical solution. In view of the difficulty in analytical calculation of the temperature field of surrounding rock with phase change, this paper, based on the apparent heat capacity method, intends to derive the transient heat transfer equation with phase change, and establishes the heat transfer model of surrounding rock of tunnel in cold regions. Firstly, the average thermal diffusion coefficient is defined based on the idea of equal accumulated temperature, and it overcomes the problem that the classical heat transfer theory cannot directly solve the continuous temperature field of “unfrozen-freezing-frozen” state. Secondly, the solution method of the transient surrounding rock temperature field of the tunnel in cold region considering phase change is given by using the separation of variables and Fourier integral transformation method. Then, the parameter sensitivity analysis of the analytical solution and the engineering case verification are carried out. Finally, the simplified form of the analytical solution of the temperature field with phase change in the surrounding rock of the tunnel in cold regions is discussed from the perspective of engineering application. The research results can provide a theoretical basis for the analysis of the temperature field of the surrounding rock of the tunnel in cold regions and its antifreeze design.

2 Apparent heat capacity method

Due to the complexity of phase change and heat transfer during rock mass freezing, the following assumptions are proposed to simplify the mathematical model: (1) The rock mass is a porous medium material, composed of rock skeleton, ice and water; (2) The model does not consider the process of water evaporation, but only considers the process of water freezing into ice and releasing latent heat.

In the process of rock mass freezing, as the thermal parameters of porous media change with time due to the freezing of water into ice and latent heat released from ice-water phase change exists, the heat transfer of porous media considering phase change is a highly nonlinear problem, which is often solved by numerical methods. In order to obtain the analytical solution of the transient temperature field of porous media considering phase change, a heat transfer model of porous media based on the apparent heat capacity method is established in this paper. In solving the phase change and heat transfer process, the method assumes that the phase change of water and ice occurs in a small temperature range, and converts the latent heat released from phase change into the apparent heat capacity in the temperature range. In this way, the phase change problem is transformed into a unidirectional nonlinear heat conduction problem in the freezing process.

The phase change process is described according to the “unfrozen-freezing-frozen” theory proposed by Konrad [28] and Yoshisuke [29]. Instead of adding a latent heat in the energy balance equation exactly when the material reaches its phase change temperature Tpc, it is assumed that the transformation occurs in a temperature interval ΔT. Within this interval, the water is transformed from liquid phase to solid phase. The corresponding volume fractions of the phases are represented by the phase transition function. As illustrated in Fig. 1, when the temperature falls below TPc−ΔT2, water completely freezes into ice; thus, δw is 0 and δi is 1. Here, δw denotes the volume fraction of liquid water, calculated as the volume of liquid water divided by the total volume of both liquid water and solid ice. Similarly, δi represents the volume fraction of solid ice, defined as the volume of solid ice divided by the total volume of both liquid water and solid ice. Conversely, when the temperature is higher than TPc+ΔT2, ice melts into water, meaning δw is 1 and δi is 0. In the freezing zone, the water is not completely frozen, and the Heaviside second-order step function δ(T) is adopted to express the volume fractions of ice and water during freezing in this study.Fig. 1 Relationship between water-ice phase volume fraction and temperature.

Fig. 1

Suppose water-ice phase change occurs and latent heat is released in the temperature range of [TPc−ΔT2, TPc+ΔT2] in the freezing zone rather than the moment when the temperature exactly meets the requirement of phase change TPc. By step function δ(T), specific enthalpy H (enthalpy per unit mass of substance (J/kg)) can be expressed by(1) H=1δiρi+δwρw(δiρiHi+δwρwHw)

According to the definition of constant pressure heat capacity, C=∂H∂T [30], the apparent heat capacity of porous media can be derived by simplifying the expression provided in Eq. (1). The resulting equation for the apparent heat capacity is:(2) Capp=1δiρi+δwρw(δiρiCi+δwρwCw)+Lⅆ12δiρi−δwρwδwρw+δiρidT

where ρ is density (kg/m3), C is heat capacity at constant pressure (J/(kg·K)), subscripts s, w, i represent rock skeleton, water and ice, respectively. L is the heat released by water of unit mass freezing into ice (J/kg), θ is volume content, Capp is apparent heat capacity (J/(m3·K)).

According to the heat transfer theory of porous media [31], the heat transfer equation of rock mass considering phase change based on the apparent heat capacity method is obtained as below,(3) (ρC)eff∂T∂t=∇(keff∇T)

where(4) (ρC)eff=θsρsCs+(1−θs)(δiρi+δwρw)Capp

(5) keff=θsks+(1−θs)(δiki+δwkw)

Where (ρC)eff is the equivalent volume heat capacity (J/(m3·K)) and keff is the equivalent heat conductivity coefficient (W/(m·K)).

3 Theoretical calculation method

To facilitate calculation, the tunnel section is equivalent to a circle, and the calculation diagram is shown in Fig. 2. It can be seen from the figure that the surrounding rock temperature field can be simplified as a one-dimensional heat conduction issue. The surface of the tunnel surrounding rock is taken as the coordinate origin and the downward direction is positive. Eq. (3) is simplified as(6) ∂T∂t=α∂2T∂x2

Where the thermal diffusion coefficient α (m2/s) is defined as the ratio of the equivalent heat conductivity coefficient keff(Eq. (5)) to the equivalent volume heat capacity (ρC)eff (Eq. (4)):(7) α=θsks+(1−θs)(δiki+δwkw)θsρsCs+(1−θs)(δiρi+δwρw)Capp

Fig. 2 Calculation model of surrounding rock temperature field of cold-region tunnel.

Fig. 2

Eq. (6) is the transient conduction equation of surrounding rock temperature field of cold-region tunnels considering phase change and Eq. (7) is the thermal diffusion coefficient considering phase change based on the apparent heat capacity method. Eq. (6) and Eq. (7) will be solved respectively in Section 3.1.

3.1 Heat transfer equation and solution

In order to solve Eq. (6), the boundary conditions and initial conditions are needed. Therein, the upper boundary condition is the air temperature inside the tunnel. According to the measured temperature inside the tunnel, the annual average temperature change function is fitted using the sine function regression method in meteorology. The fitting formula is as follow,(8) Ta=T0+Asin(ωt+φ)

Where Ta is the temperature inside the tunnel (°C), TO is the annual average temperature of the air in the tunnel (°C), A is the temperature amplitude of the air in the tunnel (°C), ω is the period of temperature change (generally 1a as a period), and φ is the initial phase and determines the starting time.

The common setting methods of the lower boundary conditions include zero geothermal flux, constant heat flux and constant temperature at corresponding depth. In view of the large deviation of the calculation results using zero geothermal flux and the difficulty of obtaining the constant temperature at the corresponding depth before design, this study selects the constant heat flux as the lower boundary condition [32]. In order to facilitate the analytical calculation, the heat flux is converted into geothermal gradient for calculation, that is,(9) T′|x=l=m

Where T′ is the partial derivative of T to x (°C/m), l2 is set to be 10 m, m is the temperature gradient (°C/m) and generally set to be 0.03 °C/m.

The initial value of surrounding rock takes the annual average temperature of surrounding rock as the upper boundary temperature, and increases linearly with the depth of surrounding rock according to the geothermal gradient, that is,(10) T|t=0=b+mx

Where b is the annual average temperature of surrounding rock (°C).

According to Eqs. (6), (8), (9), (10), the temperature field equation of surrounding rock in cold-region tunnel can be described as follows,(11) {∂T∂t=α∂2T∂x2T|x=0=Ta(t)T′|x=l=mT|t=0=b+mx

Eq. (11) is a one-dimensional heat conduction differential equation with non-homogeneous boundary conditions. In order to obtain the analytical solution, the non-homogeneous boundary shall be firstly homogenized. Suppose(12) T(x,t)=U(x,t)+V(x,t)

Substitute upper and lower boundary conditions (8) and (9) into Eq. (12), and homogenize U(x,t) boundary conditions to obtain(13) T(x,t)=U(x,t)+Ta(t)+mx

Substitute Eq. (13) into Eq. (11) and obtain the homogeneous boundary condition Eq. (14) as below,(14) {∂U∂t=α∂2U∂x2−Ta′(t)U|x=0=0U′|x=l=0U|t=0=b

According to superposition principle, the solution of Eq. (14) with homogeneous boundary condition and non-homogeneous equation with initial value being not 0 is divided into the solutions of Eq. (15) and Eq. (16). Therein, Eq. (15) contains homogeneous boundary condition and homogeneous equation with initial value being not 0 and Eq. (16) with homogeneous boundary condition and non-homogeneous equation with initial value being 0.(15) {∂U1∂t=α∂2U1∂x2U1|x=0=0U1′|x=l=0U1|t=0=b

(16) {∂U2∂t=α∂2U2∂x2−Ta′(t)U2|x=0=0U2′|x=l=0U2|t=0=0

(17) U(x,t)=U1(x,t)+U2(x,t)

Use separation variable method, Fourier series transformation and characteristic function method to solve Eq. (15) and Eq. (16) as below,(18) U1(x,t)=∑n=0∞U1n(t)sin(n+12)πlx

(19) U2(x,t)=∑n=0∞U2n(t)sin(n+12)πlx

where(20) U1n(t)=b(n+12)πe−α((n+12)πl)2t

(21) U2n(t)=−2A(n+12)π[1+α2((n+12)π)4ω2l4][sinωt+α((n+12)π)2ωl2cosωt−α((n+12)π)2ωl2e−α((n+12)πl)2t]

The first 100 terms in Eq. (20) and Eq. (21) are calculated using MATLAB. The results are then substituted into Eq. (18) and Eq. (19), followed by substituting these intermediate results into Eq. (17). Finally, the results are substituted into Eq. (13) to obtain the distribution of temperature field of cold-region tunnel surrounding rock. Wherein, the thermal diffusion coefficient α changes with time t and involves phase change. The solution of α using the integration method is presented in next section.

3.2 Derivation of thermal diffusion coefficient considering phase change

According to the “unfrozen-freezing-frozen” theoretical model, when the rock mass temperature decreases to 0 °C, the pore water begins to freeze, and the latent heat is gradually released with the temperature decreasing. How the phase transition from liquid water to solid ice takes place is described by the phase transition function. Lu [33] and Tan [34] use a smoothed step function with a continuous second derivative (Heaviside function) to describe a continuous transition from ice to water in numerical calculations. This function provides a smooth and continuous representation of the phase transition, aligning more closely with the actual physical processes. However, the use of the Heaviside function introduces complexity in analytical computations. Therefore, this paper utilizes a linear function to describe the ice-water phase change. Fig. 3 shows two types of phase transition functions, where the freezing temperature of rock mass TPc=−0.5°C and the range of freezing zone ΔT=1°C. Simplify the water phase fraction in the freezing zone to a straight line (as shown in Fig. 3), and obtain Eq. (22) as below,(22) δw(T)={0T+11T<TPc−ΔT2TPc−ΔT2≤T≤TPc+ΔT2T>TPc+ΔT2

Fig. 3 Phase transition function.

Fig. 3

Substitute Eq. (22) into Eq. (2) and obtain the apparent heat capacity as below,(23) Capp={Ci1δiρi+δwρw(δiρiCi+δwρwCw)+Lρiρw((ρw−ρi)T−ρw)2CwT<TPc−ΔT2TPc−ΔT2≤T≤TPc+ΔT2T>TPc+ΔT2

According to Eq. (23), the apparent heat capacity is the sum of an equivalent heat capacity (1δiρi+δwρw(δiρiCi+δwρwCw)) and the distribution of latent heat (Lρiρw((ρw−ρi)T−ρw)2) in freezing zone. The apparent heat capacity was calculated by substituting values δw = T+1, δi = -T, ρi = 914 kg/m3, ρw = 1000 kg/m3, Ci = 2100 J/kg/°C, Cw = 4180 J/kg/°C and L = 334000J/kg into Eq. (23), and the results are shown in Fig. 4. It is evident from the figure that the distribution of latent heat is significantly greater than equivalent heat capacity. Therefore, the apparent heat capacity can be simplified as,(24) Capp≈{CiLρiρw((ρw−ρi)T−ρw)2CwT<TPc−ΔT2TPc−ΔT2≤T≤TPc+ΔT2T>TPc+ΔT2

Fig. 4 The distribution of latent heat and equivalent heat capacity.

Fig. 4

According to Eq. (7) and Eq. (24), the thermal diffusion coefficient α of unfrozen zone, freezing zone and frozen zone can be obtained,(25) α={θsksθsks+(1−θs)(δiki+δwkw)θsρsCp,s+(1−θs)(δiρi+δwρw)Lρiρw((ρw−ρi)T−ρw)2θsks+(1−θs)kwθsρsCs+(1−θs)ρwCwT<TPc−ΔT2TPc−ΔT2≤T≤TPc+ΔT2T>TPc+ΔT2

The thermal diffusion coefficient changes with t and contains phase change, which is a strong nonlinear problem in calculation. Therefore, it is difficult to obtain its accurate analytical solution. We adopt the idea of equal accumulated temperature in meteorology to calculate the average thermal diffusion coefficient through the integration method. The specific method is as follow.

The temperature in the tunnel in cold regions varies periodically and can be fitted as a trigonometric function (as shown in Fig. 5). From Fig. 5, the road temperature is positive within the time range of 0∼t1 and t2∼365d, so the thermal diffusion coefficient of unfrozen zone is selected to calculate the temperature field during this period. Within the time range of t1∼t2, the surrounding rock begins to freeze, and freezing zone and frozen zone coexist in the surrounding rock. According to Eq. (25), the thermal diffusion coefficient of the freezing zone is far less than that of the frozen zone. At this time, it is equivalent to an insulation layer in the surrounding rock which blocks the heat transfer from the surface to the inside of the surrounding rock. Therefore, within the time range of t1∼t2, the thermal diffusion coefficient of the freezing zone is used to calculate the temperature field.Fig. 5 Temperature variation inside the tunnel.

Fig. 5

According to ∫0365α(t)dt=αr×365, the average thermal diffusion coefficient αr of surrounding rock in a freeze-thaw cycle is as follow,(26) αr=1365((3652+365πarcsin(−T0A))θsks+(1−θs)(δiki+δwkw)θsρsCs+(1−θs)(δiρi+δwρw)Lρiρw((ρw−ρi)T−ρw)2+(3652−365πarcsin(−T0A))θsks+(1−θs)kwθsρsCs+(1−θs)ρwCw)

Substitute Eq. (26) into Eq. (13) and obtain the transient analytical solution of the surrounding rock temperature field of the cold-region tunnel considering the phase change.

4 Analysis on influential factors

According to the analytical solution expression (Eq. (13), Eq. (18), Eq. (19), Eq. (26)), the temperature of tunnel surrounding rock in cold regions is related to the thermal physical parameters of surrounding rock (thermal diffusion coefficient α and the initial value of surrounding rock temperature b) and the external environmental parameters (annual temperature amplitude A and annual average temperature T0). This section analyzes the characteristics of the surrounding rock temperature field under different parameters by using the proposed analytical calculation method.

4.1 Thermal diffusion coefficient

Suppose the initial temperature of surrounding rock is 3 °C, annual temperature amplitude is 12 °C, annual average temperature is 2 °C and the variation range of thermal diffusion coefficient is from 1 × 10−7 m2/s to 9 × 10−7 m2/s, and the corresponding temperature field can be calculated as shown in Fig. 6(a) and Fig. 7(a).Fig. 6 Temperature distribution of surrounding rock at radial depth of 1 m. (a) thermal diffusion coefficient; (b) initial temperature; (c) annual temperature amplitude; (d) annual average temperature.

Fig. 6

Fig. 7 Maximum frozen circle. (a) thermal diffusion coefficient; (b) initial temperature; (c) annual temperature amplitude; (d) annual average temperature.

Fig. 7

In Fig. 5(a), with different thermal diffusion coefficients of surrounding rock, the temperature of surrounding rock at the radial depth of 1 m varies in a sine function and its amplitude increases with the increase of thermal diffusion coefficient. When the thermal diffusion coefficient is 1 × 10−7 m2/s and 9 × 10−7 m2/s, the peak temperature of surrounding rock at the radial depth of 1 m is 6.61 °C and 10.68 °C, respectively, and the valley value is −2.25 °C and −6.52 °C, respectively. Meanwhile, with the increase of thermal diffusion coefficient, the hysteresis effect of peak and valley values is getting more obvious (phase lag). When thermal diffusion coefficient is 1 × 10−7 m2/s, the surrounding rock temperature reaches peak and valley values on the 150th and 332nd day; when it is 9 × 10−7 m2/s, the surrounding rock temperature reaches peak and valley values on the 111th and 293rd day. It shows that the temperature transfer efficiency increases with the increase of thermal diffusion coefficient. In Fig. 7(a), the range of the maximum frozen circle increases with the increase of the thermal diffusion coefficient, but the increase rate is getting smaller. When the thermal diffusion coefficient is 1 × 10−7 m2/s and 9 × 10−7 m2/s, the maximum frozen circle radius is 1.47 m and 4.33 m, respectively. The thermal diffusion coefficient represents the rate at which the temperature disturbance of rock mass is transferred to a certain point by the external low temperature. When the thermal diffusion coefficient increases, the external low temperature is more easily transferred to the depth of the rock mass, resulting in an increase in the radius of the frozen circle.

4.2 Initial temperature of surrounding rock

Suppose the thermal diffusion coefficient of surrounding rock is 5 × 10−7 m2/s, annual temperature amplitude is 12 °C, annual average temperature is 2 °C and the variation range of the initial temperature of surrounding rock is 1 °C–5 °C, and the corresponding temperature field can be calculated as shown in Figs. 6(b) and Fig. 7(b).

As shown in Fig. 6(b), with different initial temperature values, the temperature change of surrounding rock at the radial depth of 1 m is almost the same, and the maximum frozen circle radius is slightly different. The radius of the frozen circle decreases with the increase of the initial temperature value of the surrounding rock (Fig. 7(b)). When the initial temperature increases from 1 °C to 5 °C, the radius of the frozen circle of the surrounding rock decreases from 3.34 m to 3.15 m, indicating the initial temperature value of surrounding rock has little influence on the analytical calculation.

4.3 Annual temperature amplitude

Suppose the thermal diffusion coefficient of surrounding rock is 5 × 10−7 m2/s, the initial temperature of surrounding rock is 3 °C, annual average temperature is 2 °C, and the variation range of annual temperature amplitude is 8 °C–16 °C, and the corresponding temperature field can be calculated as shown in Figs. 6(c) and Fig. 7(c).

As shown in Fig. 6(c), the annual temperature amplitude mainly affects the peak and valley values of the surrounding rock temperature field, and does not affect the average temperature. When the annual temperature amplitude is 8 °C, 10 °C, 12 °C, 14 °C, 16 °C, the peak value is 7.22 °C, 8.50 °C, 9.78 °C, 11.06 °C, 12.34 °C and the valley value is −3.02 °C, −4.30 °C, −5.58 °C, −6.86 °C, −8.14 °C, respectively. Clearly, the difference between adjacent peaks and adjacent valleys is the same, i.e., 1.28 °C. Although the average temperature of surrounding rock does not change with the change of annual temperature amplitude, the radius of frozen circle increases with the increase of annual temperature amplitude (Fig. 7(c)). With the same annual average temperature, the greater the temperature amplitude, the greater the amount of cold released by the external environment, and thus the frozen depth of surrounding rock increases.

4.4 Annual average temperature

Suppose the thermal diffusion coefficient of surrounding rock is 5 × 10−7 m2/s, the initial temperature of surrounding rock is 3 °C, annual temperature amplitude is 12 °C, and the variation range of annual average temperature is 0 °C–4 °C, and the corresponding temperature field can be calculated as shown in Figs. 6(d) and Fig. 7(d).

As shown in Fig. 6(d), the variation law of surrounding rock temperature field at different times has nothing to do with the annual average temperature. With the decrease of annual average temperature, the sine function of surrounding rock temperature field shifts downward as a whole. The radius of frozen circle increases with the decrease of annual average temperature, when the annual average temperature increases from 0 °C to 4 °C, the radius of frozen circle decreases from 5.03 m to 2.09 m (Fig. 7(d)). The lower the annual average temperature, the greater the amount of cold released by the external environment, leading to the increase of the frozen depth of surrounding rock.

5 Engineering case verification and analysis

5.1 Engineering overview

The Xing'anling Tunnel of Binzhou Railway is located in the northeast of Hulunbei'er City, Inner Mongolia Autonomous Region, China (Fig. 8(a)), crossing the hinterland of the middle and low mountains in the middle and south of the Greater Xing'anling, with an elevation of 900∼1210 m above sea level. The tunnel is located in the high latitude zone of Eurasia, with long and cold winter and short and warm summer. The minimum temperature is - 46.7 °C, and the initial temperature of surrounding rock is 3 °C. The mountains are often covered with thick snow, and the frozen depth is more than 3 m. The values of relevant parameters of Xing'anling Tunnel are shown in Table 1.Fig. 8 Location of Xing'anling Tunnel and temperature sensor arrangement. (a) Location of Xing'anling Tunnel; (b) Test tunnel cross-section and temperature sensor arrangement.

Fig. 8

Table 1 Relevant parameters of Xing'anling Tunnel.

Table 1Item	Thermal conductivity coefficient k/ (W/(m⋅°C))	Constant pressure heat capacity
Cp/ (J/(kg⋅°C))	Density
ρ/ (kg/m3)	Water content
θf/％	
Surrounding rock	3.2	840	2690	15	
Water	0.608	4180	1000	
Ice	2.26	2100	914	

Zhao et al. [9] carried out on-site monitoring on the surrounding rock temperature field of Xing'anling Tunnel, and arranged 15 monitoring sections symmetrically along the longitudinal direction of the tunnel. In the radial direction of the tunnel, the monitoring points are arranged at one side of the tunnel, and 11 temperature sensors are arranged at each monitoring section, as shown in Fig. 8(b). The temperature-time of all monitoring points in the 15 cross-sections of Xing'anling Tunnel varies periodically in the form of sine function and the expression of air temperature can be referred to paper [9].

5.2 Comparative analysis of temperature field

Taking the measured temperature at the depth of 0.25 m and 2.50 m of cross-section 1 as an example, the temperature parameters in Table 1 and paper [9] are substituted into Eq. (11), Eq. (15), Eq. (17), Eq. (18) and Eq. (26) to calculate the temperature distribution of surrounding rock with and without phase change parameters. The calculation results are shown in Fig. 9.Fig. 9 Surrounding rock temperature of monitoring section 1 and verification. (a) at the depth of 1 m; (b) at the depth of 2.50 m.

Fig. 9

As can be seen in Fig. 9(a), the calculated temperature curve with and without phase change at the depth of 0.25 m is basically consistent with the measured temperature curve, either in terms of distribution shape or change amplitude. The maximum and minimum values of surrounding rock measured temperature at 0.25 m is 11.72 °C and −12.28 °C, and that for calculation without phase change is 11.00 °C and −10.73 °C, with the corresponding difference of 0.72 °C and 1.55 °C, respectively. The maximum and minimum values for calculation with phase change is 10.75 °C and −10.44 °C, with the corresponding difference of 0.97 °C and 1.84 °C, respectively. It can be seen from the data that whether the phase change is considered or not, there is little difference between the calculated results and the measured values in the shallow layer. The total pore water content of the shallow rock mass is small, and the latent heat released under the action of low temperature is less. Consequently, there is no significant difference in the calculation results with and without phase change in the shallow rock mass.

It can be seen from Fig. 9(b) that the calculated temperature curve with and without phase change at the depth of 2.5 m is basically consistent with the measured temperature in distribution shape, but differs in the change amplitude. The maximum and minimum values of surrounding rock measured temperature at 2.50 m is 2.85 °C and −3.03 °C, and that for calculation without phase change is 5.98 °C and −5.47 °C, with the corresponding difference of 3.13 °C and 2.44 °C, respectively. The data shows that there is a large difference between the calculated value without phase change and the measured value. And it is because the energy of the low-temperature air is gradually absorbed by the surrounding rock during its radial propagation of the tunnel, and the pore water in the surrounding rock freezes into ice to release the latent heat, which further consumes energy, resulting in a small calculated value of the surrounding rock temperature without considering phase change. With the phase change in consideration, the calculated maximum and minimum values are 4.39 °C and −3.74 °C respectively, and the corresponding differences with the measured values are 1.54 °C and 0.71 °C respectively.

According to Eq. (23), in the absence of phase change, the apparent heat capacity in the freezing zone is equivalent to the effective heat capacity, which is 1δiρi+δwρw(δiρiCi+δwρwCw). Substituting this result into Eq. (7) yields the thermal diffusion coefficient of the freezing zone without considering phase change as θsks+(1−θs)(δiki+δwkw)θsρsCp,s+(1−θs)(δiρiCi+δwρwCw). The thermal diffusion coefficients considering and not considering phase change were calculated separately, with the results illustrated in Fig. 10. The figure shows that the range of thermal diffusion coefficient in the freezing zone considering phase change is between 5.9 × 10−8 m2/s and 8.2 × 10−8 m2/s, while the range without considering phase change is between 1.1 × 10−6 m2/s∼1.3 × 10−8 m2/s, approximately 17 times that of the phase change-considered diffusivity. According to the thermal diffusion coefficients sensitivity analysis in Section 4.1, a larger thermal diffusion coefficient corresponds to a greater temperature amplitude. In conjunction with Fig. 9(b), the temperature amplitude curve without considering phase change is larger than the actual temperature amplitude curve, whereas the temperature amplitude curve considering phase change closely matches the actual temperature amplitude curve. Therefore, the thermal diffusivity proposed in this paper, which considers phase change, is not only more accurate but also more aligned with practical conditions.Fig. 10 Thermal diffusion coefficient with and without phase change.

Fig. 10

There is still a certain deviation between the calculated result and the measured value, because the air temperature in the selected monitoring section is obtained by fitting the monitored temperature data, and it is different from the actual air temperature. Besides, the joints in the rock mass also affect the heat conduction of the rock mass, but in the analytical calculation, it is assumed that the rock mass is continuous, uniform and isotropic. Overall, the calculation error with phase change is smaller than that without phase change, which shows that the water-ice phase change has a significant impact on the distribution of surrounding rock temperature field. In addition, the error of the calculation considering phase change is within the allowable range, which indicates the correctness and necessity of establishing the calculation model considering phase change and latent heat when calculating the temperature field of tunnel surrounding rock in cold regions.

5.3 Comparative analysis of maximum frozen depth

The maximum frozen depth of surrounding rock of tunnels in cold regions is an important parameter in the design of drainage system and thermal insulation of tunnels in cold regions. The research shows that the frozen depth of surrounding rock is the maximum when the tunnel temperature reaches positive temperature from negative temperature (t=t2). In previous section, the analytical calculation formula of surrounding rock temperature has been obtained, and the maximum frozen depth can be obtained by T(x,t2) = -0.5 °C. Fig. 11(a) shows the calculated values and error values of the maximum frozen depth of the 15 monitoring sections of the Xing'anling Tunnel. It can be seen that the maximum frozen depth error calculated with phase change is far less than that calculated without phase change.Fig. 11 Frozen depth values and calculation errors. (a) phase change and non-phase change; (b) phase change and the modified Stefan formula.

Fig. 11

Notably, the maximum frozen depth error calculated at cross- sections 3, 7 and 15 considering phase change is relatively large, which is 1.4 m, 1.83 m and 1.05 m, respectively. The reasons for the errors are as follows. First, the rock mass status and composition of different cross-sections are different. The total length of Xing'anling Tunnel is 3960 m, and the 15 monitoring sections cannot be of the same rock mass. Second, measurement errors may be present. Fig. 12 shows the annual average temperature and annual temperature amplitude of different monitoring sections. Take cross-section 3 as an example for further illustration. At monitoring section 3, the annual average temperature of the surface environment is −0.45 °C, and the annual amplitude is 11.67 °C. Similarly, cross-section 1 has an annual average temperature of 0.07 °C and an annual amplitude of 11.73 °C. Compared with cross-section 1, cross-section 3 has an even lower annual minimum temperature and a longer period of negative temperature. Theoretically, this should result in a greater frozen depth for cross-section 3 than for cross-section 1. However, the measured frozen depth of cross-section 3 is 4.65 m, which is less than that of cross-section 1, which measures 5.13 m. The same discrepancy is observed in cross- sections 7 and 15.Fig. 12 Annual average temperature and annual temperature amplitude of different monitoring sections.

Fig. 12

Exclude the data with large errors of cross- sections 3, 7 and 15, adopt the same surrounding rock parameters and environmental parameters and calculate the maximum frozen depth using the modified Stefan formula. We compared the frozen depth values and calculation errors between method considering phase change and method of modified Stefan formula, and the results are shown in Fig. 11(b). The error between calculated values by our method and the measured values on cross- sections 1, 8, 9, 10, 11, 12 and 14 is less than the error between the calculated values by the modified Stefan formula and the measured values. The average error value of frozen depth obtained by our analytical solution method (0.23 m) is less than that by the modified Stefan formula (0.34 m).

To sum up, after excluding some data with large errors, the frozen depth obtained by the analytical solution of the temperature field of surrounding rock with phase change proposed in this paper is in good agreement with the field measured results, and the calculation accuracy is better than that of the modified Stefan formula. Therefore, the analytical solution method of temperature field of surrounding rock with phase change can meet the requirements of tunnel drainage system and thermal insulation design in cold regions.

6 Discussion

Based on the above research and results, the heat transfer calculation formula of tunnel surrounding rock with phase change in cold regions established in this paper can truly reflect the phenomenon of water-ice phase change, and accurately present the radial heat transfer distribution of tunnel temperature in cold regions. Compared with the analytical solution formula without phase change, the key of the analytical calculation proposed in this paper is to calculate the average thermal diffusion coefficient αr with phase change. Eq. (26) is the calculation formula of the average thermal diffusion coefficient. Its calculation process is complex and requires numerous parameters, which is not convenient for engineering application. In order to facilitate engineering application, Eq. (26) is simplified and the resulting error is further analyzed. Eq. (26) is composed of two parts and calculation results showed that the value of the first half is far less than that of the second half. Therefore, Eq. (26) can be simplified as follows(27) αr=βαe

Where β is simplified coefficient and β=1−t2−t1365 (see t1 and t2 in Fig. 4), αe is the thermal diffusion coefficient of rock mass at normal temperature (m2/s). If there is no measured value, please query the thermodynamic parameters of the soil and then substitute them into αe=θsks+θfkwθsρsCp,s+θfρwCp,w.

The thermal diffusion coefficient of surrounding rock of Xing'anling Tunnel at normal temperature is 1.1×10−6 and the values of thermal diffusion coefficient calculated by Eq. (26) and Eq. (27) are shown in Fig. 13(a). Clearly, there are minor differences of the thermal diffusion coefficient between the values calculated by the simplified formula and the original formula. The calculation error at cross-section 3 is the largest, only 2.16 %. Further, the maximum frozen depth of 15 cross-sections is calculated by using the thermal diffusion coefficients of Eq. (26) and Eq. (27) respectively, and the results are shown in Fig. 13(b). The frozen depth calculated by the simplified thermal diffusion coefficient is almost the same as that calculated by the average thermal diffusion coefficient with phase change. Cross-section 3 has the maximum error of 0.06 m and other sections have an error of less than 0.02 m. Therefore, the simplified thermal diffusion coefficient proposed in this paper can replace the average thermal diffusion coefficient with phase change in the design of tunnels in cold regions.Fig. 13 The theoretical calculation method and the simplified calculation method. (a) the average thermal diffusion coefficient; (b) the maximum frozen depth.

Fig. 13

7 Conclusions

The refined equivalent method of water-ice phase change and latent heat of tunnel surrounding rock in cold regions and the heat transfer model of surrounding rock considering phase change are established. The theoretical calculation formula of the temperature field of tunnel surrounding rock with phase change in cold regions is presented. Factors affecting the temperature of surrounding rock with phase change are deeply studied, and the engineering verification and related discussions are carried out. Main conclusions are drawn as follows.(1) Based on the apparent heat capacity method, a refined method is proposed to convert the latent heat of water-ice phase change into heat capacity in cold-region tunnels, the transient heat transfer equation considering phase change is derived, and the theoretical model of surrounding rock temperature field considering phase change in cold-region tunnels is established. Based on the idea of equal accumulated temperature, the average thermal diffusion coefficient is defined, which overcomes the problem that the classical heat transfer theory cannot directly solve the continuous temperature field of “unfrozen-freezing-frozen” state. Combined with the separation variable method and Fourier integral transformation method, the analytical formula for calculating the transient temperature field of surrounding rock of tunnels in cold regions considering phase change is presented. The model and formula can fully consider the influence of the water-ice phase change and the impact of the initial temperature of surrounding rock, thermal diffusion coefficient, annual average temperature and annual temperature amplitude.

(2) The frozen range of surrounding rock with phase change in cold-region tunnel increases with the increase of thermal diffusion coefficient and annual temperature amplitude, and decreases with the increase of annual average temperature. The peak temperature of surrounding rock increases with the increase of thermal diffusion coefficient, annual temperature amplitude and annual average temperature; The valley temperature of surrounding rock decreases with the increase of thermal diffusion coefficient and annual temperature amplitude, and increases with the increase of annual average temperature. The larger the thermal diffusion coefficient and annual temperature amplitude and the smaller the annual average temperature, the more significant the phase lag effect of the temperature fluctuation of surrounding rock. The initial value of surrounding rock temperature has a weak influence on the distribution of surrounding rock temperature field.

(3) The field engineering application and analysis show that the temperature distribution characteristics obtained by the theoretical calculation method proposed in this paper are in good agreement with the field measurement. The maximum error between the calculation method considering phase change and the measured value is 1.54 °C, which is less than the maximum error between the calculation method not considering phase change and the measured value (3.13 °C). The average error between the maximum frozen depth of surrounding rock calculated by the method proposed in this paper and that obtained by field measurement is 0.23 m, which is smaller than the error between the calculation result of the modified Stefan formula and the measured value (0.34 m).

(4) From the perspective of convenient engineering application, the simplified coefficient of thermal diffusivity is deduced and defined, and the simplified form of analytical solution of the temperature field with phase change in the surrounding rock of the tunnel in cold regions is presented. The simplified coefficient can be calculated by using the air temperature inside the cold-region tunnel, and the thermal diffusion coefficient of rock mass at normal temperature is modified. In this way, the effect of water-ice phase change can be considered in the calculation of tunnel surrounding rock temperature field in cold regions, which effectively improves the calculation speed. The maximum frozen depth obtained by the simplified form is basically consistent with the results obtained by the theoretical calculation method of the temperature field of the surrounding rock with phase change in the cold-region tunnel proposed in this paper, and the maximum error is only 0.02 m.

Funding

This research was funded by the 10.13039/501100001809 National Natural Science Foundation of China (52178388 , 52208385 ), 10.13039/501100002858 China Postdoctoral Science Foundation (2018M631114 ) and bureau-level key project of Henan Province (2020-4-17).

Data availability statement

The data used to support the findings of this study are included within the article.

CRediT authorship contribution statement

Wentao Wu: Writing – original draft, Visualization, Methodology, Formal analysis, Data curation, Conceptualization. Jiaqi Guo: Validation, Funding acquisition, Data curation. Xiaochuan Wang: Writing – review & editing, Data curation. Huanmeng Hu: Writing – review & editing, Validation. Pengyu Zhao: Writing – review & editing, Visualization, Validation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
==== Refs
References

1 Ma Q.G. Luo X.X. Lai Y.M. Niu F.J. Gao J.Q. Numerical investigation on thermal insulation layer of a tunnel in seasonally frozen regions Appl. Therm. Eng. 138 2018 280 291
2 Kang F.C. Li Y.C. Tang C.A. Numerical study on airflow temperature field in a high-temperature tunnel with insulation layer Appl. Therm. Eng. 179 2020 115654
3 Gao J.Q. Lai Y.M. Zhang M.Y. Chang D. The thermal effect of heating two-phase closed thermosyphons on the high-speed railway embankment in seasonally frozen regions Appl. Therm. Eng. 141 2018 948 957
4 Lai J.X. Wang X.L. Qiu J.L. Zhang G.Z. Chen J.X. Xie Y.L. Luo Y.B. A state-of-the-art review of sustainable energy based freeze proof technology for cold-region tunnels in China Renew. Sustain. Energy Rev. 82 2018 3554 3569
5 Chen J.X. Zhao P.Y. Luo Y.B. Deng X.H. Liu Q. Damage of shotcrete under freeze-thaw loading J. Civ. Eng. Manag. 23 2017 583 593
6 Lai Y.M. Wu Z.W. Zhu Y.L. Zhu L.N. Nonlinear analysis for the coupled problem of temperature, seepage and stress fields in cold-region tunnels Tunn. Undergr. Space Technol. 13 1998 435 440
7 Tan X.J. Chen W.Z. Tian H.M. Cao J.J. Water flow and heat transport including ice/water phase change in porous media: numerical simulation and application Cold Reg. Sci. Technol. 68 2011 74 84
8 Zhao P.Y. Chen J.X. Luo Y.B. Li Y. Chen L.J. Chen L.J. Wang C.W. Field measurement of air temperature in a cold region tunnel in northeast China Cold Reg. Sci. Technol. 171 2020 102957
9 Zhao X. Zhang H.W. Lai H.P. Yang X.H. Wang X.H. Zhao X.L. Temperature field characteristics and influencing factors on frost depth of a highway tunnel in a cold region Cold Reg. Sci. Technol. 179 2020 103141
10 Zhou X.H. Zeng Y.H. Fan L. Temperature field analysis of a cold-region railway tunnel considering mechanical and train-induced ventilation effects Appl. Therm. Eng. 100 2016 114 124
11 Zeng Y. Liu K. Zhou X. Fan L. Tunnel temperature fields analysis under the couple effect of convection-conduction in cold regions Appl. Therm. Eng. 120 2017 378 392
12 Zhou Y. Zhang X. Deng J. A mathematical optimization model of insulation layer's parameters in seasonally frozen tunnel engineering Cold Reg. Sci. Technol. 101 2014 73 80
13 Yang P. Ke J.M. Wang J.G. Chou Y.K. Zhu F.B. Numerical simulation of frost heave with coupled water freezing, temperature and stress fields in tunnel excavation Comput. Geotech. 33 2006 330 340
14 Yan Q.X. Li B.J. Zhang Y.Y. Yan J. Zhang C. Numerical investigation of heat-insulating layers in a cold region tunnel, taking into account airflow and heat transfer Appl. Sci. 7 2017 679
15 Tan X.J. Chen W.Z. Wu G.J. Yang J.P. Numerical simulations of heat transfer with ice–water phase change occurring in porous media and application to a cold-region tunnel Tunn. Undergr. Space Technol. 38 2013 170 179
16 Zhang Y.W. Song Z.P. Lai J.X. Spatial-temporal distribution of 3D temperature field of highway tunnel in seasonal frozen soil region J. Highw. Transp. Res. Dev. 37 2020 96 103
17 Zhang X.F. Lai Y.M. Yu W.B. Zhang S.J. Zhang J.Z. Forecast analysis of the refreezing of Kunlun mountain permafrost tunnel on Qing–Tibet railway in China Cold Reg. Sci. Technol. 39 2004 19 31
18 Zhan Y.X. Lu Z. Yao H.L. Numerical analysis of thermo-hydro-mechanical coupling of diversion tunnels in a seasonally frozen region J. Cold Reg. Eng. 34 2020 04020018
19 Ma Q.G. Luo X.X. Lai Y.M. Niu F.J. Gao J.Q. Numerical investigation on thermal insulation layer of a tunnel in seasonally frozen regions Appl. Therm. Eng. 138 2018 280 291
20 Li S.Y. Niu F.J. Lai Y.M. Pei W.S. Yu W.B. Optimal design of thermal insulation layer of a tunnel in permafrost regions based on coupled heat-water simulation Appl. Therm. Eng. 110 2017 1264 1273 2017
21 Zhu S. Cheng J.W. Song W.T. Borowski M. Zhang Y.J. Using seasonal temperature difference in underground surrounding rocks to cooling ventilation airflow: a conceptual model and simulation study Energy Sci. Eng. 8 2020 3457 3475
22 Lai Y.M. Liu S.Y. Wu Z.W. Yu W.B. Approximate analytical solution for temperature fields in cold regions circular tunnels Cold Reg. Sci. Technol. 34 2002 43 49
23 Zhang Y. He S.H. Li J.B. Analytic solutions for the temperature fields of a circular tunnel with insulation layer in cold region J. Glaciol. Geocryol. 31 2009 113 118
24 Xia C.C. Zhang G.Z. Xiao S. Analytical solution to temperature fields of tunnel in cold region considering lining and insulation layer Chin. J. Rock Mech. Eng. 29 2010 1767 1773
25 Lu X. Tervola P. Viljanen M. A new analytical method to solve the heat equation for a multi-dimensional composite slab J. Phys. Math. Gen. 38 2005 2873
26 Liu W.W. Feng Q. Wang C.X. Lu C.K. Xu Z.Z. Li W.T. Analytical solution for three-dimensional radial heat transfer in a cold-region tunnel Cold Reg. Sci. Technol. 164 2019 102787
27 Gao Y. Ding Y.F. Feng Y. Xia J.J. Tian X.Z. Geng J.Y. Research on the length of the anti-freezing layer of cold-region tunnels under the influence of train-induced wind Tunn. Undergr. Space Technol. 146 2024 105665
28 Konrad J.K. Effects of applied pressure on freezing soils Can. Geotech. J. 19 1982 494 505
29 Nakano Y. Quasi-steady problems in freezing soils: I. Analysis on the steady growth of an ice layer Cold Reg. Sci. Technol. 17 1990 207 226
30 Hashemi H.T. Sliepcevich C.M. A numerical method for solving two-dimensional problems of heat conduction with phase change Chem. Eng. Prog. Symp. Ser. 63 1967 34 41
31 Harlan R.L. Analysis of coupled heat-fluid transport in partially frozen soil Water Resour. Res. 9 1973 1314 1323
32 Zhu S. Wu S.Y. Cheng J.W. Li S.Y. Li M.M. An underground air-route temperature prediction model for ultra-deep coal mines Minerals 5 2015 527 545
33 Lu J.G. Wan X.S. Yan Z.R. Qiu E.X. Pirhadi N.M. Liu J.N. Modeling thermal conductivity of soils during a freezing process Heat Mass Tran. 58 2022 283 293
34 Tan X.J. Chen W.Z. Jia S.P. Lv S.P. A coupled hydro-thermal model for low temperature rock including phase change Chin. J. Rock Mech. Eng. 27 2008 1455 1461
