Numerical Simulation of Heat Transfer of Porous Rock Layers in Cold Sandy Regions

.


Introduction
Climate warming has led to the warming and thawing of permafrost at a global scale [1,2].To alleviate the road settlement caused by permafrost degradation, a PRL, which has excellent and stable heat transfer performance, is widely used to maintain the thermal stability of road embankments in cold regions, such as along the Qinghai-Tibet Railway (QTR) in China, Baikal-Amur Railway in Russia, Inuvik-Tuktoyaktuk and Yukon Highway in northern Canada, and Taylor Highway in Alaska, USA [3][4][5][6][7].However, climate change and frequent human activities have led to serious degradation and desertification of the grassland on the Qinghai-Tibet Plateau (QTP) [8].Strong wind conditions result in large amounts of sandy material deposition on the natural surface, ranging in thickness from a few centimeters to several meters [9].The construction of transportation infrastructure has altered the conditions of flow, transportation, and accumulation of near-surface aeolian sand.The embankment is buried due to the accumulation of aeolian sand [10].Further investigation reveals that there are approximately 450 km of hazardous sand sections in the Golmud-Lhasa section of the QTR, accounting for 25% of the entire section (Figure 1).Therefore, whether the PRL structure can maintain good heat-transfer performance even after being clogged by aeolian sand has always been a concern for engineers.
The PRL structure has excellent air convection and "thermal semiconductor" characteristics and has become the mainstream technology for maintaining the thermal stability of embankments in the permafrost zone of the QTP [3,7,11].In recent years, the heat transfer mechanism of PRLs has undergone comprehensive exploration using diverse technical methods, including field monitoring, laboratory model tests, and numerical simulations [12][13][14][15][16][17][18][19].As a highly permeable porous medium, a PRL has the characteristics of large pores and high pore-space connectivity.In winter, air density gradients caused by temperature changes at the bottom and top of the PRL can trigger gravity-driven internal air circulation, with cooler air flowing down to the bottom of the rock layer.In summer, the direction of the air density gradient in the porous medium is stable, and the cooler air is trapped at the bottom of the rock layer, which acts as an effective thermal insulator owing to the low thermal conductivity of the air [20,21].Although global warming has led to a rise in ground temperature, the PRL can still cool the subgrade below [22,23].However, in the areas with severe sand damage, a thick layer of aeolian sand usually accumulates on the slope and toe of the embankment, which can quickly fill the voids in the PRL [24].When the rough pores of the PRL are filled with aeolian sand, the connectivity of the high pore space is reduced, the thermal resistance of the variable equivalent thermal conductivity of the porous medium layer will cease to exist, and the convective heat transfer process is inhibited [25,26].The completely clogged PRL mainly depends on the heat conduction between solid particles, which leads to a significant reduction in the cooling effect of the PRL.Understanding the effect of aeolian sand filling on the heat transfer performance of closed PRLs is an important but unresolved issue in permafrost engineering.At present, there is a lack of scientific evidence to evaluate the impact of aeolian sand on convective heat transfer in PRLs.In this study, a two-dimensional finite-element model was used to simulate the impacts of climate warming and sand filling on the heat-transfer characteristics of a PRL and the thermal state of the underlying permafrost under different MAATs of −3.5, −4.5, and −5.5 °C.Our study aims to: (1) understand the heat transfer laws of a PRL in aeolian sand environments; (2) find out the necessary conditions under which a PRL can maintain subgrade stability under the dual effects of climate warming and aeolian sand; (3) provide a scientific basis for the evaluation of the long-term stability of a PRL applied to frozen-soil engineering in the aeolian sand environment of the QTP.

Governing Equations
In addressing permafrost engineering problems, it is often assumed that the rock layer can be considered infinitely long in the longitudinal direction [16].Therefore, the problem in this study can be reduced to a two-dimensional isotropic problem.Due to the different heat-transfer characteristics of different media, the model can be divided into porous media and soil layer zones.Based on the thermodynamic theories [27], air convection heat transfer within the PRL, coupled heat transfer for the PRL-soil system, and heat conduction and ice-water phase change in the soil layer are considered in the two-dimensional model.The sensible heat capacity method solves the phase-change problem in the soil layer.

Porous Media Zone
The PRL can be regarded as a high-permeability porous media zone with an unsteady non-isothermal percolation mode of thermal convection.The following two-dimensional differential equations describe the continuity, momentum, and energy transport in the study zones [27]. Continuity: where vx and vy are the velocity components of air (m•s − 1 ) in the directions x and y.Momentum: where ||   , µ is the kinematic viscosity of the air, ρ is the air density, P is the air pressure, k is the permeability of the porous medium, B is the Beta factor for non-Darcy flow (i.e., inertial resistance coefficient), and || is the inertial loss term.
The k and B in the PRL are calculated using the following equations: where dp is the effective average particle size of the medium, φ is the porosity of the medium, and α is a parameter related to the shape characteristics of the medium.
Since air is incompressible, its density  is a function of temperature and follows the Boussinesq approximation [21]: where β is the thermal expansion coefficient of air, and ρ0 and T0 are reference values for the density and temperature.Energy: where C * and λ * are the specific heat capacity and thermal conductivity.
Because the heat transfer between the PRL and soil is an unsteady process of phase change, according to the sensible heat capacity method, it is assumed that the phase change in each zone occurs within the temperature range of Tm ± ΔT.The C* and λ* are given as [28]: where Cf and Cu are the volumetric heat capacity of frozen and unfrozen states, respectively; λf and λu are the thermal conductivity of frozen and unfrozen states, respectively; L is the latent heat per unit volume.

Soil Layer Zone
Considering only the heat conduction of the soil skeleton and medium water as well as the phase change, and neglecting the thermal dissipation during the evaporation of the moisture in soil, the heat-transfer governing equation can be simplified to the following form: Since the PRL is in a closed state, each boundary is impermeable and can be regarded as a fixed boundary, the conditions are: The initial conditions are given as: where q, B, and n are the heat flux, solid boundary, and normal vector of each fixed boundary, respectively.Since the freeze-thaw problem is highly nonlinear, its analytical solution cannot be obtained.Therefore, numerical computational methods are used to solve the above governing equations [29,30].

Physical Model and Parameters
The physical domain chosen in this computational model consisted of a coarse rock layer with a grain size of 15 cm, aeolian sand, 2.5 m thick gravel soil, and 27 m thick strongly weathered mudstone.A 10 cm thick thermal insulation board was installed on both sides of the PRL to ensure insulation boundary conditions.The geometric dimension of the PRL was 1.5 m × 1.5 m × 1.5 m, corresponding to the field test.The strata were the same as those at the field test site [13].The porosity of the PRL was set to 0.4, and the corresponding permeability and inertial resistance factor were 1.58 × 10 −6 m 2 and 840.32 m −1 .The geometric model is shown in Figure 2. Parts I to IV were the aeolian sand-filling PRL, PRL, gravelly soil, and strongly weathered mudstone, respectively.The thermal conductivity of the sand-filling PRL is one of the important parameters in the temperature field calculation, which determines the rapidity of heat transfer between the ground and air and the change in ground temperature.At present, researchers have proposed many prediction models for the effective thermal conductivity of geotechnical materials based on mathematical and numerical modeling [31].The heat in the sandfilling PRL is transferred between solids (aeolian sand and rocks), pore water, and air, respectively.In this study, the thermal conductivity of the sand-filling PRL was calculated using a simplified model of aeolian sand developed by Lu et al. (2018) [32].The model gives a result between the upper limit obtained using the parallel flow model and the lower limit obtained using the tandem flow model.
In the unfrozen state: In the frozen state, assuming that the volume fraction of unfrozen water is zero, namely: where λs, λw, λi, and λa are the thermal conductivity of solid particles (aeolian sand and rock), water, ice, and air, respectively; A and φ are empirical parameters; xs, xw, xi, and xa denote the volume fraction of solid particles, water, ice, and air, respectively, obtained from Equation (15).
where ρd, ρs, and ρw are the density of dry sands, solid particles, and water.According to the calculation results, Table 1 summarizes the thermophysical parameters of each medium in the model.The study site is located at an altitude of 4500 m on the QTP, and the physical parameters of the air are shown in Table 2.The water content of the aeolian sand in the model was set to 3% based on field tests and laboratory experiments [33].There was no difference in the heat-transfer characteristics of the rock layer and insulation board in the frozen and unfrozen states; thus, the latent heat of freezing of these layers was also approximately zero.

Physical Domain cp (J•kg
Note: cp is the specific heat of air.

Temperature Boundary Condition
To represent the thermal interaction between the atmosphere and the ground, the Nfactor method was used to determine the empirical relation between the atmosphere and the ground surface [35].This empirical surrogate has been widely adopted to obtain temperatures for different ground surface conditions such as concrete, asphalt pavements, gravel surfaces (embankment slopes, shoulders), seasonal snow-cover, and vegetation coverage [36].A specified time-varying temperature boundary condition was set at the top surface of the PRL model.The temperature boundary conditions of the computational model were set as follows under the scenarios of MAATs of −3.5, −4.5, and −5.5 °C.
The natural ground surface (boundaries AB and CD) is expressed as: where Ta is the MAAT (°C).
As in some previous studies, climate warming is considered as the mean annual warming rate.The temperature boundary conditions of the numerical model calibration are available in reference [13], which is based on field observation data near the Hongliang River on the Tibetan plateau.The atmospheric warming due to climate change is represented by a linear growth rate of 0.052 °C/a [24].In addition, as there is no horizontal heat transfer in a semi-infinite body, the adiabatic boundary conditions are applied on lateral boundaries (AEG and DFH).A geothermal flux at the bottom boundary (GH) was applied as 0.04 W/m 2 .

Modeling Sequence and Calculation Cases
The numerical simulation was divided into two stages.In the first stage, the model used the natural ground surface temperature to generate a quasi-dynamic equilibrium of the subgrade.The initial value of each node of the natural ground was −1.5 °C.The output of the first stage was used as the initial temperature distribution for the second stage.In the second stage, the PRL area was manually activated after reaching the equilibrium state.The initial air temperature was set to the average temperature in summer considering that construction activities always occur in summer.Numerical simulations of the heat transfer of the closed PRL were performed over a simulation period of 50 years after the completion of the PRL for the MAAT scenarios of −3.5, −4.5, and −5.5 °C, respectively, subject to climate warming and sand filling.We selected 8 h as the time step.The time was accounted at the beginning of the second stage.To assess the effect of sand filling on the heat transfer of the PRL, five sand-filling thicknesses were chosen, 0, 20, 50, 80, and 120 cm.Furthermore, an additional case was executed to simulate the cooling effect of sand filling on the PRL in an un-warmed climate (Case 4).During the simulation period, the sand-filling layer was activated in the first year.The simulation cases are given in Table 3.The temperature data were recorded for each time compensation and used to generate the following simulation results.

Modeling Validation
The computational model used in this study was calibrated and validated based on the existing field test results [13].The field test area was selected near the QTR in the Hongliang River area of the QTP, which was heavily desertified.The PRL dimensions of the field test were consistent with those of the numerical model, and the sand-filling thickness was 50 cm.A thermistor was used to monitor the ground temperature, and the data were collected once a day.
Based on the numerical simulation results, the temperature and Rayleigh (Ra) number of the PRL during the period from June 2013 to June 2015 were calculated and compared with the monitoring data, as shown in Figure 3.The numerical model accurately simulated the thermal conductivity and convective heat transfer in a closed PRL for both sandfree (0 cm) and sand-filling of 50 cm cases.In addition, the monitoring results show that the convective heat transfer of the closed PRL only occurred in winter and the intensity of convection peaked in January, when the air temperature was at its lowest over the whole year.

Natural Convection Characteristics in the Closed PRL
A PRL is recognized as a porous medium material with high permeability.The Ra is an important index that determines the intensity of natural convection and pore air flow in porous media [22].When the following conditions are met, natural convection will occur in the PRL under the closed top: where C, β, and v are the volumetric heat capacity, thermal expansion factor, and kinematic viscosity, respectively; g is the gravitational acceleration; K is the medium permeability; H is the PRL thickness; ΔT is the temperature difference between the bottom and top of the PRL; and λ is the thermal conductivity of the medium.Particularly, the criterion for natural convection to occur in the PRL is Ra > 4π 2 [16].

Impact of Sand Filling on the Convection Characteristics of the PRL
To assess the impacts of sand filling on the cooling performance of closed PRLs, five cases (0, 20, 50, 80, and 120 cm) were carried out for study under the scenario of a MAAT of −3.5 °C. Figure 4a shows the time-varying ΔT of the closed PRL under the different sand-filling cases.The ΔT varies sinusoidally with time for each condition.The amplitude of ΔT fluctuations remained essentially unchanged over time.However, the mean annual ΔT decreases gradually with the increase in sand-filling thickness.In the first year of operation, the difference between the mean annual ΔT with 0 cm and 120 cm thick sand layers was 1.38 °C, and this difference gradually increased with the operation time.
According to Equation ( 18), the ΔT of the PRL determines the magnitude of the Ra number.When the ΔT of the PRL is negative, the Ra number tends to zero; conversely, when the ΔT is positive, the Ra number is positive.The corresponding temperature difference is the critical temperature difference (ΔTac) when Rac = 4π 2 .When Ra > 4π 2 , the Ra number gradually increases with the increase in the ΔT, and strong natural convection occurs inside the PRL, which enhances the convective heat transfer capacity of the bottom and top.Figure 4b presents the variation in the Ra number with time for different sandfilling cases.Without considering climate warming, the amplitude of the Ra number remained essentially constant with the increase in the operation time, but gradually decreased with the increase in the sand-filling thickness in PRLs.This is because the increase in the sand-filling thickness resulted in a decrease in the effective height of the porous media layer, forcing the ΔT of the PRL to decrease gradually, thereby increasing the difficulty of natural convection.In operation to the 10th year, the maximum values of the Ra number corresponding to the five cases were 123.38, 98.56, 68.52, 43.32, and 11.54, respectively.It can be seen that the Ra number gradually decreases with the increase in the thickness of the sand filling.Under the condition of the sand-filling thickness of 80 cm, only weak natural convection occurred within the closed PRL.When the sand-filling thickness exceeded 80 cm, the predicted Ra number was less than 4π 2 , and convection would not occur under this case.
In addition, the ΔTac and duration of natural convection were also different under each sand-filling case.Table 4 summarizes the ΔTac, natural convection period, and maximum Ra number of the four cases under different periods.The simulation results showed that the natural convection in the closed PRL was the strongest in cold seasons, which was consistent with the field monitoring results.Furthermore, the occurrence and end of natural convection were gradually delayed and advanced with the deepening of the sand thickness.In the sand-free case, the predicted duration of natural convection was usually approximately 120 days.When the sand-filling thickness increased to 80 cm, the natural convection duration decreased by approximately 75%, which was approximately 30 days on average.The velocity and direction of the air within the PRL can reflect its convective magnitude as well as a cooling effect.Research has confirmed that the direction of air movement inside the closed PRL is upward in the middle and downward on both sides, forming two vortices on the left and right [37].This is conducive to the dissipation of heat from the interior of the PRL and the introduction of cold energy from outside.In the cold season, the variation range of flow velocity within the PRL under the sand-filling 0, 20, 50, 80 and 120 cm cases was 1.56 × 10 −4 -2.63 × 10 −3 m•s −1 , 1.97 × 10 −4 -1.78 × 10 −3 m•s −1 , 8.47 × 10 −5 -7.56 × 10 −4 m•s −1 , 9.44 × 10 −6 -2.91 × 10 −4 m•s −1 , and 8.15 × 10 −6 -8.93 × 10 −5 m•s −1 , respectively.In the warm season, the variation range of flow velocity within the PRL was 6.14 × 10 −6 -1.73 × 10 −4 m•s −1 , 7.11 × 10 −6 -1.76 × 10 −4 m•s −1 , 3.54 × 10 −6 -1.43 × 10 −4 m•s −1 , 4.85 × 10 −7 -8.32 × 10 −5 m•s −1 , and 5.65 × 10 −7 -4.03 × 10 −5 m•s −1 , respectively.The simulation results showed that the air convection velocity within the PRL was larger in the cold season and smaller in the warm season, and thus the convective intensity was greater in the cold season than in the warm season.Meanwhile, driven by an unstable air density gradient attributed to the temperature difference between the top and bottom, the cooler air in the warm season was retained at the bottom of the PRL, which favored the lower underlying soil temperature.In addition, different sand-filling thicknesses influence the thermal state of the underlying permafrost by affecting the strength of air convection within the PRL.With the deepening of the sand-filling thicknesses, the air variation within the PRL decreased significantly, leading to a weakening of the natural convection strength and a reduction in the cooling performance.

Comprehensive Impact of Climate Warming and Sand Filling on the Convection Characteristics of the Closed PRL
To investigate the comprehensive impact of climate warming and sand filling on the convective heat transfer characteristics in the closed PRL, five case (0, 20, 50, 80, 120 cm) studies were carried out under three MAAT scenarios (MAAT=−3.5,−4.5, and −5.5 °C) at a warming rate of 0.052 °C per year.Figure 5 shows the variation in the Ra number over the 50-year simulation period under the three MAAT conditions affected by climate warming.It can be seen from the figure that the Ra number in each period was expected to gradually decrease as the sand-filling thickness increased.In the first year, the maximum Ra number corresponding to the three MAATs decreased by 81.32, 81.75, and 85.13, respectively, as the sand-filling thickness increased from 0 cm to 80 cm.Furthermore, the maximum Ra number also decreased gradually with the increase in the operation year.By the 50th year, the maximum Ra number for the three MAATs in the sand-free case decreased by 14%, 5.6%, and 2.9%, respectively.The maximum Ra number decreased by 11.1%, 4.8%, and 3.3% for the three MAATs in the 80 cm sand-filling case, respectively.The results show that the lower the MAAT, the smaller the decrease in the Ra number with the increase in operation years.Also, the mean annual ΔT exhibited similar variation trends with the Ra number.In the case of sand-free, the mean annual ΔT corresponding to the three MAATs decreased from −2.32 °C, −1.81 °C, and −2.24 °C in the 1st year to −3.01 °C, −2.73 °C, and −2.62 °C in the 50th year.Simultaneously, the amplitude of ΔT in each period also decreased slightly, sloping from 11.54 °C, 11.82 °C, and 11.12 °C in the first year to 9.36 °C, 11.26 °C, and 10.41 °C in the 50th year.
The duration of natural convection occurring in the PRL also changed dramatically due to the combined impacts of aeolian sand and climate warming.When the thickness of the sand filling was less than 50 cm, there was strong natural convection in the PRL throughout the simulation span.However, due to the impact of climate change, the duration of natural convection gradually decreased with the increase in clogging thickness and time of operation.Especially in the case of the sand-filling thickness reaching 80 cm, the Ra number gradually decreased over time, and the PRL finally lost its cooling effect when the Ra number was less than Rac.The prediction results indicated that the cooling effect of the closed PRL disappeared after the 30 th and 40 th years under scenarios with MAATs of −3.5 °C and −4.5 °C, respectively.In contrast, for a MAAT = −5.5 °C, although natural convection could still occur throughout the simulation period, the duration of natural convection per year was reduced from 27 days in the 1 st year to 11 days in the 50 th year.Furthermore, for the case of 120 cm sand-filling thickness under three MAATs, the ΔT was far less than the ΔTac that allowed natural convection to occur, and no natural convection occurred during the 50-year simulation period.

Variation in the Permafrost Table
Figures 6-8 show the changes in the permafrost table against time for two cases of sand-free and sand-filling of 80 cm in different MAAT scenarios.The predicted results showed that the permafrost table for both cases reflected a gradual downward trend against time due to the influence of climate warming and sand filling.In the first 20 years, the permafrost table beneath the PRL was higher than the initial permafrost table, which indicated that the PRL cooled the subgrade during this period and could weaken the degradation of permafrost caused by climate warming.Under the scenario of MAAT = −5.5 °C, in the case of sand-free, the permafrost table in the 10 th year was 0.56 m higher than the initial permafrost table, with the artificial permafrost table entering the PRL.During the subsequent simulation period, the permafrost table gradually increased and continued to penetrate deep into the subgrade.By the 30 th year of operation, the artificial permafrost table under both cases was already lower than the initial permafrost table, indicating that the cooling performance of rock layers could not resist the impact of climate warming on permafrost.Since then, the permafrost subgrade was in a continuous degradation state, and the difference between the artificial permafrost table and the natural permafrost table beneath the PRL gradually decreased due to climate warming.At a MAAT of −3.5 °C, the differences between the cases of sand-free and sand-filling thickness of 80 cm were 0.49 m and 0.33 m.At MAATs of −4.5 °C and −5.5 °C, the differences between both cases were 0.64 m and 0.41 m, and 0.76 m and 0.5 m, respectively.However, this difference will be further narrowed in the next 20 years.Between the 30 th and the 50 th year, for the case of sand-free, the reduction rates of the artificial permafrost table corresponding to the three MAATs were 0.049 m•a −1 , 0.037 m•a −1 , and 0.035 m•a −1 , respectively.Under the case of 80 cm sand filling, the reduction rates of the artificial permafrost table were 0.042 m•a −1 , 0.026 m•a −1 , and 0.025 m•a −1 , respectively.In addition, the permafrost table of the sand-free case was lesser than that of the sand-filling cases, indicating that sand filling has a negative impact on the cooling performance of the PRL.Overall, with the decrease in MAAT, the impact of aeolian sand accumulation on the permafrost table in the lower part of the PRL was weakened, the cooling performance of the PRL became significant, and the duration in which the artificial permafrost table was higher than the natural permafrost table was prolonged.Therefore, closed PRLs are more suitable for cooling reinforcement measures of the subgrade in permafrost regions with low temperatures to maximize the cooling performance of natural convective heat transfer.3.3.2.Variation in Heat Flux of the Shallow Soil Layer Beneath the PRL Implementing the PRL alters the heat exchange between the atmosphere and the subgrade, thus increasing the heat extraction of the underlying permafrost.To accurately estimate the cooling effect of the PRL on the thermal state of the soil layer, the heat fluxes at depths of 1.5 to 2.0 m beneath the PRL are estimated using Equation ( 19) according to heat transfer theory [12].
where λ is the thermal conductivity of the soil, Δz is the soil layer thickness, and T1.5 and T2.0 are the soil temperatures at the depths of 1.5 m and 2.0 m.The mean annual heat budget for the subgrade was obtained by integrating the heat flux over time.The thermal budget of the soil layer under different operating cases is shown in Figure 9.
The annual heat budget of the soil layer was mainly affected by the thickness of the sand filling.Under the three MAAT conditions, the annual heat budget was negative in the cases of sand-filling thicknesses of less than or equal to 50 cm, indicating that the soil layer was always in a release heat state, and the PRL still played a cooling effect during the operating period.However, the annual heat budget tended to increase gradually over time due to climate warming.In the MAAT of −3.5 °C scenario, the thickness of sand filling reached 80 cm, and the thermal budget turned from negative to positive, with the heat absorption gradually increasing with time.In the MAAT of −4.5 °C scenario, the underlying permafrost was exothermic in the first 10 years under the case of 80 cm sand-filling thickness, the heat budget became positive in the 20th year, and the heat budget increased during the subsequent operating time, reaching 439.43 kJ in the 50th year.In the MAAT of −5.5 °C scenario, due to the low MAAT, the exothermic state of the underlying permafrost in the 80 cm sand-filling case lasted longer, with the heat budget becoming positive after the 30th year.In addition, in the 120 cm sand-filling case, the heat budgets corresponding to the three MAATs were all positive, indicating that the active cooling performance of the PRL failed to play a role.Simultaneously, the continuous heat accumulation in the shallow layer inevitably leads to a temperature increase in the deep permafrost layer, which promotes the degradation of permafrost.

Thermal Changes in the Deep Soil
Figure 10 shows the mean annual ground temperature (MAGT) at a depth of 15 m beneath the centerline of the PRL over the simulation periods.Due to climate warming, the MAGT at 15 m for the three MAAT scenarios increased significantly with time.It should be noted that the MAGT of the deep layer with sand filling was higher than that in the sand-free case.After 30 years of operation, the temperature differences at 15 m in the two cases (sand-free and 80 cm sand filling) were 0.014 °C, 0.089 °C, and 0.036 °C, respectively.In addition, the cooling effect of the PRL in the first 10 years could weaken the permafrost warming caused by climate warming, and essentially maintained the thermal state of the deep permafrost beneath the PRL under the warming rate of 0.052 °C•a −1 .In the first 10 years, the warming rates of MAGT for the two cases were only 0.013 °C•a −1 and 0.018 °C•a −1 .During subsequent simulation periods, the MAGT at 15 m exhibited a more pronounced warming trend.From the 20th to 50th year, in the sand-free case, the corresponding temperature rise amplitudes of the three MAATs were 0.174 °C, 0.899 °C, and 1.298 °C, and with warming rates of 0.006 °C•a −1 , 0.03 °C•a −1 , and 0.043 °C•a −1 , respectively.In the case of the sand-filling of 80 cm, the three MAAT scenarios had temperature rise amplitudes of 0.166 °C, 0.844 °C, and 1.306 °C, with warming rates of 0.005 °C•a −1 , 0.028 °C•a −1 , and 0.044 °C•a −1 , respectively.Overall, the changes in deep ground temperature under the PRL were mainly affected by climate warming and MAAT, while the aeolian sand had an insignificant effect on the deep ground temperature.

Discussion
At present, the application of a PRL in the QTR encounters challenges attributed to climate warming and sand damage.In response to this issue, this paper evaluated the cooling performance of a closed PRL under a cold sandy environment based on convective heat-transfer characteristics, permafrost table, heat budget, and MAAT.The temperature boundary conditions of the model were determined based on measured data in the Hongliang River region of the QTP.Owing to the differences in the climatic conditions (warming rate, solar radiation, and rainfall) and geological conditions (vegetation coverage, soil type, MAGT, and ice content) in different regions, the laws of the adherent layer are also different [36,[38][39][40].Furthermore, the parameters of air in the PRL assumed in this study are constant, but they are moderately affected by changes in air pressure, temperature, and humidity.Therefore, they should be assigned different values at different MAGTs and altitudes.Meanwhile, most studies use a linear ground surface temperature increase rate to evaluate the effects of climate change induced temperature increase on permafrost warming and thawing [29,41].Therefore, the predicted values of climate warming rates used in this study do not affect the regularity of the results.
This study has confirmed that climate warming, MAAT, and sand filling have a significant impact on the convective heat-transfer characteristics of the closed PRL.Engineering activities will disturb the original surface energy balance, and the construction of the PRL will raise the permafrost table to a certain extent.We found that the rate of change of the permafrost table beneath the closed PRL was different for three MAAT scenarios.With a decrease in MAAT, the effect of sand filling on the permafrost table beneath the PRL was weakened, the difference between the artificial permafrost table and natural permafrost table increased, and the cooling performance of the closed PRL became significant.Under the MAAT of −5.5 °C scenario, in the first 10 years of operation, the artificial permafrost table entered the rock layer (Figure 8).As indicated, a closed PRL is more suitable for cooling reinforcement measures of the subgrade in the low-temperature permafrost region to maximize the cooling performance of natural convection heat transfer.In addition, aeolian sand accumulation and clogging greatly alters the surface albedo, where the moisture content of the aeolian sand affects the amount of heat entering the subgrade.Field surveys show that the moisture content of aeolian sand varies from 0.1 to 15% in different areas of the QTP, and its thermal conductivity varies from 0.3 to 2.4 W•m −1 •°C −1 [8,32].The dry sand layer can play the role of heat insulation, hinder the heat transfer between the atmosphere and the subgrade, and promote the warming of the underlying permafrost.On the contrary, the high thermal conductivity of wet sand is conducive to more heat extraction in winter.This shows that the thermophysical properties of aeolian sand play an essential role in the thermal state of permafrost beneath the PRL.Therefore, it is necessary to establish a surface energy balance model and hydro-thermal coupling model to study the serviceability of a PRL with serious sand damage in future research.

Conclusions
A coupled heat transfer model for a soil-PRL system was developed.The heat transfer characteristics of PRLs in cold sandy environments were analyzed via numerical simulation and revealed the necessary conditions under which a PRL can maintain subgrade stability under the dual effects of climate warming and aeolian sand.The main conclusions are as follows.
(1) The accuracy of the numerical model was verified using field tests.Natural convection within the closed PRL occurred only in cold seasons, and the convection strength was related to the effective convection height of the rock layer.As the thickness of sand filling increased, the Tac allowing natural convection to occur increased, and the Ra number decreased, which caused the weakening of the duration and intensity of natural convection.
(2) Under a warming scenario of 0.052 °C•a −1 , the cooling performance of a PRL can offset the adverse impacts of climate warming and raise the permafrost table during the first 20 years of operation.However, the cooling performance of the PRL diminishes with the increase in the operation year, and the underlying permafrost continues to degrade over the next several decades.
(3) A closed PRL is more suitable for cooling measures of the subgrade in permafrost regions with colder MAATs.In the context of climate change and sand damage, the cooling effect of a PRL on the permafrost can no longer meet the long-term requirements.

Figure 3 .
Figure 3. Model calibration using field test results from June 2013 to June 2015: (a) Temperature variation at depths of 0.5 m and 1.5 m in the PRL, (b) Ra numbers.

Figure 4 .
Figure 4.The time-varying ΔT and Ra number of the closed PRL under the different sand-filling cases: (a) Temperature difference between the bottom and top of the PRL, (b) Ra number.

Figure 9 .
Figure 9. Heat budget of the soil layer between depths of 1.5 and 2.0 m beneath the PRL: (a) MAAT = −3.5 °C, (b) MAAT = −4.5 °C, (c) MAAT = −5.5 °C.Note: when q > 0, it means heat absorption from the shallow layer; conversely, q < 0 means heat release from the deep layer of permafrost to the shallow layer.

Figure 10 .
Figure 10.Variations in MAGT at a depth of 15 m beneath the natural ground surface over a period of 50 years under three MAAT (−3.5 °C, −4.5 °C, −5.5 °C) scenarios.

Table 1 .
Thermophysical parameters of the media in model.

Table 3 .
List of simulated cases.

Table 4 .
Critical temperature difference, duration, and maximum Ra number of the closed PRL in different periods.