# Thermal-Moisture Dynamics at Different Underlying Surfaces in Permafrost Regions of the Central Tibetan Plateau by considering the Effect of Rainfall

Authors: Bingbing Lei; Bin Wang; Mingli Zhang; Zhixiong Zhou; Guang Li; Guodong Yue
Journal: Geofluids
 (2022)

---

## Abstract

The Qinghai-Tibetan Plateau (QTP) has undergone an increase in rainfall and a drastic alteration in the moisture-heat regime in active layers and engineering. To investigate the water and heat responses of natural ground and engineering to rainfall, the differences in energy on the ground surface and the thermal-moisture dynamics of different permafrost underlying surfaces were discussed. Based on the meteorological data in 2013 observed at the Beiluhe observation station, three types of underlying surfaces (i.e., natural ground, asphalt pavement, and gravel pavement) were selected to compare the differences in energy balance at the ground surface, water-energy transport process, and coupling mechanism in active layers under rainfall conditions by a coupled water vapor-heat model of the unsaturated frozen soil. The results show that the asphalt pavement greatly increases the surface net radiation and soil surface heat flux, decreases the surface evaporation latent heat, and cuts off the moisture migration between the atmosphere and the active layer. The gravel pavement significantly increases the surface evaporation latent heat to lower the soil surface heat flux, and the amount of moisture in shallow soil is strongly influenced by rainfall and evaporation. Therefore, the moisture migration and accumulation under the asphalt pavement are dominated by the water vapor flux under thermal gradients, whereas the liquid water under the water potential gradients is the major source of moisture migration under the gravel pavement. The heat transfer in the shallow active layer is dominated by heat conduction. The effect of heat conduction, water vapor migration, and phase transition on the soil temperature is evident for the asphalt pavement, while the impact of liquid water migration on the shallow soil temperature for the gravel pavement is significant in the thawing period. As a result, the soil temperature relationship between different underlying surfaces is asphalt pavement>gravel pavement>natural ground. The thickness of the active layer gradually decreases. Although rainfall infiltration promotes the liquid water convection of the gravel pavement, the decrease in heat flux is less than the increase in thermal conductivity. In general, the construction of asphalt pavement and gravel pavement accelerates the degradation of permafrost. The results can provide theoretical and simulated guidance for the stability prediction and analysis of various underlying surfaces in the central QTP where rainfall is increasing.


---

## Body

## 1. Introduction

The formation of permafrost is attributed to the long-term exchange of moisture-heat between the active layers and the atmosphere [1]. The hydrothermal properties and distribution of permafrost are influenced by climate changes and human activities [2, 3]. Climatic changes strongly affect the active layer, which is a common underlying surface widely distributed in permafrost regions of the QTP [4]. Changes in underlying surfaces from engineering affect the exchange process of energy and moisture at the boundary and hence the original energy and mass balance in active layers [5, 6]. The hydrological and ecological environment of permafrost and the freeze-thaw process, as well as the permafrost stability and engineering, are disrupted by the variation in active layers [7, 8]. Rainfall, an important supplemental water in the QTP, significantly influences the energy budget of underlying surfaces and the hydrothermal state in shallow active layers [4, 9]. Under the climate background of increasing rainfall [10, 11], it is vital to realize the thermal-moisture dynamics of different underlying surfaces under the effect of rainfall.The distribution of moisture-heat in permafrost under different underlying surfaces has been widely investigated, including different natural surfaces and engineering surfaces. Based on the monitoring of different land surfaces in the QTP, it was found that the difference in water retention capacity and thermal conductivity in the surface layer of soil results in different responses of the active layer to soil temperature [12]. The coarse-grained soil affects the hydrothermal properties of the active layer by controlling the soil water retention capacity and thermal conductivity [10]. The organic soil, with its low heat conduction and high heat capacity properties, reduces permafrost perturbation by decreasing the reaction of soil thermal conditions to solar radiation and ambient temperature [13]. Compared with engineering surfaces and natural ground, the buildup of heat in the active layer of asphalt pavement with strong heat absorption capacity (low albedo) causes the permafrost to melt [14]. The research on hard pavement found that the accumulation of moisture in the active layer caused by changes in the temperature field increases the accumulation of soil heat, which is detrimental to the stabilization of permafrost [6, 15]. Moreover, it was found that the heat dissipation in active layers is affected by the air humidity, temperature, and heat capacity by the observation of the near-surface water-thermal state [5]. The hydraulic properties of the pavement structure significantly change the infiltration and evaporation of water and thus affect the hydrothermal distribution of the active layer [16]. The current knowledge about the moisture migration and energy transport in various permafrost underlying surfaces is inadequate due to the complexity of hydrothermal characteristics and coupling mechanisms in the freeze-thaw process. The surface thermal-moisture exchange and soil water and heat transport mechanisms must be studied to fully comprehend the issue.The hydrothermal distribution of active layers is intimately associated with the moisture-heat exchange of the ground-atmosphere boundary. In the research on the thermal boundary, researchers proposed an empirical and statistical approach, a direct regression of the observed surface temperature [17]. Other methods were also used to estimate the ground surface temperature, such as the N-factors method [12], the adherent layer theory [18], and the thermal orbit regression method [19]. The use of the above method still requires abundant on-site observed data to improve its reliability and accuracy in the study of boundaries. The differences and changes in sensible heat, latent heat, moisture content, and solar radiation on different underlying surfaces are evident in fact [9, 15]. Therefore, the boundary moisture-heat exchange process of different underlying surfaces is considered by the principle of ground surface energy balance. The characteristics of the moisture boundary are more complicated than those of the thermal boundary because of the shallow active layers that are strongly influenced by the environment. According to some researchers, the moisture boundary is a saturated water condition or a boundary determined by the freezing characteristic curve under the water budget constraints [6, 20]. The water migration law in partially frozen soil was investigated using numerical simulation with a constant water flux boundary [21]. The flux affected by the boundary pressure head difference was considered in the study of the rainfall infiltration mechanism of the seasonally frozen soil slope [22]. The above methods are based on a series of assumptions about the characteristics of research objects. It is unreasonable to ignore the evaporation in the central QTP where the single rainfall is small and the rainfall frequency is high. Therefore, it is essential to consider the effect of rainfall and evaporation on the studies of the moisture boundary of different underlying surfaces.The thermal-moisture dynamics of different underlying surfaces can be calculated and analyzed using numerical models that integrate theoretical approaches and measured values. Chen et al. studied the temperature difference of frozen soil between the sand-filled and sand-covered gravel layers under climate warming conditions through numerical simulation [23]. The freeze-thaw process of permafrost under snow cover and organic matter soil cover was studied by a numerical model driven by meteorological conditions [13]. The numerical model, CoupModel, was used to find that the vegetation cover and soil properties are the main reasons for the difference in the water and thermal transfer in permafrost [24]. The hydrothermal response characteristics of the mixed layer of sand and rock under rainfall conditions were numerically studied [25]. However, in the above computational models, the mechanism of water transfer in soil and the moisture-heat exchange at the boundary are not fully considered. Virtually, the liquid and vaporous water migration in soil is obvious under large soil temperature difference conditions in the QTP [15, 26]. The water budget of the active layer, which is affected by rainfall and evaporation, alters the energy balance and the hydrologic and thermal properties of shallow soil [4, 9]. Furthermore, the changes in soil hydraulic and thermal parameters during the freeze-thaw process are vital for determining numerical simulation outcomes [4, 12, 27].In view of the above, the aerodynamic approach was used to estimate the energy balance of the ground surface between the frozen soil and the atmosphere, which is regarded as the thermal boundary condition. The role of rainfall and evaporation was considered in the moisture boundary theory. A numerical model established by our team that has been well validated and applied was used in the study [4, 27, 28]. The typical underlying surfaces, such as asphalt pavement (AP), gravel pavement (GP), and natural ground (NG) of an alpine meadow in the central QTP, were taken as research objects. The meteorological data were introduced into the coupled model to calculate the differences in the ground surface energy components, the hydrothermal dynamics in the active layers, and the permafrost table of the three underlying surfaces. The results can provide theoretical guidance for the subgrade stability and the stability of permafrost under different underlying surfaces under the background of humidification.


## 2. Methods

### 2.1. Mathematical-Physical Model

#### 2.1.1. Hydrothermal Coupling Model

The liquid water flow in unsaturated soil has a positive correlation with hydraulic conductivity. It is governed by the potential energy of soil matrix, temperature, and gravity. The process obeys the Richards function [29], namely,
(1)ql=qlh+qlT=−Klh∇h+y−KlT∇T,where ql is the liquid water fluxes (m/s).qlh and qlT are the liquid water fluxes under the effect of soil matrix and temperature potential (m/s), respectively. Klh and KlT are the liquid water hydraulic conductivities of hydraulic pressure (m/s) and temperature potential (m2/(K·s)), respectively. h is the hydraulic pressure (m), and y is the gravitational potential (m). T is the temperature in soil (K).In contrast to the mechanism of liquid water flow, it is noted that the vaporous water transport in unsaturated soil is affected by other potential energy except the gravity potential. The process obeys Fick’s law [30]:
(2)qv=qvh+qvT=−Kvh∇h−KvT∇T,where qv is the vaporous water fluxes (m/s).qvh and qvT are the vaporous water fluxes of soil matrix and temperature potential (m/s), respectively. Kvh and KvT are the vapor hydraulic conductivities of hydraulic pressure (m/s) and temperature potential (m2/(K·s)).According to the mass conservation law, the fluid flow in and out of the soil per unit of time keeps the balance. The liquid and vaporous water during the soil freeze-thaw process can be described as follows [4, 27]:
(3)∂θ∂t=∇∙Klhh+y+KlT∇T+∇∙Kvh∇h+KvT∇T−ρvρw∂∂tn−θl−θi,(4)θ=θl+ρiρwθi,where θ is the total volumetric content of water and t is time (s). ρw, ρv, and ρi are the densities of liquid, vaporous, and solid phases of water in soil (kg/m3). n is the soil porosity. θl and θi are the volume contents of liquid and vaporous phases of water (m3/m3). The detailed hydraulic properties can be referred to in the research by Zhang et al. [27, 28]. The right third term in Equation (3) is the variation of equivalent vapor.The following energy conservation equation provides the governing equation for energy transport in a variably saturated, stiff porous material [4, 27]:
(5)Cp∂T∂t=∇·λp∇T−ql∇ClT−qv∇CvT−Lwρw∇qv+Lfρi∂θi∂t−Lwρw∂θv∂t,where λp is the heat conduction coefficient of soil (W/(m·K)). Cl, Cv, and Cp are the volumetric thermal capacity of liquid water, vaporous water, and soil (J/(kg·K)), respectively. Lw and Lf are the latent heat of soil water evaporation and ice-water phase transition (J/kg), respectively. θv is the volumetric content of vaporous water (m3/m3). The left term in Equation (5) is the variation of temperature per unit of time. The right terms in Equation (5) are the heat conduction, the convective heat of the liquid and vaporous water, the vapor diffusion heat, and the phase transition heat of ice-liquid and liquid-vapor.The volumetric thermal capacity and heat conduction coefficient of frozen soil are comprehensive parameters affected by thermal parameters and the contents of soil components. It can be expressed using the following formulas [31]:
(6)Cp=Cs1−n+Clθl+Ciθi,(7)λp=λs1−nλlθlλiθiλaθv,where Ci and Cs are the volumetric thermal capacity of the ice and soil matrix, respectively (J/(m3·K)). λs, λl, λi, and λa are the heat conduction coefficients of the soil matrix, liquid, solid, and vaporous water (W/(m·K)).The volumetric liquid water and ice content in the soil during the freeze-thaw process can be calculated by the relationship between soil temperature and equivalent volumetric water content [31], namely,
(8)θumax=a×T−273.15−b,(9)θl=θumax,T<273.15andθl>θumax,θ,otherwise,(10)θi=ρwρiθ−θl,where a and b are two soil characteristic parameters [31]. θumax is the maximum unfrozen water content of frozen soil below the freezing point (m3/m3).The V-G model [32] was applied to establish the relationship between the matrix potential and the liquid water in unsaturated soil:
(11)Se=θl−θrθs−θr=11+aψhωm,h<0,1,h≥0,(12)Klh=KsSel1−1−Se1/mm2∕I,h<0,Ks,h≥0,(13)KlT=KlhhGWT1γ0dγdT,where Se is the soil saturation. θs and θr are the liquid water contents in the saturated and residual states (m3/m3), respectively. aψ is a soil intake reciprocal value. ω and m (=1−1/ω) are empirical parameters. Ks is the hydraulic conductivity in a saturated state (m/s). I is the ice resistance coefficient. GWT is a gain factor, and γ is the soil water surface tension (J/m2), γ0 is the surface tension at 25°C (J/m2).The fundamental equation used to calculate the diffusion of water vapor in soil is Fick’s diffusion law. The diffusion coefficients of water vapor under the water potential gradient of soil and temperature gradient are expressed as follows [33]:
(14)Kvh=DaρwρvsΗrΜgRT,(15)KvT=DaρwηΗrdρvsdT,(16)η=9.5+3θvθs−8.5exp−1+2.6fcθl+θiθs4,where Da is the vapor diffusivity in soil (m2/s). ρvs is the density of saturated water vapor (kg/m3). Ηr is the relative humidity. Μ is the molecular weight of water (kg/mol). g is the gravitational acceleration (m/s2). R is the universal gas constant (J/(mol·K)). η is the enhancement factor of water vapor [34]. fc is the proportion of clay in the soil.


#### 2.1.2. Hydrothermal Flux Boundary

The aerodynamic principle at the near-surface is an effective method to characterize the thermal-moisture exchange process between the underlying surface and the atmosphere caused by meteorological conditions such as radiation, rainfall, and evaporation. The ground surface energy balance equation is as follows [35]:
(17)G=Rn−Hs−LwE,where G is the soil surface heat flux (W/m2). Rn, Hs, and LwE are the fluxes of net radiation, sensible heat, and latent heat on the ground surface (W/m2), respectively. The Rn can be written as follows [36]:
(18)Rn=Rns+Rnl=1−αRns↓+σTa40.605+1.77ha−σεsTs41063600,where Rns and Rnl are the surface net shortwave and net longwave radiation (W/m2), respectively. Rns↓ is the solar shortwave radiation (W/m2). α is the albedo on the ground. σ is the Stefan-Boltzmann constant (J/(m2·s·K4)). Ta and Ts are the temperatures of air and ground, respectively (K). ha is the effective volumetric water content of air (kg/m3). εs is the surface emissivity.Soil color has a significant impact on the ground albedo. The pure black soil can absorb almost all visible solar radiation. Its albedo is the smallest. Inversely, the albedo of pure white soil is the largest. The water around the soil particles increases the absorption route of sunlight, which makes the surface reflectivity decrease gradually with the increase of soil moisture. The ground albedo is defined as a function of the surface type and water content [37]:
(19)α=αs+∆αg,(20)∆αg=0.11−0.4θl,where αs is the saturated soil albedo and ∆αg is the increment of albedo induced by dehumidification and drying of the soil. The ground albedo of different soil colors and moisture is in Table 1 [38].Table 1
Surface albedo of different surface types and water contents.
ParameterWavelengthColor (from white (1) to black (8))12345678Saturated soil<0.7μm0.120.110.100.090.080.070.060.05≥0.7μm0.240.220.200.180.160.140.120.10Dry soil<0.7μm0.230.220.200.180.160.140.120.10≥0.7μm0.460.440.400.360.320.280.240.20The surface evaporation water flux (kg/(m2·s)) is as follows:
(21)E=ρvss−ρvara+rs,where ρvss and ρva are the vapor densities in the near-surface soil and the atmosphere, respectively (kg/m3). ra is the aerodynamic resistance to heat transfer (s/m). rs is the surface soil resistance to vapor (s/m).The surface sensible heating flux is as follows:(22)Hs=CaTs−Tara,where Ca is the volumetric thermal capacity of the near-surface air (J/(m3·K)). The G obtained from formulas (17)–(22) was used as the upper thermal boundary.In the Beiluhe area, there is relatively little snowfall. The influence of snowfall on the ground can be overlooked [27, 28]. The soil surface water flux (m/s) can be calculated by the E, ρw, and precipitation rate (Prain, m/s):
(23)Q=Eρw−Prain.


### 2.2. Model Implementation and Validation

The above mathematical-physical model formed by some high-order nonlinear differential equations was solved by the COMSOL Multiphysics 4.3a finite element software. COMSOL has powerful multiphysics modeling capabilities, making it easy for users to use the various built-in physics modules. Moreover, users can insert partial differential equations (PDEs) and set how to establish the connection between multiphysics equations for the realization of sophisticated physical processes. The governing equations in the model are Equations (3) and (5), with temperature and total volumetric water content as the essential variables, which can be calculated by the PDE tool of COMSOL. The standard form of a generalized PDE is as follows:
(24)ea∂2U∂t2+da∂U∂t+∇∙Γ=f,where ea is a damping coefficient. da is a mass coefficient. Γ is a flux vector, and f is a source term. We can input the governing equations of the above mathematical model with reference to the standard form of PDE, where U is a variable. The θ and T can be input, and the remaining items can be input parameters, partial derivatives, etc. The above equations were solved by invoking the COMSOL transient solver, which solves a possibly nonlinear system of equations at each time step through a set of iterative techniques based on Newton’s method. The iterative technique of Newton’s method allows the function and its Jacobian matrix (derivative) to be computed at each time step. If the solver is poorly converged, it will reduce the requested time step and try to compute the result. The iterative calculation is repeated until a relatively accurate result is obtained within the set maximum number of Newton’s iterations. Therefore, we can improve the model convergence and get an accurate result through the internal setting of the nonlinear solver in advance. The steps of the backward difference method (BDF) in the solver were set to “Strict,” which helps to get more accurate solutions within a finite time step. Discretization in time varies between a minimum time step of 60 s and a maximum time step of 900 s. The allowable relative tolerance is 10-3. Meanwhile, a smaller mesh of spatial discretization, namely, 1 cm, was used.The reliability of the aforesaid numerical model was confirmed by comparing the expected and measured outcomes. MBE (mean bias error) andd (agreement index) were used as the key evaluation indicators. The simulated temperatures and volumetric water contents were compared to the measured values at the test site from October 1, 2012, to October 1, 2013. As shown in Table 2, the processed statistics indicate that the verification results of the temperature at 5, 15, 25, 35, 75, 150, 160, 170, 200, and 220 cm depths are d>0.802 and MBE<1.204°C. Meanwhile, the result is d>0.853 and MBE<1.353% for the water content [4, 28]. A deviation of 2~3°C from the measured soil temperature is regarded as acceptable [35]. The d index of the estimated value steadily approaches 1.0, which is more favorable for enhancing the matching degree between the anticipated and measured values [35]. The simulated results correlate well with the observed findings, and the error control is good. It indicates that the model is applicable to computing the moisture-heat condition in the shallow soil of various underlying surfaces.Table 2
Summary statistics for differences between the measured and simulated values at different depths [39].
Depth (cm)Soil temperatureLiquid water contentsd (-)MBE (°C)d (-)MBE (%)50.9620.6160.9430.267150.9680.7560.9770.051250.9730.3930.8551.353350.976-0.0520.9991.031750.9570.4180.8530.7181500.917-0.5730.8910.8081600.881-1.0360.9160.4131700.826-1.2040.9020.3112000.8020.9580.9780.1752200.862-0.7550.8931.288


### 2.3. Model Application

The model has been used in different application scenarios and engineering environments. The reliability is high [27, 28]. The NG, AP, and GP are some of the typical underlying surfaces in the central QTP. It is of great engineering and practical significance to carry out extensive research. Therefore, we selected several underlying surfaces near the frozen soil observation station of Beiluhe as the research objects, which allows us to use the same meteorological factors as the model-driven data in the numerical models of different underlying surfaces.


#### 2.3.1. Research Location

The place of research lies in the permafrost test section in the Beiluhe Basin in the central QTP (34.845°N, 92.940°E, 4620 m a.s.l.) (Figure 1). The NG is an alpine meadow with a high content of silty clay and less vegetation cover [4, 15]. In the test section, there are two kinds of artificially covered underlying surfaces: the AP and GP. The soil hydraulic and thermal parameters on the different underlying surfaces were determined by the field-tested data. The upper limit of permafrost is 2 to 3 m, and the annual mean ground temperature is about -1.8 to -0.4°C [28]. The average air temperature for the time window between 2004 and 2013 was −3.4°C [4]. In 2013, the yearly average air temperature was -3.31°C, with the maximum and lowest monthly temperatures of 8.86°C (August) and -14.14°C (January), respectively. The analysis of meteorological data reveals that the rainfall is characterized by multiple times and small amounts, and the rainfall mainly occurs in May-September of each year. There are fewer snow-covered events on the ground surface in cold seasons, generally no more than a week [27, 28].Figure 1
Schematic location of the study site (from LocaSpaceViewer).


#### 2.3.2. Model Setup and Forcings

The three kinds of typical underlying surfaces (i.e., AP, GP, and NG) near the frozen soil observation station of Beiluhe were selected as the research objects. A model with a depth of 15 m was established. The downward coordinate is negative (Figure2). It was reduced to a one-dimensional model since the test location is flat and there is no visible hydrothermal lateral advection. For these three different underlying surface models, the soils within 1 m below the ground surface were set as GP (0.30 m gravel packing and 0.70 m subsand), AP (0.10 m asphalt surface layer, 0.20 m sandy gravel, and 0.70 m clayey silt), and NG (1.0 m clayey silt), respectively. The properties of soils below 1.0 m depth are the same. The thermal-moisture parameters required for the calculation of the coupled model are shown in Table 3. The data in Ref. [31] was obtained from a large number of experimental results on different soils in different research sites, including the QTP.Figure 2
Vertical soil-layer distributions at the research location.Table 3
Soil-layer thermal-moisture parameters of the underlying surface [31].
Soil typenλs (W/(m·K))Cs (J/(m3·K))abaψ (1/m)ωθr (m3/m3)θs (m3/m3)Ks (m/s)fcAsphalt surface0.151.595014.410.56021.562.000.010.072.3E-110.00Gravel surface0.303.07902.040.355.302.000.010.253.2E-60.00Clayey silt0.482.584010.670.573.281.540.010.442.4E-70.15Silty clay0.372.010206.900.472.602.300.020.353.3E-80.20Ice-rich soil0.603.1102011.810.572.602.300.020.552.5E-60.00Mudstone0.203.512009.300.522.301.200.020.251.2E-80.00The models were driven by the measured meteorological data in 2013 from the literature [4, 27, 28]. The measured data for air temperature, rainfall, air moisture, solar radiation, and wind regime during 0.5 h were input into the models. Equations (17) and (23) were used to obtain the upper thermal and water flux boundaries. According to the historical geothermal data from Beiluhe, it is reasonable to take the temperature gradient (0.03°C/m) at the 15 m depth as a lower thermal boundary [31]. The lower water boundary can be viewed as zero flux due to the perennial freezing at 15 m depth. Measured soil temperatures and water content were used as initial values. The 1st of January 2013 was selected as the initial time of the simulation, and the initial soil water content and temperature values for each node were interpolated or extrapolated using the measured data. The temperatures and water content of the soil below 2.2 m were determined from the deep hole detection data in Beiluhe. The models have the same time step, relative tolerance, and mesh size as the model used for validation.


## 2.1. Mathematical-Physical Model

### 2.1.1. Hydrothermal Coupling Model

The liquid water flow in unsaturated soil has a positive correlation with hydraulic conductivity. It is governed by the potential energy of soil matrix, temperature, and gravity. The process obeys the Richards function [29], namely,
(1)ql=qlh+qlT=−Klh∇h+y−KlT∇T,where ql is the liquid water fluxes (m/s).qlh and qlT are the liquid water fluxes under the effect of soil matrix and temperature potential (m/s), respectively. Klh and KlT are the liquid water hydraulic conductivities of hydraulic pressure (m/s) and temperature potential (m2/(K·s)), respectively. h is the hydraulic pressure (m), and y is the gravitational potential (m). T is the temperature in soil (K).In contrast to the mechanism of liquid water flow, it is noted that the vaporous water transport in unsaturated soil is affected by other potential energy except the gravity potential. The process obeys Fick’s law [30]:
(2)qv=qvh+qvT=−Kvh∇h−KvT∇T,where qv is the vaporous water fluxes (m/s).qvh and qvT are the vaporous water fluxes of soil matrix and temperature potential (m/s), respectively. Kvh and KvT are the vapor hydraulic conductivities of hydraulic pressure (m/s) and temperature potential (m2/(K·s)).According to the mass conservation law, the fluid flow in and out of the soil per unit of time keeps the balance. The liquid and vaporous water during the soil freeze-thaw process can be described as follows [4, 27]:
(3)∂θ∂t=∇∙Klhh+y+KlT∇T+∇∙Kvh∇h+KvT∇T−ρvρw∂∂tn−θl−θi,(4)θ=θl+ρiρwθi,where θ is the total volumetric content of water and t is time (s). ρw, ρv, and ρi are the densities of liquid, vaporous, and solid phases of water in soil (kg/m3). n is the soil porosity. θl and θi are the volume contents of liquid and vaporous phases of water (m3/m3). The detailed hydraulic properties can be referred to in the research by Zhang et al. [27, 28]. The right third term in Equation (3) is the variation of equivalent vapor.The following energy conservation equation provides the governing equation for energy transport in a variably saturated, stiff porous material [4, 27]:
(5)Cp∂T∂t=∇·λp∇T−ql∇ClT−qv∇CvT−Lwρw∇qv+Lfρi∂θi∂t−Lwρw∂θv∂t,where λp is the heat conduction coefficient of soil (W/(m·K)). Cl, Cv, and Cp are the volumetric thermal capacity of liquid water, vaporous water, and soil (J/(kg·K)), respectively. Lw and Lf are the latent heat of soil water evaporation and ice-water phase transition (J/kg), respectively. θv is the volumetric content of vaporous water (m3/m3). The left term in Equation (5) is the variation of temperature per unit of time. The right terms in Equation (5) are the heat conduction, the convective heat of the liquid and vaporous water, the vapor diffusion heat, and the phase transition heat of ice-liquid and liquid-vapor.The volumetric thermal capacity and heat conduction coefficient of frozen soil are comprehensive parameters affected by thermal parameters and the contents of soil components. It can be expressed using the following formulas [31]:
(6)Cp=Cs1−n+Clθl+Ciθi,(7)λp=λs1−nλlθlλiθiλaθv,where Ci and Cs are the volumetric thermal capacity of the ice and soil matrix, respectively (J/(m3·K)). λs, λl, λi, and λa are the heat conduction coefficients of the soil matrix, liquid, solid, and vaporous water (W/(m·K)).The volumetric liquid water and ice content in the soil during the freeze-thaw process can be calculated by the relationship between soil temperature and equivalent volumetric water content [31], namely,
(8)θumax=a×T−273.15−b,(9)θl=θumax,T<273.15andθl>θumax,θ,otherwise,(10)θi=ρwρiθ−θl,where a and b are two soil characteristic parameters [31]. θumax is the maximum unfrozen water content of frozen soil below the freezing point (m3/m3).The V-G model [32] was applied to establish the relationship between the matrix potential and the liquid water in unsaturated soil:
(11)Se=θl−θrθs−θr=11+aψhωm,h<0,1,h≥0,(12)Klh=KsSel1−1−Se1/mm2∕I,h<0,Ks,h≥0,(13)KlT=KlhhGWT1γ0dγdT,where Se is the soil saturation. θs and θr are the liquid water contents in the saturated and residual states (m3/m3), respectively. aψ is a soil intake reciprocal value. ω and m (=1−1/ω) are empirical parameters. Ks is the hydraulic conductivity in a saturated state (m/s). I is the ice resistance coefficient. GWT is a gain factor, and γ is the soil water surface tension (J/m2), γ0 is the surface tension at 25°C (J/m2).The fundamental equation used to calculate the diffusion of water vapor in soil is Fick’s diffusion law. The diffusion coefficients of water vapor under the water potential gradient of soil and temperature gradient are expressed as follows [33]:
(14)Kvh=DaρwρvsΗrΜgRT,(15)KvT=DaρwηΗrdρvsdT,(16)η=9.5+3θvθs−8.5exp−1+2.6fcθl+θiθs4,where Da is the vapor diffusivity in soil (m2/s). ρvs is the density of saturated water vapor (kg/m3). Ηr is the relative humidity. Μ is the molecular weight of water (kg/mol). g is the gravitational acceleration (m/s2). R is the universal gas constant (J/(mol·K)). η is the enhancement factor of water vapor [34]. fc is the proportion of clay in the soil.


### 2.1.2. Hydrothermal Flux Boundary

The aerodynamic principle at the near-surface is an effective method to characterize the thermal-moisture exchange process between the underlying surface and the atmosphere caused by meteorological conditions such as radiation, rainfall, and evaporation. The ground surface energy balance equation is as follows [35]:
(17)G=Rn−Hs−LwE,where G is the soil surface heat flux (W/m2). Rn, Hs, and LwE are the fluxes of net radiation, sensible heat, and latent heat on the ground surface (W/m2), respectively. The Rn can be written as follows [36]:
(18)Rn=Rns+Rnl=1−αRns↓+σTa40.605+1.77ha−σεsTs41063600,where Rns and Rnl are the surface net shortwave and net longwave radiation (W/m2), respectively. Rns↓ is the solar shortwave radiation (W/m2). α is the albedo on the ground. σ is the Stefan-Boltzmann constant (J/(m2·s·K4)). Ta and Ts are the temperatures of air and ground, respectively (K). ha is the effective volumetric water content of air (kg/m3). εs is the surface emissivity.Soil color has a significant impact on the ground albedo. The pure black soil can absorb almost all visible solar radiation. Its albedo is the smallest. Inversely, the albedo of pure white soil is the largest. The water around the soil particles increases the absorption route of sunlight, which makes the surface reflectivity decrease gradually with the increase of soil moisture. The ground albedo is defined as a function of the surface type and water content [37]:
(19)α=αs+∆αg,(20)∆αg=0.11−0.4θl,where αs is the saturated soil albedo and ∆αg is the increment of albedo induced by dehumidification and drying of the soil. The ground albedo of different soil colors and moisture is in Table 1 [38].Table 1
Surface albedo of different surface types and water contents.
ParameterWavelengthColor (from white (1) to black (8))12345678Saturated soil<0.7μm0.120.110.100.090.080.070.060.05≥0.7μm0.240.220.200.180.160.140.120.10Dry soil<0.7μm0.230.220.200.180.160.140.120.10≥0.7μm0.460.440.400.360.320.280.240.20The surface evaporation water flux (kg/(m2·s)) is as follows:
(21)E=ρvss−ρvara+rs,where ρvss and ρva are the vapor densities in the near-surface soil and the atmosphere, respectively (kg/m3). ra is the aerodynamic resistance to heat transfer (s/m). rs is the surface soil resistance to vapor (s/m).The surface sensible heating flux is as follows:(22)Hs=CaTs−Tara,where Ca is the volumetric thermal capacity of the near-surface air (J/(m3·K)). The G obtained from formulas (17)–(22) was used as the upper thermal boundary.In the Beiluhe area, there is relatively little snowfall. The influence of snowfall on the ground can be overlooked [27, 28]. The soil surface water flux (m/s) can be calculated by the E, ρw, and precipitation rate (Prain, m/s):
(23)Q=Eρw−Prain.


## 2.1.1. Hydrothermal Coupling Model

The liquid water flow in unsaturated soil has a positive correlation with hydraulic conductivity. It is governed by the potential energy of soil matrix, temperature, and gravity. The process obeys the Richards function [29], namely,
(1)ql=qlh+qlT=−Klh∇h+y−KlT∇T,where ql is the liquid water fluxes (m/s).qlh and qlT are the liquid water fluxes under the effect of soil matrix and temperature potential (m/s), respectively. Klh and KlT are the liquid water hydraulic conductivities of hydraulic pressure (m/s) and temperature potential (m2/(K·s)), respectively. h is the hydraulic pressure (m), and y is the gravitational potential (m). T is the temperature in soil (K).In contrast to the mechanism of liquid water flow, it is noted that the vaporous water transport in unsaturated soil is affected by other potential energy except the gravity potential. The process obeys Fick’s law [30]:
(2)qv=qvh+qvT=−Kvh∇h−KvT∇T,where qv is the vaporous water fluxes (m/s).qvh and qvT are the vaporous water fluxes of soil matrix and temperature potential (m/s), respectively. Kvh and KvT are the vapor hydraulic conductivities of hydraulic pressure (m/s) and temperature potential (m2/(K·s)).According to the mass conservation law, the fluid flow in and out of the soil per unit of time keeps the balance. The liquid and vaporous water during the soil freeze-thaw process can be described as follows [4, 27]:
(3)∂θ∂t=∇∙Klhh+y+KlT∇T+∇∙Kvh∇h+KvT∇T−ρvρw∂∂tn−θl−θi,(4)θ=θl+ρiρwθi,where θ is the total volumetric content of water and t is time (s). ρw, ρv, and ρi are the densities of liquid, vaporous, and solid phases of water in soil (kg/m3). n is the soil porosity. θl and θi are the volume contents of liquid and vaporous phases of water (m3/m3). The detailed hydraulic properties can be referred to in the research by Zhang et al. [27, 28]. The right third term in Equation (3) is the variation of equivalent vapor.The following energy conservation equation provides the governing equation for energy transport in a variably saturated, stiff porous material [4, 27]:
(5)Cp∂T∂t=∇·λp∇T−ql∇ClT−qv∇CvT−Lwρw∇qv+Lfρi∂θi∂t−Lwρw∂θv∂t,where λp is the heat conduction coefficient of soil (W/(m·K)). Cl, Cv, and Cp are the volumetric thermal capacity of liquid water, vaporous water, and soil (J/(kg·K)), respectively. Lw and Lf are the latent heat of soil water evaporation and ice-water phase transition (J/kg), respectively. θv is the volumetric content of vaporous water (m3/m3). The left term in Equation (5) is the variation of temperature per unit of time. The right terms in Equation (5) are the heat conduction, the convective heat of the liquid and vaporous water, the vapor diffusion heat, and the phase transition heat of ice-liquid and liquid-vapor.The volumetric thermal capacity and heat conduction coefficient of frozen soil are comprehensive parameters affected by thermal parameters and the contents of soil components. It can be expressed using the following formulas [31]:
(6)Cp=Cs1−n+Clθl+Ciθi,(7)λp=λs1−nλlθlλiθiλaθv,where Ci and Cs are the volumetric thermal capacity of the ice and soil matrix, respectively (J/(m3·K)). λs, λl, λi, and λa are the heat conduction coefficients of the soil matrix, liquid, solid, and vaporous water (W/(m·K)).The volumetric liquid water and ice content in the soil during the freeze-thaw process can be calculated by the relationship between soil temperature and equivalent volumetric water content [31], namely,
(8)θumax=a×T−273.15−b,(9)θl=θumax,T<273.15andθl>θumax,θ,otherwise,(10)θi=ρwρiθ−θl,where a and b are two soil characteristic parameters [31]. θumax is the maximum unfrozen water content of frozen soil below the freezing point (m3/m3).The V-G model [32] was applied to establish the relationship between the matrix potential and the liquid water in unsaturated soil:
(11)Se=θl−θrθs−θr=11+aψhωm,h<0,1,h≥0,(12)Klh=KsSel1−1−Se1/mm2∕I,h<0,Ks,h≥0,(13)KlT=KlhhGWT1γ0dγdT,where Se is the soil saturation. θs and θr are the liquid water contents in the saturated and residual states (m3/m3), respectively. aψ is a soil intake reciprocal value. ω and m (=1−1/ω) are empirical parameters. Ks is the hydraulic conductivity in a saturated state (m/s). I is the ice resistance coefficient. GWT is a gain factor, and γ is the soil water surface tension (J/m2), γ0 is the surface tension at 25°C (J/m2).The fundamental equation used to calculate the diffusion of water vapor in soil is Fick’s diffusion law. The diffusion coefficients of water vapor under the water potential gradient of soil and temperature gradient are expressed as follows [33]:
(14)Kvh=DaρwρvsΗrΜgRT,(15)KvT=DaρwηΗrdρvsdT,(16)η=9.5+3θvθs−8.5exp−1+2.6fcθl+θiθs4,where Da is the vapor diffusivity in soil (m2/s). ρvs is the density of saturated water vapor (kg/m3). Ηr is the relative humidity. Μ is the molecular weight of water (kg/mol). g is the gravitational acceleration (m/s2). R is the universal gas constant (J/(mol·K)). η is the enhancement factor of water vapor [34]. fc is the proportion of clay in the soil.


## 2.1.2. Hydrothermal Flux Boundary

The aerodynamic principle at the near-surface is an effective method to characterize the thermal-moisture exchange process between the underlying surface and the atmosphere caused by meteorological conditions such as radiation, rainfall, and evaporation. The ground surface energy balance equation is as follows [35]:
(17)G=Rn−Hs−LwE,where G is the soil surface heat flux (W/m2). Rn, Hs, and LwE are the fluxes of net radiation, sensible heat, and latent heat on the ground surface (W/m2), respectively. The Rn can be written as follows [36]:
(18)Rn=Rns+Rnl=1−αRns↓+σTa40.605+1.77ha−σεsTs41063600,where Rns and Rnl are the surface net shortwave and net longwave radiation (W/m2), respectively. Rns↓ is the solar shortwave radiation (W/m2). α is the albedo on the ground. σ is the Stefan-Boltzmann constant (J/(m2·s·K4)). Ta and Ts are the temperatures of air and ground, respectively (K). ha is the effective volumetric water content of air (kg/m3). εs is the surface emissivity.Soil color has a significant impact on the ground albedo. The pure black soil can absorb almost all visible solar radiation. Its albedo is the smallest. Inversely, the albedo of pure white soil is the largest. The water around the soil particles increases the absorption route of sunlight, which makes the surface reflectivity decrease gradually with the increase of soil moisture. The ground albedo is defined as a function of the surface type and water content [37]:
(19)α=αs+∆αg,(20)∆αg=0.11−0.4θl,where αs is the saturated soil albedo and ∆αg is the increment of albedo induced by dehumidification and drying of the soil. The ground albedo of different soil colors and moisture is in Table 1 [38].Table 1
Surface albedo of different surface types and water contents.
ParameterWavelengthColor (from white (1) to black (8))12345678Saturated soil<0.7μm0.120.110.100.090.080.070.060.05≥0.7μm0.240.220.200.180.160.140.120.10Dry soil<0.7μm0.230.220.200.180.160.140.120.10≥0.7μm0.460.440.400.360.320.280.240.20The surface evaporation water flux (kg/(m2·s)) is as follows:
(21)E=ρvss−ρvara+rs,where ρvss and ρva are the vapor densities in the near-surface soil and the atmosphere, respectively (kg/m3). ra is the aerodynamic resistance to heat transfer (s/m). rs is the surface soil resistance to vapor (s/m).The surface sensible heating flux is as follows:(22)Hs=CaTs−Tara,where Ca is the volumetric thermal capacity of the near-surface air (J/(m3·K)). The G obtained from formulas (17)–(22) was used as the upper thermal boundary.In the Beiluhe area, there is relatively little snowfall. The influence of snowfall on the ground can be overlooked [27, 28]. The soil surface water flux (m/s) can be calculated by the E, ρw, and precipitation rate (Prain, m/s):
(23)Q=Eρw−Prain.


## 2.2. Model Implementation and Validation

The above mathematical-physical model formed by some high-order nonlinear differential equations was solved by the COMSOL Multiphysics 4.3a finite element software. COMSOL has powerful multiphysics modeling capabilities, making it easy for users to use the various built-in physics modules. Moreover, users can insert partial differential equations (PDEs) and set how to establish the connection between multiphysics equations for the realization of sophisticated physical processes. The governing equations in the model are Equations (3) and (5), with temperature and total volumetric water content as the essential variables, which can be calculated by the PDE tool of COMSOL. The standard form of a generalized PDE is as follows:
(24)ea∂2U∂t2+da∂U∂t+∇∙Γ=f,where ea is a damping coefficient. da is a mass coefficient. Γ is a flux vector, and f is a source term. We can input the governing equations of the above mathematical model with reference to the standard form of PDE, where U is a variable. The θ and T can be input, and the remaining items can be input parameters, partial derivatives, etc. The above equations were solved by invoking the COMSOL transient solver, which solves a possibly nonlinear system of equations at each time step through a set of iterative techniques based on Newton’s method. The iterative technique of Newton’s method allows the function and its Jacobian matrix (derivative) to be computed at each time step. If the solver is poorly converged, it will reduce the requested time step and try to compute the result. The iterative calculation is repeated until a relatively accurate result is obtained within the set maximum number of Newton’s iterations. Therefore, we can improve the model convergence and get an accurate result through the internal setting of the nonlinear solver in advance. The steps of the backward difference method (BDF) in the solver were set to “Strict,” which helps to get more accurate solutions within a finite time step. Discretization in time varies between a minimum time step of 60 s and a maximum time step of 900 s. The allowable relative tolerance is 10-3. Meanwhile, a smaller mesh of spatial discretization, namely, 1 cm, was used.The reliability of the aforesaid numerical model was confirmed by comparing the expected and measured outcomes. MBE (mean bias error) andd (agreement index) were used as the key evaluation indicators. The simulated temperatures and volumetric water contents were compared to the measured values at the test site from October 1, 2012, to October 1, 2013. As shown in Table 2, the processed statistics indicate that the verification results of the temperature at 5, 15, 25, 35, 75, 150, 160, 170, 200, and 220 cm depths are d>0.802 and MBE<1.204°C. Meanwhile, the result is d>0.853 and MBE<1.353% for the water content [4, 28]. A deviation of 2~3°C from the measured soil temperature is regarded as acceptable [35]. The d index of the estimated value steadily approaches 1.0, which is more favorable for enhancing the matching degree between the anticipated and measured values [35]. The simulated results correlate well with the observed findings, and the error control is good. It indicates that the model is applicable to computing the moisture-heat condition in the shallow soil of various underlying surfaces.Table 2
Summary statistics for differences between the measured and simulated values at different depths [39].
Depth (cm)Soil temperatureLiquid water contentsd (-)MBE (°C)d (-)MBE (%)50.9620.6160.9430.267150.9680.7560.9770.051250.9730.3930.8551.353350.976-0.0520.9991.031750.9570.4180.8530.7181500.917-0.5730.8910.8081600.881-1.0360.9160.4131700.826-1.2040.9020.3112000.8020.9580.9780.1752200.862-0.7550.8931.288


## 2.3. Model Application

The model has been used in different application scenarios and engineering environments. The reliability is high [27, 28]. The NG, AP, and GP are some of the typical underlying surfaces in the central QTP. It is of great engineering and practical significance to carry out extensive research. Therefore, we selected several underlying surfaces near the frozen soil observation station of Beiluhe as the research objects, which allows us to use the same meteorological factors as the model-driven data in the numerical models of different underlying surfaces.


### 2.3.1. Research Location

The place of research lies in the permafrost test section in the Beiluhe Basin in the central QTP (34.845°N, 92.940°E, 4620 m a.s.l.) (Figure 1). The NG is an alpine meadow with a high content of silty clay and less vegetation cover [4, 15]. In the test section, there are two kinds of artificially covered underlying surfaces: the AP and GP. The soil hydraulic and thermal parameters on the different underlying surfaces were determined by the field-tested data. The upper limit of permafrost is 2 to 3 m, and the annual mean ground temperature is about -1.8 to -0.4°C [28]. The average air temperature for the time window between 2004 and 2013 was −3.4°C [4]. In 2013, the yearly average air temperature was -3.31°C, with the maximum and lowest monthly temperatures of 8.86°C (August) and -14.14°C (January), respectively. The analysis of meteorological data reveals that the rainfall is characterized by multiple times and small amounts, and the rainfall mainly occurs in May-September of each year. There are fewer snow-covered events on the ground surface in cold seasons, generally no more than a week [27, 28].Figure 1
Schematic location of the study site (from LocaSpaceViewer).


### 2.3.2. Model Setup and Forcings

The three kinds of typical underlying surfaces (i.e., AP, GP, and NG) near the frozen soil observation station of Beiluhe were selected as the research objects. A model with a depth of 15 m was established. The downward coordinate is negative (Figure2). It was reduced to a one-dimensional model since the test location is flat and there is no visible hydrothermal lateral advection. For these three different underlying surface models, the soils within 1 m below the ground surface were set as GP (0.30 m gravel packing and 0.70 m subsand), AP (0.10 m asphalt surface layer, 0.20 m sandy gravel, and 0.70 m clayey silt), and NG (1.0 m clayey silt), respectively. The properties of soils below 1.0 m depth are the same. The thermal-moisture parameters required for the calculation of the coupled model are shown in Table 3. The data in Ref. [31] was obtained from a large number of experimental results on different soils in different research sites, including the QTP.Figure 2
Vertical soil-layer distributions at the research location.Table 3
Soil-layer thermal-moisture parameters of the underlying surface [31].
Soil typenλs (W/(m·K))Cs (J/(m3·K))abaψ (1/m)ωθr (m3/m3)θs (m3/m3)Ks (m/s)fcAsphalt surface0.151.595014.410.56021.562.000.010.072.3E-110.00Gravel surface0.303.07902.040.355.302.000.010.253.2E-60.00Clayey silt0.482.584010.670.573.281.540.010.442.4E-70.15Silty clay0.372.010206.900.472.602.300.020.353.3E-80.20Ice-rich soil0.603.1102011.810.572.602.300.020.552.5E-60.00Mudstone0.203.512009.300.522.301.200.020.251.2E-80.00The models were driven by the measured meteorological data in 2013 from the literature [4, 27, 28]. The measured data for air temperature, rainfall, air moisture, solar radiation, and wind regime during 0.5 h were input into the models. Equations (17) and (23) were used to obtain the upper thermal and water flux boundaries. According to the historical geothermal data from Beiluhe, it is reasonable to take the temperature gradient (0.03°C/m) at the 15 m depth as a lower thermal boundary [31]. The lower water boundary can be viewed as zero flux due to the perennial freezing at 15 m depth. Measured soil temperatures and water content were used as initial values. The 1st of January 2013 was selected as the initial time of the simulation, and the initial soil water content and temperature values for each node were interpolated or extrapolated using the measured data. The temperatures and water content of the soil below 2.2 m were determined from the deep hole detection data in Beiluhe. The models have the same time step, relative tolerance, and mesh size as the model used for validation.


## 2.3.1. Research Location

The place of research lies in the permafrost test section in the Beiluhe Basin in the central QTP (34.845°N, 92.940°E, 4620 m a.s.l.) (Figure 1). The NG is an alpine meadow with a high content of silty clay and less vegetation cover [4, 15]. In the test section, there are two kinds of artificially covered underlying surfaces: the AP and GP. The soil hydraulic and thermal parameters on the different underlying surfaces were determined by the field-tested data. The upper limit of permafrost is 2 to 3 m, and the annual mean ground temperature is about -1.8 to -0.4°C [28]. The average air temperature for the time window between 2004 and 2013 was −3.4°C [4]. In 2013, the yearly average air temperature was -3.31°C, with the maximum and lowest monthly temperatures of 8.86°C (August) and -14.14°C (January), respectively. The analysis of meteorological data reveals that the rainfall is characterized by multiple times and small amounts, and the rainfall mainly occurs in May-September of each year. There are fewer snow-covered events on the ground surface in cold seasons, generally no more than a week [27, 28].Figure 1
Schematic location of the study site (from LocaSpaceViewer).


## 2.3.2. Model Setup and Forcings

The three kinds of typical underlying surfaces (i.e., AP, GP, and NG) near the frozen soil observation station of Beiluhe were selected as the research objects. A model with a depth of 15 m was established. The downward coordinate is negative (Figure2). It was reduced to a one-dimensional model since the test location is flat and there is no visible hydrothermal lateral advection. For these three different underlying surface models, the soils within 1 m below the ground surface were set as GP (0.30 m gravel packing and 0.70 m subsand), AP (0.10 m asphalt surface layer, 0.20 m sandy gravel, and 0.70 m clayey silt), and NG (1.0 m clayey silt), respectively. The properties of soils below 1.0 m depth are the same. The thermal-moisture parameters required for the calculation of the coupled model are shown in Table 3. The data in Ref. [31] was obtained from a large number of experimental results on different soils in different research sites, including the QTP.Figure 2
Vertical soil-layer distributions at the research location.Table 3
Soil-layer thermal-moisture parameters of the underlying surface [31].
Soil typenλs (W/(m·K))Cs (J/(m3·K))abaψ (1/m)ωθr (m3/m3)θs (m3/m3)Ks (m/s)fcAsphalt surface0.151.595014.410.56021.562.000.010.072.3E-110.00Gravel surface0.303.07902.040.355.302.000.010.253.2E-60.00Clayey silt0.482.584010.670.573.281.540.010.442.4E-70.15Silty clay0.372.010206.900.472.602.300.020.353.3E-80.20Ice-rich soil0.603.1102011.810.572.602.300.020.552.5E-60.00Mudstone0.203.512009.300.522.301.200.020.251.2E-80.00The models were driven by the measured meteorological data in 2013 from the literature [4, 27, 28]. The measured data for air temperature, rainfall, air moisture, solar radiation, and wind regime during 0.5 h were input into the models. Equations (17) and (23) were used to obtain the upper thermal and water flux boundaries. According to the historical geothermal data from Beiluhe, it is reasonable to take the temperature gradient (0.03°C/m) at the 15 m depth as a lower thermal boundary [31]. The lower water boundary can be viewed as zero flux due to the perennial freezing at 15 m depth. Measured soil temperatures and water content were used as initial values. The 1st of January 2013 was selected as the initial time of the simulation, and the initial soil water content and temperature values for each node were interpolated or extrapolated using the measured data. The temperatures and water content of the soil below 2.2 m were determined from the deep hole detection data in Beiluhe. The models have the same time step, relative tolerance, and mesh size as the model used for validation.


## 3. Results and Discussion

### 3.1. Surface Energy Budgets

The variation of the monthly average for the four surface energy components is illustrated in Figure3. The following is the relationship between the Rn in May to September: AP>GP>NG. The variation of Rn is obvious during the period, and the Rn of the AP in May is 9.82 W/m2 (7.4%) and 10.07 W/m2 (7.6%) larger than those of the GP and NG. However, the Rn of the NG is larger than that of the GP in other months. The maximum difference is 2.87 W/m2 (Figure 3(a)). It can be seen from Equation (22) that the Hs is positively correlated with the temperature difference between the ground and the air. The Hs has the following relationship from June to September: AP>NG>GP. The difference of Hs is obvious for the three underlying surfaces, showing that in July, the Hs of the AP is 46.89 W/m2 (118.6%) and 33.96 W/m2 (64.7%) higher than those of the GP and NG, respectively. However, the relationship between GP and NG shows the opposite result in other months. The GP has a maximum of 7.47 W/m2 higher Hs than the NG (Figure 3(b)). The porosity and saturated hydraulic conductivity of GP are much larger, and the water infiltration and evaporation are much stronger than those of the NG and AP under the same rainfall condition. Therefore, the LwE decreases significantly from the GP to the AP. The LwE difference in the three underlying surfaces is evident in rainy months. The monthly average LwE of the AP, NG, and GP in July is 29.4, 62.3, and 81.7 W/m2, respectively. The LwE of the AP decreases by 52.8% and 64.1% compared with those of the NG and GP. During the frozen period, the evaporation of the ground surface is tiny, and the differences between LwE in the three underlying surfaces are small (Figure 3(c)). In general, the difference in the surface energy balance caused by the different underlying surfaces is larger in summer and autumn and smaller in others.Figure 3
Monthly average of surface energy components at different underlying surfaces.
(a)
Rn(b)
Hs(c)
LwE(d)
GThe annual average value of the surface energy components is shown in Figure4. The Rn of the GP and the AP increases by 0.1 and 7.3 W/m2 compared with that of the NG, and the increase rates are 0.12% and 8.98%, respectively. The Hs of the AP increases by 18.8 and 19.1 W/m2, i.e., 35.4% and 36.2% higher than those of the NG and GP, respectively. The LwE of the AP is 12.5 W/m2 (48.4%) and 13.0 W/m2 (49.4%) smaller than those of the NG and GP, respectively. The G of the AP increases by 0.9 and 1.2 W/m2, i.e., 36.0% and 54.5% higher than those of the NG and GP, respectively. Therefore, compared with the NG, the AP can increase Rn and Hs and reduce LwE to increase G on the ground surface. The GP mainly reduces the G by increasing LwE. This is consistent with the observation results for the AP from Zhang et al. in the QTP [15, 16]. It is noted that the increase or decrease of Hs of the AP is a relative value compared with that of the NG. In the above theoretical model, the Hs was determined by the aerodynamic theory. However, the Hs was calculated by the friction velocity and the characteristic scale of potential temperature in the literature [15]. Moreover, it is observed that the sensible heat of the AP is greater than that of the NG in individual months of the warm season. It is believed that it is related to the distribution difference in rainfall [40].Figure 4
The annual average value of energy components at different underlying surfaces.


### 3.2. Soil Water Transport

It is noted that a positive or negative water flux value only indicates that water is migrating up or down. Therefore, the absolute value of the parameter is uniformly used for comparative analysis in this paper. Figure5(a) shows the monthly average changes of qlh at 5 cm depth of different underlying surfaces. During the frozen period (i.e., January–February and November–December), the soil transmissibility coefficient is so low that the migration of liquid water is difficult, and the differences between qlh of the different underlying surfaces are not obvious. The qlh at the near-surface gradually migrate down towards the deep soil as the subsoil thaws and rainfall increases. In July, the qlh at 5 cm depth are −2.47×10−11m/s (AP), −4.48×10−9m/s (GP), and −1.63×10−9m/s (NG), respectively. In other words, compared with the NG, the AP can significantly reduce the qlh, while the GP is the opposite. In October, the liquid water migrates upward to the freezing front as the active layers gradually freeze from the near-surface to the permafrost table. Combining the Ks in Table 3, the qlh in October are 3.35×10−11m/s (AP), 2.07×10−9m/s (GP), and 3.12×10−10m/s (NG), respectively. It demonstrates that liquid water can migrate upward at a faster pace under the GP, while the AP is the opposite. As depth increases, the difference in qlh between the different underlying surfaces gradually decreases (Figure 5(b)). Meanwhile, the qvh is 4 orders of magnitude less than that of the qlh. Figures 5(c) and 5(d) show that they can be ignored.Figure 5
Soil water migration in different underlying surfaces.
(a)
Monthly mean ofqlh(b)
Annual mean ofqlh(c)
Monthly mean ofqvh(d)
Annual mean ofqvh(e)
Monthly mean ofqlT(f)
Annual mean ofqlT(g)
Monthly mean ofqvT(h)
Annual mean ofqvTThe monthly average variation of theqlT and qvT is shown in Figures 5(e)–5(h). During the frozen phase, the qlT and qvT are low. During the melting phase, they are high. The variation of qlT is similar to that of qlh but less an order of magnitude than that of qlh. The variation of qvT is more obvious in the different underlying surfaces. In July, it shows that the qvT at 5 cm depth are −2.69×10−10m/s (AP), −1.55×10−10m/s (GP), and −2.05×10−10m/s (NG), respectively. The qlT at 5 cm depth in July are −5.31×10−12m/s (AP), −1.24×10−11m/s (GP), and −1.09×10−10m/s (NG), respectively. Therefore, the predominant mode of transfer under a temperature gradient is revealed to be water vapor transport. The variation of qlT in the three underlying surfaces is not obvious. In a word, compared with the NG, the water transport of the AP is dominated by the qvT, while the GP is dominated by the qlh.


### 3.3. Soil Water Content

Figure6(a) illustrates that the monthly average of θ under the three underlying surfaces is significant during the thawing period (i.e., from April to September). The θ of AP, NG, and GP increases by 0.013, 0.048, and 0.053, respectively, and the monthly average increased amplitudes are 2.1% (AP), 7.6% (NG), and 11.5% (GP), respectively. However, during the frozen period, the θ of AP in December is 158.5% and 12.3% higher than those of the GP and NG, respectively. Moreover, the annual average value at 5 cm depth of the AP is 4.0% and 43.9% higher than the ones of the NG and GP, respectively. For the GP, the water content in soil at 25 cm or much deeper depth gradually increases, while the soil water in the 5 to 75 cm depth below the AP initially declines and then increases (Figure 6(b)). That is, the good connectivity between the soil and the atmosphere in the active layers of the NG and GP promotes the exchange of moisture and heat at the boundary, and the transmissibility coefficient of the GP is larger than that of the NG. However, the AP is a top-closed system. Therefore, rainfall infiltration and evaporation have a significant impact on the water content of the NG and GP but little influence on the AP. Meanwhile, it also indicates the existence of a moisture accumulation phenomenon under the closed pavement. It is related to the vaporous water that migrates upward to the lower layers of the pavement under the effect of the soil temperature field of the AP [16, 41]. It is consistent with the results of the previous section. However, it is not a good thing for cold region engineering. The research shows that the high aquifer (average volumetric water content is 22–26%) of the water vapor phase transition formed under the AP promotes the accumulation of heat in the soil [15].Figure 6
Soil water content under different underlying surfaces.
(a)
Monthly mean ofθ(b)
Annual mean ofθ


### 3.4. Soil Energy Transfer

Figure7 gives the monthly average and annual average variation of the energy transport components under the different underlying surfaces. During the frozen period, the effect of water transport on energy transfer is weak, and most soil energy is transferred through heat conduction. The differences in soil energy transport components of the three underlying surfaces are small. During the thawing period, the difference in water transport caused by different underlying surfaces is obvious, which results in a significant difference in energy transfer. Taking the change in energy components at 5 cm depth of the three underlying surfaces in July as an example (Figures 7(a)–7(h)), the monthly average convective heat flux of liquid water is 0.04 W/m2 (AP), 2.08 W/m2 (NG), and 5.18 W/m2 (GP), respectively. The monthly average value in the GP is 5.14 W/m2 larger than that of the AP and 3.1 W/m2 larger than that of the NG. The convective heat flux of water vapor in the AP in July is 0.08 W/m2 larger than that of the NG and 0.14 W/m2 larger than that of the GP. The vapor diffusion latent heat of the AP is 0.16 and 0.28 W/m2 larger than those of the NG and GP, respectively. For the annual mean of the heat flux at 5 cm depth, the heat conduction, convective heat of vaporous water, and vapor diffusion latent heat under the three underlying surfaces show a relationship of AP>NG>GP, and the annual maximum difference is 0.21 W/m2, 0.036 W/m2, and 0.069 W/m2, respectively. The annual average liquid water convective heat flux shows AP<NG<GP, and the annual maximum difference is 0.28 W/m2. Moreover, the energy transfer differences between the underlying surfaces gradually decrease with the increase of depth.Figure 7
Soil heat flux under different underlying surfaces.
(a)
Monthly mean of heat conduction(b)
Annual mean of heat conduction(c)
Monthly mean of liquid water heat convection(d)
Annual mean of liquid water heat convection(e)
Monthly mean of vaporous water heat convection(f)
Annual mean of vaporous water heat convection(g)
Monthly mean of vapor diffusion latent heat(h)
Annual mean of vapor diffusion latent heat(i)
Annual mean of water flux(j)
Annual mean of heat fluxIt is shown in Figures7(i) and 7(j) that the qlh under the GP increases obviously in the rainfall period, which results in a significant increase in the liquid water heat convection and a decrease in the heat conduction. However, the soil heat conduction has a strong effect under the AP with the increase of qvT. The impact of the water vapor transfer heat flux and water vapor diffusion latent heat flux under the AP is stronger as well, compared with the GP and NG. The research shows that the vapor is condensed into liquid water when the surface temperature of the AP drops to a relatively low level (e.g., rainfall time). Moreover, the migrated quantity of vapor is considerable under the AP, and it can sometimes exceed the amount of liquid water [15, 42]. The process will significantly increase the heat exchange of water vapor. It is also conducive to the heat conduction process because the liquid water produced by phase change also increases the thermal properties of soil (e.g., the heat conduction coefficient).


### 3.5. Soil Temperature

Figure8(a) shows the monthly mean temperature at 5 cm depth soil. Except in July and August, the soil temperatures of the three underlying surfaces in other months have the following relationship: AP>GP>NG. Comparing the annual average soil temperature changes at different depths (i.e., 5, 25, 50, and 75 cm), the soil temperatures of the three underlying surfaces have the following relationship: AP>GP>NG. The soil temperature differences on different underlying surfaces decrease with the increase of depth.Figure 8
The temperature variation in soil under different underlying surfaces.
(a)
Monthly mean of soil temperature(b)
Annual mean of soil temperatureTo analyze the thermal properties caused by the change of the underlying surface and the influence of variations in the soil heat transfer rate on soil temperature, the relationship between the soil heat conduction fluxqT, the thermal conductivity λp, and temperature gradient ∇T is introduced as follows:
(25)qT=−λp∇T.From the AP to the GP, theG gradually decreases, but the impact of λp and rainfall infiltration on the heat convection strengthens. Based on the ground surface energy balance, the G of the AP is 36.0% larger than that of the NG, and the λs is 40.0% lower than that of the NG. It causes the soil temperature gradient under the AP to be similar to that under the NG. However, because the AP absorbs more heat, the surface temperature is much higher than that of the NG, which leads to the soil temperature at the same depth being higher than that of the NG. A high aquifer formed under the AP also promotes the accumulation of heat. The GP has a 12.0% reduction in the G and a 20.0% increase in the λs compared to the NG. It is obvious that the increased magnitude of λs is much greater than the decreased degree of G, which causes the soil temperature gradient of the GP to gradually decrease (Equation (25)). It means that the decrease in the soil temperature gradient under the GP leads to a higher soil temperature at the same depth. Combined with the water transport (Figure 5) and energy transfer (Figure 7), the liquid water convective effect carries heat from the GP to the interior of the soil under the effect of rainwater infiltration. The whole process accelerates the thickening and degradation of the active layer in permafrost.


### 3.6. Permafrost Table and Stability

Figure9 presents the relationship between rainfall and evaporation on the three underlying surfaces. The evaporation on the ground surface accordingly increases from the AP to the GP. The annual evaporation amounts in the 1.0 m thickness of the soils are 168, 309, and 314 mm for the AP, NG, and GP under the same rainfall conditions in 2013 (an annual rainfall amount of 318 mm). It shows that the AP with higher water-blocking properties and lower water conductivity makes the effect of rainfall infiltration and evaporation the minimum. However, a higher water content formed by the vapor migration exacerbates the heat accumulation under the AP. The evaporation rate in the GP is large, and rainfall does not accumulate in the soil. However, the convective heat transfer is strengthened by the rainfall infiltration through the GP with large porosity and high hydraulic conductivity. Additionally, the decrease in moisture in the shallow active layers decreases the generated latent heat during the water-ice phase transition process. Moreover, compared with the AP and the NG, the decrease of G for the GP is smaller than that of the increased magnitude of the λs, which results in a drop in the thermal gradient and an increase in temperature. Comprehensive functions accelerate the increase in the thickness of permafrost active layers. The upper bounds are 268 cm (NG), 302 cm (GP), and 331 cm (AP), respectively. In general, the construction of AP and GP accelerates the degradation process of permafrost.Figure 9
Surface evaporation of different underlying surfaces and rainfall.
(a)
Monthly evaporation and rainfall(b)
Annual evaporationIn a word, the AP increases soil temperature by the thermal absorption of the pavement and the heat accumulation created by the moisture accumulation. However, the apparent liquid water convective heat exchange under rainfall and evaporation and the high thermal conductivity decrease the temperature gradient of the active layer and hence increase perturbations to the permafrost under the GP. It is similar to the thermal-moisture regime of the shallow sandy soil under the influence of rainfall [27]. It is noted that it is not contradictory between the increased thermal convection caused by the rainfall infiltration and the cooling effect of precipitation on the permafrost. It is related to the comprehensive response of soil temperature to various effects such as energy transfer, evaporation, rainfall, and phase transition. The GP promotes rainfall infiltration to bring more heat into the deeper soil but dissipates the internal heat through the water vapor phase transition and convective transport, which depend on significant evaporation after rainfall. Moreover, it was found that the natural gravel underlying surface is helpful to the stability of frozen soil [6]. However, the hydrothermal characteristics of the engineering gravel pavement are not favorable to permafrost stability. The research on the sand-filled stone layers found that rainfall boosts the liquid water convective heat transfer and decreases the soil heat transfer flux, which is conducive to reducing soil temperature. But the properties of the sand-filled stone layer are major factors in the thermal disturbance of frozen soil [25]. The sandy soil with strong permeability and low water retention capacity increases the convective heat transfer caused by rainfall, which is unfavorable for the thermal stability of permafrost [43]. The difference is related to the thermophysical properties of the soil. A similar problem also exists in the study of AP [14–16]. Therefore, it is necessary to analyze the hydrothermal state from the perspective of the interaction between rainfall and pavement engineering characteristics.The interactive effect between climate change and engineering on the stability of engineering and permafrost is evident, such as pavement cracking, uplift, deformation, collapse, geological disasters, and even ecological deterioration. In this paper, the synthetic effect between the rainfall infiltration process and the change in thermal conductivity on the variation of heat transfer lacks consideration. Moreover, the model was only used to analyze the thermal-moisture dynamics of individual study fields in the QTP, and the calculation time was only in 2013. Therefore, under the conditions of increasing extreme rainfall events, the long-term effect analysis of the soil thermal-moisture and the permafrost stability of different underlying surfaces caused by rainfall infiltration will be a research direction. The permafrost water vapor-thermal coupling analysis method based on the ground surface energy balance can provide a reference for the prediction of impacts of the change of the underlying surface on permafrost stability under a wetting background.


## 3.1. Surface Energy Budgets

The variation of the monthly average for the four surface energy components is illustrated in Figure3. The following is the relationship between the Rn in May to September: AP>GP>NG. The variation of Rn is obvious during the period, and the Rn of the AP in May is 9.82 W/m2 (7.4%) and 10.07 W/m2 (7.6%) larger than those of the GP and NG. However, the Rn of the NG is larger than that of the GP in other months. The maximum difference is 2.87 W/m2 (Figure 3(a)). It can be seen from Equation (22) that the Hs is positively correlated with the temperature difference between the ground and the air. The Hs has the following relationship from June to September: AP>NG>GP. The difference of Hs is obvious for the three underlying surfaces, showing that in July, the Hs of the AP is 46.89 W/m2 (118.6%) and 33.96 W/m2 (64.7%) higher than those of the GP and NG, respectively. However, the relationship between GP and NG shows the opposite result in other months. The GP has a maximum of 7.47 W/m2 higher Hs than the NG (Figure 3(b)). The porosity and saturated hydraulic conductivity of GP are much larger, and the water infiltration and evaporation are much stronger than those of the NG and AP under the same rainfall condition. Therefore, the LwE decreases significantly from the GP to the AP. The LwE difference in the three underlying surfaces is evident in rainy months. The monthly average LwE of the AP, NG, and GP in July is 29.4, 62.3, and 81.7 W/m2, respectively. The LwE of the AP decreases by 52.8% and 64.1% compared with those of the NG and GP. During the frozen period, the evaporation of the ground surface is tiny, and the differences between LwE in the three underlying surfaces are small (Figure 3(c)). In general, the difference in the surface energy balance caused by the different underlying surfaces is larger in summer and autumn and smaller in others.Figure 3
Monthly average of surface energy components at different underlying surfaces.
(a)
Rn(b)
Hs(c)
LwE(d)
GThe annual average value of the surface energy components is shown in Figure4. The Rn of the GP and the AP increases by 0.1 and 7.3 W/m2 compared with that of the NG, and the increase rates are 0.12% and 8.98%, respectively. The Hs of the AP increases by 18.8 and 19.1 W/m2, i.e., 35.4% and 36.2% higher than those of the NG and GP, respectively. The LwE of the AP is 12.5 W/m2 (48.4%) and 13.0 W/m2 (49.4%) smaller than those of the NG and GP, respectively. The G of the AP increases by 0.9 and 1.2 W/m2, i.e., 36.0% and 54.5% higher than those of the NG and GP, respectively. Therefore, compared with the NG, the AP can increase Rn and Hs and reduce LwE to increase G on the ground surface. The GP mainly reduces the G by increasing LwE. This is consistent with the observation results for the AP from Zhang et al. in the QTP [15, 16]. It is noted that the increase or decrease of Hs of the AP is a relative value compared with that of the NG. In the above theoretical model, the Hs was determined by the aerodynamic theory. However, the Hs was calculated by the friction velocity and the characteristic scale of potential temperature in the literature [15]. Moreover, it is observed that the sensible heat of the AP is greater than that of the NG in individual months of the warm season. It is believed that it is related to the distribution difference in rainfall [40].Figure 4
The annual average value of energy components at different underlying surfaces.


## 3.2. Soil Water Transport

It is noted that a positive or negative water flux value only indicates that water is migrating up or down. Therefore, the absolute value of the parameter is uniformly used for comparative analysis in this paper. Figure5(a) shows the monthly average changes of qlh at 5 cm depth of different underlying surfaces. During the frozen period (i.e., January–February and November–December), the soil transmissibility coefficient is so low that the migration of liquid water is difficult, and the differences between qlh of the different underlying surfaces are not obvious. The qlh at the near-surface gradually migrate down towards the deep soil as the subsoil thaws and rainfall increases. In July, the qlh at 5 cm depth are −2.47×10−11m/s (AP), −4.48×10−9m/s (GP), and −1.63×10−9m/s (NG), respectively. In other words, compared with the NG, the AP can significantly reduce the qlh, while the GP is the opposite. In October, the liquid water migrates upward to the freezing front as the active layers gradually freeze from the near-surface to the permafrost table. Combining the Ks in Table 3, the qlh in October are 3.35×10−11m/s (AP), 2.07×10−9m/s (GP), and 3.12×10−10m/s (NG), respectively. It demonstrates that liquid water can migrate upward at a faster pace under the GP, while the AP is the opposite. As depth increases, the difference in qlh between the different underlying surfaces gradually decreases (Figure 5(b)). Meanwhile, the qvh is 4 orders of magnitude less than that of the qlh. Figures 5(c) and 5(d) show that they can be ignored.Figure 5
Soil water migration in different underlying surfaces.
(a)
Monthly mean ofqlh(b)
Annual mean ofqlh(c)
Monthly mean ofqvh(d)
Annual mean ofqvh(e)
Monthly mean ofqlT(f)
Annual mean ofqlT(g)
Monthly mean ofqvT(h)
Annual mean ofqvTThe monthly average variation of theqlT and qvT is shown in Figures 5(e)–5(h). During the frozen phase, the qlT and qvT are low. During the melting phase, they are high. The variation of qlT is similar to that of qlh but less an order of magnitude than that of qlh. The variation of qvT is more obvious in the different underlying surfaces. In July, it shows that the qvT at 5 cm depth are −2.69×10−10m/s (AP), −1.55×10−10m/s (GP), and −2.05×10−10m/s (NG), respectively. The qlT at 5 cm depth in July are −5.31×10−12m/s (AP), −1.24×10−11m/s (GP), and −1.09×10−10m/s (NG), respectively. Therefore, the predominant mode of transfer under a temperature gradient is revealed to be water vapor transport. The variation of qlT in the three underlying surfaces is not obvious. In a word, compared with the NG, the water transport of the AP is dominated by the qvT, while the GP is dominated by the qlh.


## 3.3. Soil Water Content

Figure6(a) illustrates that the monthly average of θ under the three underlying surfaces is significant during the thawing period (i.e., from April to September). The θ of AP, NG, and GP increases by 0.013, 0.048, and 0.053, respectively, and the monthly average increased amplitudes are 2.1% (AP), 7.6% (NG), and 11.5% (GP), respectively. However, during the frozen period, the θ of AP in December is 158.5% and 12.3% higher than those of the GP and NG, respectively. Moreover, the annual average value at 5 cm depth of the AP is 4.0% and 43.9% higher than the ones of the NG and GP, respectively. For the GP, the water content in soil at 25 cm or much deeper depth gradually increases, while the soil water in the 5 to 75 cm depth below the AP initially declines and then increases (Figure 6(b)). That is, the good connectivity between the soil and the atmosphere in the active layers of the NG and GP promotes the exchange of moisture and heat at the boundary, and the transmissibility coefficient of the GP is larger than that of the NG. However, the AP is a top-closed system. Therefore, rainfall infiltration and evaporation have a significant impact on the water content of the NG and GP but little influence on the AP. Meanwhile, it also indicates the existence of a moisture accumulation phenomenon under the closed pavement. It is related to the vaporous water that migrates upward to the lower layers of the pavement under the effect of the soil temperature field of the AP [16, 41]. It is consistent with the results of the previous section. However, it is not a good thing for cold region engineering. The research shows that the high aquifer (average volumetric water content is 22–26%) of the water vapor phase transition formed under the AP promotes the accumulation of heat in the soil [15].Figure 6
Soil water content under different underlying surfaces.
(a)
Monthly mean ofθ(b)
Annual mean ofθ


## 3.4. Soil Energy Transfer

Figure7 gives the monthly average and annual average variation of the energy transport components under the different underlying surfaces. During the frozen period, the effect of water transport on energy transfer is weak, and most soil energy is transferred through heat conduction. The differences in soil energy transport components of the three underlying surfaces are small. During the thawing period, the difference in water transport caused by different underlying surfaces is obvious, which results in a significant difference in energy transfer. Taking the change in energy components at 5 cm depth of the three underlying surfaces in July as an example (Figures 7(a)–7(h)), the monthly average convective heat flux of liquid water is 0.04 W/m2 (AP), 2.08 W/m2 (NG), and 5.18 W/m2 (GP), respectively. The monthly average value in the GP is 5.14 W/m2 larger than that of the AP and 3.1 W/m2 larger than that of the NG. The convective heat flux of water vapor in the AP in July is 0.08 W/m2 larger than that of the NG and 0.14 W/m2 larger than that of the GP. The vapor diffusion latent heat of the AP is 0.16 and 0.28 W/m2 larger than those of the NG and GP, respectively. For the annual mean of the heat flux at 5 cm depth, the heat conduction, convective heat of vaporous water, and vapor diffusion latent heat under the three underlying surfaces show a relationship of AP>NG>GP, and the annual maximum difference is 0.21 W/m2, 0.036 W/m2, and 0.069 W/m2, respectively. The annual average liquid water convective heat flux shows AP<NG<GP, and the annual maximum difference is 0.28 W/m2. Moreover, the energy transfer differences between the underlying surfaces gradually decrease with the increase of depth.Figure 7
Soil heat flux under different underlying surfaces.
(a)
Monthly mean of heat conduction(b)
Annual mean of heat conduction(c)
Monthly mean of liquid water heat convection(d)
Annual mean of liquid water heat convection(e)
Monthly mean of vaporous water heat convection(f)
Annual mean of vaporous water heat convection(g)
Monthly mean of vapor diffusion latent heat(h)
Annual mean of vapor diffusion latent heat(i)
Annual mean of water flux(j)
Annual mean of heat fluxIt is shown in Figures7(i) and 7(j) that the qlh under the GP increases obviously in the rainfall period, which results in a significant increase in the liquid water heat convection and a decrease in the heat conduction. However, the soil heat conduction has a strong effect under the AP with the increase of qvT. The impact of the water vapor transfer heat flux and water vapor diffusion latent heat flux under the AP is stronger as well, compared with the GP and NG. The research shows that the vapor is condensed into liquid water when the surface temperature of the AP drops to a relatively low level (e.g., rainfall time). Moreover, the migrated quantity of vapor is considerable under the AP, and it can sometimes exceed the amount of liquid water [15, 42]. The process will significantly increase the heat exchange of water vapor. It is also conducive to the heat conduction process because the liquid water produced by phase change also increases the thermal properties of soil (e.g., the heat conduction coefficient).


## 3.5. Soil Temperature

Figure8(a) shows the monthly mean temperature at 5 cm depth soil. Except in July and August, the soil temperatures of the three underlying surfaces in other months have the following relationship: AP>GP>NG. Comparing the annual average soil temperature changes at different depths (i.e., 5, 25, 50, and 75 cm), the soil temperatures of the three underlying surfaces have the following relationship: AP>GP>NG. The soil temperature differences on different underlying surfaces decrease with the increase of depth.Figure 8
The temperature variation in soil under different underlying surfaces.
(a)
Monthly mean of soil temperature(b)
Annual mean of soil temperatureTo analyze the thermal properties caused by the change of the underlying surface and the influence of variations in the soil heat transfer rate on soil temperature, the relationship between the soil heat conduction fluxqT, the thermal conductivity λp, and temperature gradient ∇T is introduced as follows:
(25)qT=−λp∇T.From the AP to the GP, theG gradually decreases, but the impact of λp and rainfall infiltration on the heat convection strengthens. Based on the ground surface energy balance, the G of the AP is 36.0% larger than that of the NG, and the λs is 40.0% lower than that of the NG. It causes the soil temperature gradient under the AP to be similar to that under the NG. However, because the AP absorbs more heat, the surface temperature is much higher than that of the NG, which leads to the soil temperature at the same depth being higher than that of the NG. A high aquifer formed under the AP also promotes the accumulation of heat. The GP has a 12.0% reduction in the G and a 20.0% increase in the λs compared to the NG. It is obvious that the increased magnitude of λs is much greater than the decreased degree of G, which causes the soil temperature gradient of the GP to gradually decrease (Equation (25)). It means that the decrease in the soil temperature gradient under the GP leads to a higher soil temperature at the same depth. Combined with the water transport (Figure 5) and energy transfer (Figure 7), the liquid water convective effect carries heat from the GP to the interior of the soil under the effect of rainwater infiltration. The whole process accelerates the thickening and degradation of the active layer in permafrost.


## 3.6. Permafrost Table and Stability

Figure9 presents the relationship between rainfall and evaporation on the three underlying surfaces. The evaporation on the ground surface accordingly increases from the AP to the GP. The annual evaporation amounts in the 1.0 m thickness of the soils are 168, 309, and 314 mm for the AP, NG, and GP under the same rainfall conditions in 2013 (an annual rainfall amount of 318 mm). It shows that the AP with higher water-blocking properties and lower water conductivity makes the effect of rainfall infiltration and evaporation the minimum. However, a higher water content formed by the vapor migration exacerbates the heat accumulation under the AP. The evaporation rate in the GP is large, and rainfall does not accumulate in the soil. However, the convective heat transfer is strengthened by the rainfall infiltration through the GP with large porosity and high hydraulic conductivity. Additionally, the decrease in moisture in the shallow active layers decreases the generated latent heat during the water-ice phase transition process. Moreover, compared with the AP and the NG, the decrease of G for the GP is smaller than that of the increased magnitude of the λs, which results in a drop in the thermal gradient and an increase in temperature. Comprehensive functions accelerate the increase in the thickness of permafrost active layers. The upper bounds are 268 cm (NG), 302 cm (GP), and 331 cm (AP), respectively. In general, the construction of AP and GP accelerates the degradation process of permafrost.Figure 9
Surface evaporation of different underlying surfaces and rainfall.
(a)
Monthly evaporation and rainfall(b)
Annual evaporationIn a word, the AP increases soil temperature by the thermal absorption of the pavement and the heat accumulation created by the moisture accumulation. However, the apparent liquid water convective heat exchange under rainfall and evaporation and the high thermal conductivity decrease the temperature gradient of the active layer and hence increase perturbations to the permafrost under the GP. It is similar to the thermal-moisture regime of the shallow sandy soil under the influence of rainfall [27]. It is noted that it is not contradictory between the increased thermal convection caused by the rainfall infiltration and the cooling effect of precipitation on the permafrost. It is related to the comprehensive response of soil temperature to various effects such as energy transfer, evaporation, rainfall, and phase transition. The GP promotes rainfall infiltration to bring more heat into the deeper soil but dissipates the internal heat through the water vapor phase transition and convective transport, which depend on significant evaporation after rainfall. Moreover, it was found that the natural gravel underlying surface is helpful to the stability of frozen soil [6]. However, the hydrothermal characteristics of the engineering gravel pavement are not favorable to permafrost stability. The research on the sand-filled stone layers found that rainfall boosts the liquid water convective heat transfer and decreases the soil heat transfer flux, which is conducive to reducing soil temperature. But the properties of the sand-filled stone layer are major factors in the thermal disturbance of frozen soil [25]. The sandy soil with strong permeability and low water retention capacity increases the convective heat transfer caused by rainfall, which is unfavorable for the thermal stability of permafrost [43]. The difference is related to the thermophysical properties of the soil. A similar problem also exists in the study of AP [14–16]. Therefore, it is necessary to analyze the hydrothermal state from the perspective of the interaction between rainfall and pavement engineering characteristics.The interactive effect between climate change and engineering on the stability of engineering and permafrost is evident, such as pavement cracking, uplift, deformation, collapse, geological disasters, and even ecological deterioration. In this paper, the synthetic effect between the rainfall infiltration process and the change in thermal conductivity on the variation of heat transfer lacks consideration. Moreover, the model was only used to analyze the thermal-moisture dynamics of individual study fields in the QTP, and the calculation time was only in 2013. Therefore, under the conditions of increasing extreme rainfall events, the long-term effect analysis of the soil thermal-moisture and the permafrost stability of different underlying surfaces caused by rainfall infiltration will be a research direction. The permafrost water vapor-thermal coupling analysis method based on the ground surface energy balance can provide a reference for the prediction of impacts of the change of the underlying surface on permafrost stability under a wetting background.


## 4. Conclusion

Referring to the meteorological data from Beiluhe in the QTP, the thermal-moisture dynamics in permafrost active layers under different underlying surfaces (asphalt pavement, gravel pavement, and natural ground) under the rainfall conditions were simulated. There are several major conclusions as follows:(1)
Compared with the natural ground, the asphalt pavement increases the surface net radiation and sensible heat flux and reduces the surface evaporation latent heat to increase the soil surface heat flux. The surface heat flux of the gravel pavement decreases mainly owing to the increase in the evaporation latent heat(2)
Asphalt pavements with a closed surface and low hydraulic conduction cut off the moisture exchange between the atmosphere and the active layer. The moisture migration and accumulation under the asphalt pavement are dominated by the temperature gradient vapor flux. The moisture of the gravel pavement is greatly affected by rainfall and evaporation, and the liquid water under the water potential gradient is primarily what causes the moisture migration beneath gravel pavements with high hydraulic conductivity(3)
The temperature change in the shallow active layer below the asphalt pavement is reflected in the four soil energy components, in which the soil heat conduction is most evident. The effects of water vapor transfer heat flux and water vapor diffusion latent heat flux are also important. However, the impact of liquid water migration from rainfall infiltration on thermal transport is significant in the shallow active layer of the gravel pavement(4)
Compared with the natural ground, the greater surface heat fluxes of the asphalt pavement and a higher aquifer under the pavement promote heat accumulation, which causes the active layer to become thicker. The liquid water convection of the gravel pavement is evident, and the decrease in heat fluxes is less than the increase in thermal conductivity. In general, the construction of asphalt pavements and gravel pavements accelerates the degradation of permafrost



---

*Source: 2230372-2022-07-28.xml*