Numerical analysis for permafrost temperature field in the short term of permafrost subgrade filling

Purpose – It isof greatsignificancetostudythe influenceofsubgradefillingonpermafrosttemperaturefield in permafrost area for the smooth construction and safe operation of railway. Design/methodology/approach – Thepaperbuildsupthemodelforthehydrothermalcouplingcalculation ofpermafrostusingfiniteelementsoftwareCOMSOLtostudyhowpermafrosttemperaturefieldchangesintheshorttermaftersubgradefilling,onwhichbasisitproposesthemethodofcalculationfortheconcavedistortionoffreezingfrontinthesubgrade-coveredarea. Findings – The results show that the freezing front below the subgrade center sinks due to the thermal effect of subgrade filling, which will trigger hydrothermal erosion in case of sufficient moisture inflows, leading to the thawing settlement or the cracking of the subgrade, etc. The heat output of soil will be hindered the most in case of July filling, in which case the sinking and the distortion of the freezing front is found to be the most severe, which the recovery of the permafrost temperature field, the slowest, constituting the most unfavorable working condition. The concave distortion of the freezing front in the subgrade area increases with the increase in temperature difference between the filler and ground surface, the subgrade height, the subgrade width and the volumetric thermal capacityoffiller, while decreases with the increase of the thermal conductivity of filler. Therefore, the filler chose for engineering project shall be of small volumetric thermal capacity, low initial temperature and high thermal conductivity whenever possible. Originality/value – The concave distortion of the freezing front under different working conditions at different times after filling can be calculated using the method proposed.


Introduction
With an average elevation of more than 4,500 m, 70% of Qinghai-Tibet Plateau is covered by permafrost (Cheng, 2003a, b;Yang, Meng, Han, Li, & Cai, 2018). As a special type of soil, permafrost is highly sensitive to temperature changes and displays a close relation between its physical and mechanical properties and soil temperature. Any engineering activities, including subgrade construction, on permafrost will inevitably disrupt the original balance between the permafrost and the environment, as well as the thermal state of permafrost field, undermining the Permafrost temperature field subgrade stability and that of other structures (Cheng, 2003a, b;Mi, Zhao, Yang, Qu, & Wu, 2017;Tian, Zhang, Mu, & Liu, 2014;Zhang, 2000). An accumulative length of 632 km of Qinghai-Tibet Railway is built on permafrost, while the building of Sichuan-Tibet Railway also faces the same issue of subgrade stability. Therefore, it is of great significance to study the influence of subgrade filling on the temperature field of permafrost to ensure the building progress and the operation safety of both railways. Both field monitoring and testing are important approaches to study the influence of subgrade filling on permafrost temperature field. Wang, Mi, Wei, and Wu (2010) analyzed the monitoring data of eleven subgrade sections of Qinghai-Tibet Railway in a span of 7 years and found that the artificial permafrost table of the subgrade center rose with time, which was conducive to the subgrade stability. Niu, Ma, and Wu (2011) investigated the monitoring results of 55 sections along Qinghai-Tibet Railway and found that the regular subgrade in some warm permafrost areas displayed descending upper limit and poor thermal stability. That aside, the thermal disturbance and hydrothermal erosion at the shallow ground caused by the construction were prone to trigger thermal thawing disasters, leading to subgrade settlement and other defects. Mu, Ma, Niu, Li, and Wang (2014) analyzed the monitoring data of the ground temperature at a regular subgrade section of Qinghai-Tibet Railway from 2005 to 2010 and found that the thermal stability of the subgrade permafrost was closely related to its average temperature, and the thermal stability of frozen soil under the subgrade was poor in high temperature permafrost regions. Zhang, Wen, Xue, Chen, and Li (2016) found that the permafrost table for the sunny slope descended at the initial operation stage. The water seepage not only accelerated the temperature rise and deformation process of the subgrade but also increased the non-uniformity of the temperature and moisture fields of the subgrade. Zhao, Cheng, Han, and Cai (2019) also conducted relevant research. However, existing studies generally focus on the influence of subgrade filling on the long-term thermal stability of permafrost, while only few examine how the temperature field changes in the short term. In addition, it is quite difficult to monitor the short-term concave distortion of the freezing front. It is said that the thermal effect of subgrade filling results in short-term concave distortion in the freezing front of the subgrade and even continuous defect in case of sufficient moisture inflows, leading to thawing settlement or cracking of the subgrade, etc. Therefore, in-depth study is necessary.
In recent years, a large number of scholars have used numerical simulation to study the permafrost-related issues, the key to which is a reasonable and effective hydrothermal coupling model. Harlan (1973) proposed the concept of hydrothermal coupling based on the moisture migration theory and established a hydrothermal coupling model with permeability coefficient and ice content incorporated as parameters for the first time. Xu, Wang, and Zhang (2001) provided an empirical formula for the content of unfrozen water based on a large number of experimental studies. Liu and Tian (2002), and Tian, Liu, Qian, and He (2002) simulated and analyzed the influence of latent heat amid phase transition using the enthalpy of soil state as a variable. Song, Tong, Luo, and Li (2019) proposed a water-gas migration control equation with overall consideration given to the influence of hydraulic and thermal gradients based on Philip's non-isothermal water-heat-gas coupled migration model. In addition, Bai, Li, Tian, and Fang (2015), and Hu, Wang, Chang, Liu, and Lu (2020) also modified and improved the hydrothermal coupling model of permafrost and applied it to numerical simulation.
Built on the control equation for hydrothermal coupling of permafrost, the paper introduces the coefficient-type partial differential equation module concluded by the finite element software COMSOL and establishes the model for hydrothermal coupling analysis of permafrost, so as to study how the permafrost temperature field changes in the short term after subgrade filling.

Control equation for hydrothermal coupling of permafrost
(1) Control equation for moisture field Both the frozen and thawed soil contain a certain amount of unfrozen water. According to Philip's non-isothermal water-heat-gas coupled migration model (Philip & De Vries, 1957) moisture migration is mainly caused by hydraulic and thermal gradients. The control equation for water-gas migration proposed in the reference (Song et al., 2019) fails to consider the migration of vaporous water, and the control equation for the migration of liquid water in unsaturated soil can be expressed as where θ u and θ ithe volume fraction of unfrozen water and ice, respectively; ttime, in s; ρ ithe density of ice, taken as 900 kg$m À3 ; ρ wthe density of water, taken as 1 000 kg$m À3 ; K hthe hydraulic conductivity caused by the hydraulic gradient under isothermal conditions, in m$s À1 ; K Tthe hydraulic conductivity caused by the thermal gradient under non-isothermal conditions, in m 2 $K À1 $s À1 ; hthe hydraulic head corresponding to the matric suction, in m; h ghydraulic head caused by gravity, i.e. elevation head, in m; Tthe absolute temperature, in K; ∇the differential operator.
The effective saturation S is defined as θu − θr θs − θr , and its relations with the hydraulic head corresponding to the matric suction h can be expressed as It needs to be noted that α, n and m are parameters of the model to be determined by tests. With α measured in m −1 , θ s and θ r are saturated volumetric moisture content and residual volumetric moisture content, respectively. Given that the ice resistance in the ice-water mixing zone amid the migration, K h can be expressed as where k sthe permeability coefficient of soil.

Permafrost temperature field
The formula for hydraulic conductivity caused by the thermal gradient under non-isothermal conditions K T (Saito, Simunek, Scanlon, & Mohanty, 2006) is (2) Control equation for temperature field The paper waives the influence of heat consumption in water evaporation during soil freezing and thawing and takes the latent heat of fusion as a heat source. According to the Fourier's law, the control equation for the transient temperature field of unsaturated soil can be expressed as follows: where Cthe volumetric thermal capacity of soil, in J$(m 3 $K) À1 ; λthe thermal conductivity of soil, in W$m À1 $K À1 ; Lthe latent heat of fusion, taken as 334.5 kJ$kg À1 .
Both the volumetric thermal capacity and thermal conductivity of permafrost are determined by soil skeleton, unfrozen water, ice and porous gas. For the sake of the calculation efficiency and convergence, the paper provides only the thermodynamic parameters of the freezing and thawing soil, in other words where C f and C ufthe volumetric thermal capacity of frozen and thawed soil, in J$m À3 $K À1 ; λ f and λ ufthe thermal conductivity of frozen and thawed soil, in W$m À1 $K À1 ; T fthe freezing temperature of soil, in K (3) Unfrozen water content The control equations for the moisture and temperature fields contain three unknown quantities, namely, the volume fractions of unfrozen water, ice and the temperature. To solve a hydrothermal coupling problem, it is still necessary to introduce an equation targeting the three quantities. According to the reference (Hu et al., 2020), the relationship between ice content and the unfrozen water content can be expressed as follows: where Bthe empirical parameter.

Verification of control equation for hydrothermal coupling of permafrost
The paper applies the finite element software COMSOL to simulate the soil column thawing test under closed conditions in the reference (Xu et al., 2001) to verify the effectiveness of the control equation for the hydrothermal calculation of permafrost. The thawing test in the reference (Xu et al., 2001) applies the soil column made of Lanzhou silty soil records 10 cm in diameter, 10 cm in height, 22.3% in initial volumetric moisture content and À1 8C in initial temperature. During the test, the bottom temperature of the soil column is kept at À1 8C, while the top, at 1 8C.
The paper uses COMSOL to build up a one-dimensional soil column model with a height of 10 cm, while the soil is simulated using a quadratic Lagrangian element, in which case the initial values and boundary conditions of the moisture and temperature fields are identical to those in the reference (Xu et al., 2001). Table 1 shows the parameters of the silty soil simulated in accordance with the references (Bai et al., 2015;Song et al., 2019). The moisture content distribution and temperature distribution of the specimens obtained from the test and numerical simulation in the reference (Xu et al., 2001) are shown in Figures 1 and 2.
It can be seen from Figures 1 and 2 that the simulation results are basically consistent with the measured values of the test. The moisture content remains generally constant in the upper half of the soil column, or to be more specific smaller than the initial value. This can be explained by the downward migration of liquid water under the combined action of hydraulic   Permafrost temperature field and thermal gradients after the thawing of the upper frozen soil. On the contrary, the moisture content in the lower half of the soil column increases first and then decreases with the decrease in temperature. The moisture content peaks as the large hydraulic gradient and the fast moisture migration near the freezing temperature, as well as the resistance of ice make it difficult for liquid water to continue downward. Therefore, water accumulates in large amounts near the freezing temperature, constituting moisture content peak, while the soil column temperature displays a generally linear distribution, in line with the actual situation. It can be seen that the control equation for hydrothermal calculation of permafrost proposed can accurately simulate the changes of moisture migration and temperature field in frozen soil, therefore the proposed approach can be used to simulate the subgrade filling issue in permafrost regions.

Modeling
In reference to the subgrade section of the Beiluhe test section on Qinghai-Tibet Railway (Hu et al., 2020;Philip & De Vries, 1957), the paper builds up the subgrade model using the finite element software COMSOL where the foundation width and depth are deemed at 50 and 20 m, respectively, subgrade width and height, 8 and 4 m, and slope ratio fixed at 1:1.5. To improve the calculation efficiency, half of the symmetrical model is taken for calculation, as shown in Figure 3. The subgrade and foundation are simulated using triangular elements. The key parameters are shown in Table 2.
The temperature boundary of the model is determined by relevant monitoring data of the Beiluhe test section of Qinghai-Tibet Railway in the references (Niu et al., 2011;Bian, 2012;Sun, 2018). The left and right boundaries of the model temperature field are adiabatic boundaries, and the Dirichlet boundary is introduced for the bottom (taken as À0.8 8C). The top temperature T top changes with time, i.e.
Foundation soil 3 3 10 À7 1 800 22.5 1 800 2 000 1.2 1.0 Subgrade filler 3 3 10 À5 2 000 12.0 1 800 2 000 1.2 1.0 Source(s): Authors own work  Given that the projects in plateau and cold areas are generally carried out in the warm seasons (between April and October), it can be concluded that the corresponding time when t 5 0 s in Equation (9)  The initial values of the moisture and temperature fields of the foundation soil are arrived in calculation. In the absence of subgrade, the initial temperature and the volume fraction of the unfrozen water of the foundation soil are À0.8 8C and 0.12 respectively. At this point, the boundary conditions described above are used to calculate the initial values of the temperature and moisture fields of the foundation soil after 50 a.
After transportation, compaction and other processes, the filler temperature is often higher than the average surface temperature, and for safety consideration, the figure can be as high as 10 8C higher than the average surface after subgrade filling. The volume fraction of unfrozen water of the filler is 0.12. The paper focuses on the thermal effect of subgrade filling in the short term, and the time period is taken as 90 d.
5. Permafrost temperature field in the short term after subgrade filling Subgrade filling exerts significant impact on the temperature field of permafrost in the short term, in which case the influencing factors mainly include the initial conditions, subgrade dimensions and thermodynamic properties of fillers. Initial conditions include the time of subgrade filling and the temperature difference between the filler and ground surface. For the foundation, the surface temperature continues to rise between April and July, during which the foundation soil absorbs heat from the surface, while as the surface temperature gradually decreases from July to October, the foundation soil releases heat to the surface. Therefore, the three key time nodes -April, July and Octoberare selected for the calculation to study the impact of filling time. Subgrade dimensions include subgrade height and width. Thermodynamic properties of the filler include volumetric thermal capacity and thermal conductivity. Figures 4 and 5 show the temperature fields 30 and 60 d after subgrade filling in April, July and October, respectively. The solid red line represents the freezing fronts, i.e. the 0 8C isotherms.

Influence of filling time
It can be seen from Figures 4 and 5 that the temperature peak appears in the center of the subgrade 30 d after filling, indicating that the heat of the filler continues to transfer to the atmosphere and the foundation; 60 d after filling, the temperature of the filled subgrade gradually decreased from the ground surface to the foundation in April, consistent with the temperature distribution of natural foundation, indicating that the heat of the subgrade filler has been basically dissipated. In July and October, the warm areas of the filled subgrade are found at the subgrade center, the temperature of which is significantly higher than that of the natural foundation, meaning that the heat of filler continues to transmit to the atmosphere and the foundation.
The ice resistance makes the permeability coefficient of frozen soil significantly lower than that of the thawed soil. The water obtained from the thawing of frozen soil and that from the surface into the foundation tends to accumulate near the freezing front, making it quite Permafrost temperature field difficult to infiltrate downward. Under the influence of gravity, the water flows and accumulates at the depression along the freezing front. The water infiltrated from the surface records high in temperature. As the water meets at the depression, the frozen soil starts to thaw under thermal erosion, leading to the aggravated depression and the sinking of the freezing front, which in turn translates into more water accumulation. In practice, the thermal effect of subgrade filling leads to the sinking of the freezing front in the subgrade area. Even in case of minor distortion, that is, the freezing front is only a few centimeters away from the natural freezing front, it may cause continuous concave distortion under sufficient moisture inflows, resulting in thawing settlement, cracking and other defects of the subgrade. Therefore, compared with the absolute elevation of the freezing front, the RS concave distortion after subgrade filling emerges as a more important concern, or in other words the elevation difference between the freezing front of the subgrade area and that of the natural foundation. Figure 6 shows the development of the elevation difference with the change of time. Positive values indicate that the freezing front is higher than the natural foundation after subgrade filling, while negative values, lower, which indicates that the freezing front is subject to concave distortion.
It can be seen from Figure 6(a) that when the filling is carried out in April, the concave distortion starts at the freezing front at the subgrade center, with the maximum depth of 20 cm. The distortion gradually decreased with time, resulting in convex distortion in the end. The slope toe continued to bulge upward after filling. The convex distortion of the freezing front at the slope toe helps keep external water from flowing down into the subgrade, which is conducive from an engineering viewpoint. However, such water may flow down the slope surface and accumulate outside the convex slope toe, causing hydrothermal erosion. It can be seen from Figure 6(b) that when the filling is carried out in July, the concave distortion of the permafrost in the subgrade area continued to develop with the elevation difference at the subgrade center increasing from 0 to 25 cm. Only 90 d after filling is slight convex distortion observed at the subgrade slope toe. It needs to be noted that summer time generally features abundant rainfall. In case of sufficient surface water, high temperatures and excessive heat, hydrothermal erosion is likely to take place at permafrost area subject to concave distortion, one of the greatest safety concerns for projects.
It can be seen from Figure 6(c) that when the filling is carried out in October, the concave distortion of the permafrost in the subgrade area continued to aggregate with the elevation difference at the subgrade center increasing from 0 to 38 cm. No convex distortion is observed at the slope toe. Given the insufficient surface water, low temperature and the fact that the surface temperature drops quickly below 0 8C after filling, hydrothermal erosion is less likely Figure 6. Temporal development of elevation difference of freezing front in subgrade area RS to occur, though the concave distortion of the permafrost is most severe in the same period after filling. Figure 7 shows the temporal development of elevation difference after filling in April, July and October. It can be seen that the development pattern observed in April is quite different from those in July and October. In April, the freezing front at the subgrade center first drops and then increases after filling, while in July and October, it continues to drop, that is, the concave distortion in the subgrade area continues to advance. In case of July, the absolute value of the elevation difference shows an approximate logarithmic relation with time; while for October, the elevation difference increases sharply with the change of time until a closed thawed pit was formed in the subgrade area. The subsequent studies will focus on the filling in July. Figures 8 and 9 show the concave distortion curves of the freezing front in the subgrade area and the elevation difference between the subgrade center and the natural ground 90 d after filling when the temperature difference between the filler and ground surface varies.

Influence of temperature difference between filler and ground surface
It can be seen from Figure 8 that the concave distortion of the freezing front intensifies with the rising temperature difference between the filler and ground surface. When the temperature  difference is 4 8C, obvious convex distortion of the freezing front is observed at the slope toe, which disrupts the accumulation of surrounding water at the bottom of the subgrade center. It can be seen from Figure 9 that the elevation difference between the subgrade center and the natural ground increases with the rising of temperature differences. Figures 10 and 11 show the concave distortion curves of the freezing front in the subgrade area and the elevation difference between the subgrade center and the natural ground 90 d after filling in case of different subgrade heights.

Influence of subgrade dimensions
It can be seen from Figure 10 that the concave distortion of the freezing front in the subgrade area intensifies as the subgrade height increases.
It can be seen from Figure 11 that the elevation difference between the subgrade center and the natural ground increases with the rising of subgrade height. Figures 12 and 13 show the concave distortion of the freezing front in the subgrade area and the elevation difference between the subgrade center and the natural ground 90 d after filling in case of different subgrade width.

RS
It can be seen from Figure 12 that the concave distortion of the freezing front in the subgrade area develops with the increase of subgrade width.
It can be seen from Figure 13 that the elevation difference between the subgrade center and the natural ground intensifies with the increase of subgrade width. Figures 14 and 15 show the concave distortion curves of the freezing front in the subgrade area and the elevation difference between the subgrade center and the natural ground 90 d after filling in case of different thermal conductivity of filler (in the frozen state).

Influence of thermodynamic properties of filler
It can be seen from Figure 14 that the concave distortion of the freezing front in the subgrade area decreases with the increase of the thermal conductivity of filler. When the conductivity reaches 2.0 W$m À1 $8C À1 , obvious convex distortion can be observed at the slope toe, which disrupts the accumulation of surrounding water at the bottom of the subgrade center.  It can be seen from Figure 15 that the elevation difference between the subgrade center and the natural ground decreases as the thermal conductivity of filler goes up. Figures 16 and 17 show the concave distortion curves of the freezing front in the subgrade area and the elevation difference between the subgrade center and the natural ground 90 d after filling in case of different volumetric thermal capacity of filler (in the frozen state).
It can be seen from Figure 16 that the concave distortion of the freezing front in the subgrade area advances with the increase of the volumetric thermal capacity of filler.
It can be seen from Figure 17 that the elevation difference between the subgrade center and the natural ground goes up with the increase of the volumetric thermal capacity of filler.
6. Calculation method of freezing front shortly after subgrade filling Based on the above research results, a simple and practical approximate calculation method of the concave distortion curve of the freezing front is presented. This method can quickly estimate the deformation of freezing front according to some parameters and provide basic data for calculating the deformation of subgrade. This method is described in detail below.
The concave distortion curve of the freezing front can be approximated as an inverted trapezoid, as shown in Figure 18.

RS
The curve can be divided into three segments, namely AB, BC and CD. Segment AB refers to the central part of the subgrade. The elevation difference in this area is the most significant and generally maintains constant, immune to the change of the horizontal distance. Segment BC refers to the status of the part from the subgrade shoulder to the subgrade slope toe, the elevation difference of which decreases rapidly in an approximately linear manner. Segment CD refers to the subgrade slope toe and adjacent areas, and the elevation difference in this area is largely zero. Therefore, the elevation difference Δh of the freezing front can be expressed as where Δh maxthe elevation difference of the freezing front at the subgrade center, i.e. OA in Figure 18; W 1the width of the segment in which the elevation difference is constant at the subgrade center, i.e. AB in Figure 18; W 2the width of the segment in which the elevation difference changes, i.e. OC in Figure 18; Wthe distance between the subgrade center and the slope toe, i.e. OD in Figure 18.
The specific calculation of the above parameters is as follows.
(1) Δh max is the approximate logarithmic with time, which increases with the rise of subgrade height, subgrade width, volumetric thermal capacity of filler as well as the temperature difference between filler and ground surface, while decreases with the increase in the thermal conductivity of filler. According to the calculation results in Figures 7-11, 13, and 15, the fitting results Δh max are as follows: Δh max ¼ ð11:3 ln t f À 26:22Þ½1 þ 0:128ðΔT À 10Þ½1 þ 0:992ðh À 4Þ½1 þ 0:0313 ðW subgrade À 8Þ½1 þ 0:342ð0:001C f À 1:8Þ½1 À 0:672ðλ f À 1:2Þ where t fthe time duration after filling, in d; ΔTthe temperature difference between the filler and ground surface during filling, in 8C; hsubgrade height, in m; W subgradesubgrade width, in m; C fthe volumetric thermal capacity of frozen filler, in J$m À3 $K À1 ; λ fthe thermal conductivity of frozen filler, in W$m À1 $K À1 .
(2) W is determined by subgrade height and width, and calculated by the following equation: Permafrost temperature field (3) W 1 is determined by subgrade height and width, and calculated by the following equation: W 1 ¼ 2 þ 0:75ðh À 4Þ þ 0:5ðW subgrade À 8Þ (13) (4) W 2 increases with the increase of the subgrade height and width, while decreases with the increase of the thermal conductivity of filler. With the calculation results of Figure 14 fitted, the calculation equation is obtained as follows.
W 2 ¼ 9 þ 0:75ðh À 4Þ þ 0:5ðW subgrade À 8Þ À 2:5ðλ f À 1:2Þ According to the method for the calculation of the concave distortion of the freezing front in the subgrade area and for the determination of the key parameters proposed in this paper, the deformation of freezing front at any time can be quickly estimated. This can provide important data support for calculating subgrade deformation.

Conclusions
(1) Shortly after the subgrade is filled, the thermal effect of subgrade filling makes the freezing front under the center of subgrade concave distortion, which causes hydrothermal erosion in case of sufficient water inflow leading to continuous concave distortion. This in turn is likely to cause thawing settlement, cracking and other defects of the subgrade, thus proper waterproofing shall be done in the process.
(2) Compared with April and October, July presents as the most unfavorable engineering conditions, as the filling impedes heat output from soil to the atmosphere the most. In this scenario, the freezing front in the subgrade area displays the most severe concave distortion, and it takes the longest time for the temperature field of frozen soil to recover. Therefore, the filling should be carried out during the low temperature season as far as possible.
(3) The concave distortion of the freezing front in the subgrade area continues with the increase of the temperature difference between the filler and ground surface, subgrade height, subgrade width, and volumetric thermal capacity of filler, while decreases with the increase of the thermal conductivity of filler. Therefore, the filler used in practical engineering shall have been of small volumetric thermal capacity, low initial temperature and large thermal conductivity if possible.
(4) The paper proposes the method for the calculation of the concave distortion of the freezing front in the subgrade area and for the determination of the key parameters. The proposal helps conclude the concave distortion curves of the freezing front under different working conditions.