# Long-Term Thermal Regimes of Subgrade under a Drainage Channel in High-Altitudinal Permafrost Environment

Authors: Hong Yu; Hongwu Han; Wei Ma; Zekun Ding; Ling Chen
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

In permafrost regions, construction of a channel involves a large amount of excavation activities and changes to surface water body, which can exert great impacts on the thermal regimes of permafrost underlying. In this paper, a coupled mathematical model of heat and moisture transfer was constructed for freeze-thaw soils to investigate the long-term thermal regimes of subgrade beneath a drainage channel built on the Qinghai-Tibet Plateau. Based on the numerical simulations, the thermal regimes of the subgrade both in warm and cold seasons were analyzed within a period of 50 years, as well as the impact of the widths of the channel. The results showed that the channel excavation and flowing water within could lead to a significant underlying permafrost degradation. During the first 30 years, the permafrost beneath the channel mainly experienced a rapid downward degradation. After that, the lateral thermal erosion of the flowing water led to a rapid permafrost degradation beneath the slope of the channel. In cold seasons, the shallow ground beneath the channel would not refreeze due to the flowing water and the thaw bulb actually expanded throughout the year. For the channel with a bottom width of 15 m, the thaw bulb beneath the channel could expand laterally to the natural ground about 10 m far away from the slope shoulder of channel till the 50th year. With different widths, the long-term thermal regimes of the subgrade beneath the channels differed considerably and the maximum difference was at the slope toe of the embankment. With the numerical simulated results, it is recommended that a channel built on permafrost should be wide-and-shallow rather than narrow-and-deep if conditions permit.


---

## Body

## 1. Introduction

Permafrost is defined as ground (soil or rock including ice) with a temperature at or below 0°C for at least two consecutive years [1, 2]. A quarter of the Northern Hemisphere and 17% of the exposed land surface of earth are underlain by permafrost [3]. The distribution of permafrost on earth can be of altitudinal and latitudinal character [4]. In China, the permafrost area accounts for 22% of the total land area, which is mainly distributed on the Qinghai-Tibet Plateau (QTP), the Northeastern China, and the high mountain areas in the Northwestern China. With an average elevation of more than 4000 m, the QTP has the largest distribution of altitudinal permafrost on earth [5]. Compared with latitudinal permafrost in the Northern America and Siberia, the permafrost on the QTP is characterized with high temperatures and thermal unstable, thus being more sensitive to climatic changes and human activities [6–8].The QTP is also known as the Asian Water Tower. It is one of the three major distribution areas of lakes in China, and its lake area accounts for more than 50% of the total in China [9, 10]. In the context of global warming, the climate warming and wetting on the QTP was significant during the past decades, which led to significant changes in lake number and area [11–14]. The recent estimations showed that, since the mid-1990s, the area, level, and volume of lakes on the QTP showed a continuous increase [15, 16]. In a study by Liu et al. [17], a rapid expansion since 2000 was observed for lakes in the endorheic basin on the QTP, particularly in the Hoh Xil region with a continuous permafrost distribution. The rapid expansion of lakes could cause flood disasters, which may exert significant threatens to the production and living of local communities, as well as major engineering infrastructures nearby [18–20].In 2011, the Zonag Lake in the central of the Hoh Xil regions burst suddenly, which was an extreme case of lake evolution on the QTP. After the outburst, a deep and wide gully developed at the breaking point and the flood entered the downstream lakes including the Kusai Lake, Haiding Nor Lake, and Yanhu Lake. At present, the four lakes are already connected together, and the first three ones lost their storage functions, making the Yanhu Lake become a tail-end-lake [21–24]. From 2011 to 2019, the surface area of the Zonag Lake decreased from 269 to 150 km2, while the surface area of Yanhu Lake increased by almost 3 times from 73 to 209 km2 [5, 25]. Considering the rapid changes, some researchers conducted predictions and concluded that an overflow would occur to the Yanhu Lake in a minimum of one to two years [21]. At about 10 km downstream the Yanhu Lake, there are the Qinghai-Tibet Railway, Qinghai-Tibet Highway, Qinghai-Tibet Power Transmission Line, Golmud-Lhasa Oil Pipeline, as well as several communication cables. It is well-known that all of them are of great importance for social and economic ties between the Tibet Autonomous Region and the inner China [26–28]. If the Yanhu Lake burst in an uncontrollable condition, these major linear infrastructures would be in danger and even be destroyed by the outburst flood.To prevent this projected flood disaster, a large-scale drainage project was conducted in 2019 between the potential breaking point of the Yanhu Lake and the Qingshui River at the downstream [9, 29]. When serving as subgrade of man-made infrastructure, the physical and mechanical properties of permafrost are closely related to its thermal regimes [30–33]. Meanwhile, along with seasonal freezing and thawing in the active layer, significant (differential) frost heave and thaw subsidence would occur at the ground surface, which can cause series of damages to structures built upon [34]. Thus, the thermal regime of permafrost is a key variable that determines the bearing capacity and deformation of foundations and the long-term stability and integrity of the structures upon [35–38]. With regard to evolution of permafrost thermal regimes beneath hydraulic structures, the related research studies were scarce as there has been almost no large-scale hydraulic engineering conducted in continuous permafrost regions. In related fields, research studies on thermal interactions between thermokarst lake and underlying permafrost were conducted both in latitudinal and altitudinal permafrost environments. In 1978, a large tundra lake was drained to study the aggradation of permafrost into newly exposed lake-bottom sediments, and the active-layer thickness, snow depth, minimum soil temperatures, near-surface ground ice, soil heave, and permafrost temperatures had been measured for over 20 years [39]. Lin et al. [40] conducted filed observations on water and ground temperatures beneath and around a 2 m deep thermokarst lake on the QTP and found that the mean annual ground temperatures beneath the thermokarst lake were more 5°C higher than those in the surrounding terrain at the same depths. In a study by Niu et al. [41], characteristics of about 250 thermokarst lakes on the QTP were investigated and their influences on permafrost were evaluated based on field measurements in ground temperatures around a lake in the Beiluhe Basin. Through using numerical methods, Wen et al. [42] and Ling and Zhang [43] simulated the temporal and spatial variation of ground temperature and talik thickness beneath an expanding thermokarst lake on the QTP and the Alakan Arctic, respectively. In a study by Wen et al. [42], the lake was assumed to expand linearly with time both in radius and depth, while in a study by Ling and Zhang [43], the lake was assumed to expand linearly with time only in radius. Li et al. [44] conducted a moisture-heat coupled numerical simulation to investigate the permafrost degradation beneath a thermokarst pond with consideration of climate warming. Compared with thermokarst lakes, however, the thermal impact of a drainage channel to underlying permafrost may be more significant and complex due to the excavation of the channel and the running water within the channel. With regard to channels built in cold regions, Li et al. [45] established a moisture-heat-mechanic couple mathematic model and investigated the freeze-thaw influence on the channel using numerical simulations. Considering the effects of the seepage water, Zhang et al. [46, 47] investigated the thermal regimes of the permafrost beneath the dike with various antiseepage measures.In this paper, the drainage channel of the Yanhu Lake was taken as an example and the thermal interactions among the drainage channel, the drainage water body, and the underlying permafrost subgrade were investigated through field observations and numerical simulations. With the Harlan model and Richard equation, a coupled mathematic model for heat and moisture transfer in freezing-thawing soils was established. Then, a numerical model was constructed based on the engineering geotechnical investigation conducted in the field and the dimensions of the drainage channel. The numerical model was validated with field observations carried out at the drainage channel. Then, the long-term evolution of thermal regimes of permafrost beneath the drainage with three bottom widths was investigated after a validation through using numerical simulations. It is hoped that the results of this study can provide scientific references for predicting and warning of the safe operation of the drainage channel.


## 2. Mathematical Model and Governing Equations

### 2.1. Liquid Water Flows

Without considering the effects of water vapor, the equivalent volume of water contentθ in freeze-thaw soil based on the law of mass conservation can be expressed as follows:(1)∂θ∂t=∂θu∂t+ρiρw∂θi∂t=∇qlh,where θ is the equivalent volume of water content (m3·m−3); θu is the volume of unfrozen water content (m3·m−3); θi is the volume of ice content (m3·m−3); ρi is the density of ice (kg·m−3); ρw is the density of water (kg·m−3); t is the time; and qlh is the liquid water flux density which can represent liquid flows due to a pressure head gradient (m·s−1).The migration of liquid water under potential gradient in frozen soil is similar to that in unfrozen unsaturated soil based on the Harlan model and can be described by Richard equation [48, 49]. Then, it is mainly influenced by the factors including water potential gradient and hydraulic conductivity of soils. The values of water potential gradient and hydraulic conductivity in unfrozen and frozen soils differ considerably, while the regulation of water migration can still be assumed to Darcy’s law [50]. Only considering the potential gradient, the flux density of liquid water can be expressed as follows [51, 52]:(2)qlh=−Klh∇h+y,where y is the vertical coordinate (m); h is the pressure head (m); and Klh is the water conductivity coefficient of liquid water under the action of soil-water potential gradient (m·s−1). Then, considering the equivalent volume water content θ as a dependent variable, the mass conservation equation of liquid water in unfrozen and frozen soils can be written as follows [53]:(3)∂θ∂t=∇·Klh∇h+y.The relationship between the amount of water and energy in the soil can be reflected by soil-water characteristic curve (SWCC). Then, SWCC also can represent the relationship between matric, volume water content, and saturation in frozen and unfrozen soils [52]. In this study, the van Genuchten model and Mualem model are adopted to describe the hydraulic properties of unsaturated freeze-thaw soil [53, 54]. The hydraulic properties can be written as follows:(4)h=−Se1/−m−11/nα,(5)Se=θl−θrθs−θr=11+αhnm,h<0,1,h≥0,(6)K=KsSel1−1−Se1/mm2,h<0,Ks,h≥0,where Se is the effective saturation; Ks is the saturation water conductivity coefficient (m·s−1); θl, θs, and θr are liquid water content, saturated liquid water content, and residual water content (m3·m−3); α is the derivative of the soil intake value (m−1); m = 1–1/n; and n and l are experience parameters, and Mualem suggested that l could be determined as 0.5 [54].


### 2.2. Heat Transfer

Compared with the heat conduction with phase change of ice to water, the energy released by the heat convection was very small and could be neglected during the heat transfer analysis of freeze-thaw soil [55–57]. Then, the governing equations of heat transfer in freeze-thaw soil can be written as follows [49]:(7)Cm∂T∂t=∂∂xλm∂T∂x+∂∂yλm∂T∂y,where Cm is the equivalent volume heat capacity and λm is the equivalent thermal conductivity. According to the method of sensible heat capacity, equivalent volume heat capacity Cm and equivalent thermal conductivity λm in freeze-thaw soils can be written as follows [55, 58]:(8)Cm=Cf,T<Tf−ΔT,Cu+Cf2+Ls2ΔT,Tf−ΔT≤T≤Tf+ΔT,Cu,T>Tf+ΔT,(9)λm=λf,T<Tf−ΔT,λu−λf2T−Tm−ΔT+λf,Tf−ΔT≤T≤Tf+ΔT,λu,T>Tf+ΔT,where Tf ± ∆T is the temperature range of phase change; Cu and λu are the volume heat capacity and thermal conductivity of unfrozen soil; Cf and λf are the volume heat capacity and thermal conductivity of frozen soil; and Ls is the latent heat of phase change in per unit volume. While the temperature ranges of phase change in different soil are difficult to gain, the heat transfer equation adopted in this study is written as follows [49]:(10)C∂T∂t=∂∂xλ∂T∂x+∂∂yλ∂T∂y+L∂θi∂tρi,where C is the volume heat capacity; λ is the thermal conductivity; and L is the latent heat of freezing of liquid water (approximately 3.34 × 105 J·kg−1).After soil freezing, not all liquid water can be transferred to solid ice. The rest of liquid water is called as unfrozen water. Then, the unfrozen water content is related to factors including the temperature, pressure, water salinity, mineralogy, soil specific surface area, and soil surface chemistry [59–62]. Based on previous work and related theories, the empirical expression as follows is used to determine the maximum unfrozen water content in the freezing process [49, 50, 52]:(11)θumax=aT−273.15−b,where a and b are parameters related to the soil properties. The volume liquid water content can be determined by temperature and equivalent volume [50, 52]:(12)θl=θ,others,θumax,T<Tfandθ>θumax,where Tf is the freezing point of saturated soil (0°C). Previous studies showed that the freezing point of soil is not a fixed value, and the ice grows only when the water content exceeds θumax [63].


## 2.1. Liquid Water Flows

Without considering the effects of water vapor, the equivalent volume of water contentθ in freeze-thaw soil based on the law of mass conservation can be expressed as follows:(1)∂θ∂t=∂θu∂t+ρiρw∂θi∂t=∇qlh,where θ is the equivalent volume of water content (m3·m−3); θu is the volume of unfrozen water content (m3·m−3); θi is the volume of ice content (m3·m−3); ρi is the density of ice (kg·m−3); ρw is the density of water (kg·m−3); t is the time; and qlh is the liquid water flux density which can represent liquid flows due to a pressure head gradient (m·s−1).The migration of liquid water under potential gradient in frozen soil is similar to that in unfrozen unsaturated soil based on the Harlan model and can be described by Richard equation [48, 49]. Then, it is mainly influenced by the factors including water potential gradient and hydraulic conductivity of soils. The values of water potential gradient and hydraulic conductivity in unfrozen and frozen soils differ considerably, while the regulation of water migration can still be assumed to Darcy’s law [50]. Only considering the potential gradient, the flux density of liquid water can be expressed as follows [51, 52]:(2)qlh=−Klh∇h+y,where y is the vertical coordinate (m); h is the pressure head (m); and Klh is the water conductivity coefficient of liquid water under the action of soil-water potential gradient (m·s−1). Then, considering the equivalent volume water content θ as a dependent variable, the mass conservation equation of liquid water in unfrozen and frozen soils can be written as follows [53]:(3)∂θ∂t=∇·Klh∇h+y.The relationship between the amount of water and energy in the soil can be reflected by soil-water characteristic curve (SWCC). Then, SWCC also can represent the relationship between matric, volume water content, and saturation in frozen and unfrozen soils [52]. In this study, the van Genuchten model and Mualem model are adopted to describe the hydraulic properties of unsaturated freeze-thaw soil [53, 54]. The hydraulic properties can be written as follows:(4)h=−Se1/−m−11/nα,(5)Se=θl−θrθs−θr=11+αhnm,h<0,1,h≥0,(6)K=KsSel1−1−Se1/mm2,h<0,Ks,h≥0,where Se is the effective saturation; Ks is the saturation water conductivity coefficient (m·s−1); θl, θs, and θr are liquid water content, saturated liquid water content, and residual water content (m3·m−3); α is the derivative of the soil intake value (m−1); m = 1–1/n; and n and l are experience parameters, and Mualem suggested that l could be determined as 0.5 [54].


## 2.2. Heat Transfer

Compared with the heat conduction with phase change of ice to water, the energy released by the heat convection was very small and could be neglected during the heat transfer analysis of freeze-thaw soil [55–57]. Then, the governing equations of heat transfer in freeze-thaw soil can be written as follows [49]:(7)Cm∂T∂t=∂∂xλm∂T∂x+∂∂yλm∂T∂y,where Cm is the equivalent volume heat capacity and λm is the equivalent thermal conductivity. According to the method of sensible heat capacity, equivalent volume heat capacity Cm and equivalent thermal conductivity λm in freeze-thaw soils can be written as follows [55, 58]:(8)Cm=Cf,T<Tf−ΔT,Cu+Cf2+Ls2ΔT,Tf−ΔT≤T≤Tf+ΔT,Cu,T>Tf+ΔT,(9)λm=λf,T<Tf−ΔT,λu−λf2T−Tm−ΔT+λf,Tf−ΔT≤T≤Tf+ΔT,λu,T>Tf+ΔT,where Tf ± ∆T is the temperature range of phase change; Cu and λu are the volume heat capacity and thermal conductivity of unfrozen soil; Cf and λf are the volume heat capacity and thermal conductivity of frozen soil; and Ls is the latent heat of phase change in per unit volume. While the temperature ranges of phase change in different soil are difficult to gain, the heat transfer equation adopted in this study is written as follows [49]:(10)C∂T∂t=∂∂xλ∂T∂x+∂∂yλ∂T∂y+L∂θi∂tρi,where C is the volume heat capacity; λ is the thermal conductivity; and L is the latent heat of freezing of liquid water (approximately 3.34 × 105 J·kg−1).After soil freezing, not all liquid water can be transferred to solid ice. The rest of liquid water is called as unfrozen water. Then, the unfrozen water content is related to factors including the temperature, pressure, water salinity, mineralogy, soil specific surface area, and soil surface chemistry [59–62]. Based on previous work and related theories, the empirical expression as follows is used to determine the maximum unfrozen water content in the freezing process [49, 50, 52]:(11)θumax=aT−273.15−b,where a and b are parameters related to the soil properties. The volume liquid water content can be determined by temperature and equivalent volume [50, 52]:(12)θl=θ,others,θumax,T<Tfandθ>θumax,where Tf is the freezing point of saturated soil (0°C). Previous studies showed that the freezing point of soil is not a fixed value, and the ice grows only when the water content exceeds θumax [63].


## 3. A Drainage Channel Built in High-Altitude Permafrost Environment

### 3.1. Study Area and Drainage Channel

In this study, the drainage channel built for Yanhu Lake, in the Hoh Xil region on the QTP, was taken as an example. In the area, the elevation ranges from 4460 to 4500 m a.s.l. The meteorological data were collected from the Wu Daoliang weather station, which is about 70 km far away from the channel. The collected data showed that, in the period from 1965 to 2019, the main annual air temperature ranged from −6.5 to 3.7°C and the main annual precipitation ranged from 136 to 480 mm. Since the 1980s, the climate warming and wetting at the region was significant.A standard whether station was set up at the channel in 2020. Figure1 shows the observed air temperatures during the monitoring period from 24 May to 10 November in 2020. In the period, the maximum air temperature reached to 12°C, which appeared in mid-August. Combined with the collected air temperatures from the Wu Daoliang weather station, the air temperature at the projected area can be fitted using a sinusoidal function as follows:(13)To=−3.02+11.6sin2πt365+0.3,where t is the time (d).Figure 1
Observed air temperatures at the channel from May 24 to November 10, 2020.


### 3.2. Ground Temperature Observation at the Channel

To observe thermal regimes of permafrost subgrade beneath the channel, three ground temperature boreholes including natural ground (NG), slope shoulder (SS), and slope toe (ST) of the channel were drilled after excavation of the channel (Figure2). The three boreholes were 20 m in depth. Considering the symmetry of the channel, the three boreholes were all arranged at one side of the channel. Within each borehole, 32 ground temperature sensors were installed with a spacing of 0.5 m between 0–10 m depth and of 1 m below 10 m depth. The measurement range of the sensor is from −40 to 40°C, with an accuracy of ±0.05°C. All the sensors were connected to a CR6 data collector, and the observation frequency of ground temperatures was 6 hours.Figure 2
Ground temperature boreholes drilled at one side of the channel (unit: m).


## 3.1. Study Area and Drainage Channel

In this study, the drainage channel built for Yanhu Lake, in the Hoh Xil region on the QTP, was taken as an example. In the area, the elevation ranges from 4460 to 4500 m a.s.l. The meteorological data were collected from the Wu Daoliang weather station, which is about 70 km far away from the channel. The collected data showed that, in the period from 1965 to 2019, the main annual air temperature ranged from −6.5 to 3.7°C and the main annual precipitation ranged from 136 to 480 mm. Since the 1980s, the climate warming and wetting at the region was significant.A standard whether station was set up at the channel in 2020. Figure1 shows the observed air temperatures during the monitoring period from 24 May to 10 November in 2020. In the period, the maximum air temperature reached to 12°C, which appeared in mid-August. Combined with the collected air temperatures from the Wu Daoliang weather station, the air temperature at the projected area can be fitted using a sinusoidal function as follows:(13)To=−3.02+11.6sin2πt365+0.3,where t is the time (d).Figure 1
Observed air temperatures at the channel from May 24 to November 10, 2020.


## 3.2. Ground Temperature Observation at the Channel

To observe thermal regimes of permafrost subgrade beneath the channel, three ground temperature boreholes including natural ground (NG), slope shoulder (SS), and slope toe (ST) of the channel were drilled after excavation of the channel (Figure2). The three boreholes were 20 m in depth. Considering the symmetry of the channel, the three boreholes were all arranged at one side of the channel. Within each borehole, 32 ground temperature sensors were installed with a spacing of 0.5 m between 0–10 m depth and of 1 m below 10 m depth. The measurement range of the sensor is from −40 to 40°C, with an accuracy of ±0.05°C. All the sensors were connected to a CR6 data collector, and the observation frequency of ground temperatures was 6 hours.Figure 2
Ground temperature boreholes drilled at one side of the channel (unit: m).


## 4. Numerical Simulations and Results Analysis

### 4.1. Computational Model

Previous studies showed that the boundary error would be less than 10% when the computational domain is 3–5 times of the equivalent diameter of the model [64]. In this study, the computational model was constructed based on the actual dimension of the channel. The depth of the channel was 4.7 m, and the bottom width was set as 25 m. The slope of the channel was released by 1 : 3. At present, the water depth within the channel was determined as 1 m throughout one year based on the field survey. The computation on thermal regimes of subgrade under the channel can be simplified as a 2D unsteady heat transfer problem.According to the geological survey conducted before the drainage excavation, the soil strata can be simplified as sandy gravel, silty clay, and weathered mudstone from the surface down. The thermal and hydraulic parameters of the three soil layers are listed in Table1, which were gained based on the borehole drillings and related laboratory tests [34, 49, 50, 52, 65].Table 1
Thermal and hydraulic parameters of the three soil layers.
ParametersλuλfCuCfABW·m−1 °C−1W·m−1 °C−1J·m−3°C−1J·m−3 °C−1Sandy gravel1.912.612.41 × 1061.86 × 10610.670.57Silty clay1.422.123.35 × 1062.54 × 1066.90.47Mudstone1.471.822.09 × 1061.84 × 1069.30.52ParametersαθrθsKsΡLsm−1m·s−1kg·m−3J·m−3Sandy gravel3.280.010.442.4 × 10−718002.31 × 107Silty clay2.60.020.353.3 × 10−816006.51 × 107Mudstone2.30.020.251.2 × 10−817003.77 × 107


### 4.2. Boundary Conditions

According to the meteorological data from the Wu Daoliang weather station, the long-term air temperatures in the study area showed a linear increase with the rate of 0.33°C/10a. Then, based on the adherent layer theory [62], the upper boundaries of GF and FE were set as temperature boundary:(14)T=−0.502+11.6sin2πt365+0.3+0.33t3650,where t is the time (d).As an alternative to heat budget models, linear correlations between air and water temperatures during open water periods were studied widely in previous studies [66]. Based on the analysis of a stream’s heat transport/budget equation, the air temperature shows a good linear relationship with water temperature. The relationship between the water temperature in the lake and the air temperature can be linearly fitted as follows:(15)Twt=A+BTat,where Twt and Ta(t) are the water temperature and the air temperature in the same time scale; A and B are constants; and t is the time (d).According to the observed water temperatures (equation (15)) in lakes near the study area [40] and the air temperature showed as equation (16), the two constants of A and B in equation (14) used in this study were determined as follows:(16)Tw1=5.96+5sin2πt365−0.65,(17)Ta1=−4.02+11.1sin2πt365−8.18,where t is the time (d).Without considering the impacts of artificial heat inputs, ground water inputs, stream shading, and wind sheltering, the values ofA and B in equation (14) in the same area have very slight differences [66]. Then, the water temperature in the channel can be written as equation (17). Thus, the boundary BH and BG were set as follows:(18)Tw=4.72+5.21sin2πt365+0.3+0.15t3650,where t is the time (d).The boundary condition at CD was set as heat flux with the value of 0.03 W/m2, which was gained from the geothermal gradient within the ground temperature boreholes. The boundary of DE was thermal insulation boundary, and that of BC was the symmetry boundary. With the governing equations and boundary conditions above, the problem was solved numerically using the commercial software of COMSOL Multiphysics. The spatial and temporal discretization of governing equations was carried out by using the finite method. The simulation was conducted over a time period of 50 years before the drainage channel excavation to gain the initial temperature field. After the excavation, the boundary conditions were set as the ones described above.


### 4.3. Model Validation

A comparison between field observed and numerical simulated ground temperatures at ST on October 15 in the first year after the channel excavation was conducted to validate the numerical model and its parameterization (Figure3). It can be found that the numerical simulated results agreed with the field observed well within the permafrost layer, including the depth of the permafrost table and the permafrost temperatures below. While in the active layer, there were some slight discrepancies between the numerical simulated and field observed ground temperatures. The discrepancies may relate to the simplification of the soil strata in the numerical simulation and the depth of the temperature sensors installed within the borehole. Overall, the numerical model and its parameterization in this study can be used to simulate the permafrost thermal regimes beneath the channel.Figure 3
Field observed and numerical simulated ground temperatures vs depth at ST on October 15, 1st year after the excavation of the channel.


## 4.1. Computational Model

Previous studies showed that the boundary error would be less than 10% when the computational domain is 3–5 times of the equivalent diameter of the model [64]. In this study, the computational model was constructed based on the actual dimension of the channel. The depth of the channel was 4.7 m, and the bottom width was set as 25 m. The slope of the channel was released by 1 : 3. At present, the water depth within the channel was determined as 1 m throughout one year based on the field survey. The computation on thermal regimes of subgrade under the channel can be simplified as a 2D unsteady heat transfer problem.According to the geological survey conducted before the drainage excavation, the soil strata can be simplified as sandy gravel, silty clay, and weathered mudstone from the surface down. The thermal and hydraulic parameters of the three soil layers are listed in Table1, which were gained based on the borehole drillings and related laboratory tests [34, 49, 50, 52, 65].Table 1
Thermal and hydraulic parameters of the three soil layers.
ParametersλuλfCuCfABW·m−1 °C−1W·m−1 °C−1J·m−3°C−1J·m−3 °C−1Sandy gravel1.912.612.41 × 1061.86 × 10610.670.57Silty clay1.422.123.35 × 1062.54 × 1066.90.47Mudstone1.471.822.09 × 1061.84 × 1069.30.52ParametersαθrθsKsΡLsm−1m·s−1kg·m−3J·m−3Sandy gravel3.280.010.442.4 × 10−718002.31 × 107Silty clay2.60.020.353.3 × 10−816006.51 × 107Mudstone2.30.020.251.2 × 10−817003.77 × 107


## 4.2. Boundary Conditions

According to the meteorological data from the Wu Daoliang weather station, the long-term air temperatures in the study area showed a linear increase with the rate of 0.33°C/10a. Then, based on the adherent layer theory [62], the upper boundaries of GF and FE were set as temperature boundary:(14)T=−0.502+11.6sin2πt365+0.3+0.33t3650,where t is the time (d).As an alternative to heat budget models, linear correlations between air and water temperatures during open water periods were studied widely in previous studies [66]. Based on the analysis of a stream’s heat transport/budget equation, the air temperature shows a good linear relationship with water temperature. The relationship between the water temperature in the lake and the air temperature can be linearly fitted as follows:(15)Twt=A+BTat,where Twt and Ta(t) are the water temperature and the air temperature in the same time scale; A and B are constants; and t is the time (d).According to the observed water temperatures (equation (15)) in lakes near the study area [40] and the air temperature showed as equation (16), the two constants of A and B in equation (14) used in this study were determined as follows:(16)Tw1=5.96+5sin2πt365−0.65,(17)Ta1=−4.02+11.1sin2πt365−8.18,where t is the time (d).Without considering the impacts of artificial heat inputs, ground water inputs, stream shading, and wind sheltering, the values ofA and B in equation (14) in the same area have very slight differences [66]. Then, the water temperature in the channel can be written as equation (17). Thus, the boundary BH and BG were set as follows:(18)Tw=4.72+5.21sin2πt365+0.3+0.15t3650,where t is the time (d).The boundary condition at CD was set as heat flux with the value of 0.03 W/m2, which was gained from the geothermal gradient within the ground temperature boreholes. The boundary of DE was thermal insulation boundary, and that of BC was the symmetry boundary. With the governing equations and boundary conditions above, the problem was solved numerically using the commercial software of COMSOL Multiphysics. The spatial and temporal discretization of governing equations was carried out by using the finite method. The simulation was conducted over a time period of 50 years before the drainage channel excavation to gain the initial temperature field. After the excavation, the boundary conditions were set as the ones described above.


## 4.3. Model Validation

A comparison between field observed and numerical simulated ground temperatures at ST on October 15 in the first year after the channel excavation was conducted to validate the numerical model and its parameterization (Figure3). It can be found that the numerical simulated results agreed with the field observed well within the permafrost layer, including the depth of the permafrost table and the permafrost temperatures below. While in the active layer, there were some slight discrepancies between the numerical simulated and field observed ground temperatures. The discrepancies may relate to the simplification of the soil strata in the numerical simulation and the depth of the temperature sensors installed within the borehole. Overall, the numerical model and its parameterization in this study can be used to simulate the permafrost thermal regimes beneath the channel.Figure 3
Field observed and numerical simulated ground temperatures vs depth at ST on October 15, 1st year after the excavation of the channel.


## 5. Results and Analysis

### 5.1. Long-Term Permafrost Thermal Regimes beneath the Channel in Warm Seasons

In the permafrost regions on the QTP, the maximum seasonal thaw depth generally occurs in mid-October. In the following analysis, the time point was chose to investigate the long-term permafrost thermal regimes beneath the channel. For brevity, we only took the channel with the width of 15 m as an example.Figure4 shows the thermal regimes of subgrade beneath the channel on October 15 within the 50 years after the channel excavation. It can be seen that, due to combined thermal effects of the flowing water and climate warming, the permafrost beneath the channel degraded rapidly with time went on. While in the natural ground far away from the channel, the permafrost degraded slowly due to the climate warming effect. In the 5th year, the maximum seasonal thaw depth at NG, ST, and centerline of the channel (CC) was 2.4, 9.2, and 10.5 m, respectively. A thaw bulb developed beneath the channel due to the thermal effect of flow water. The permafrost temperatures beneath the channel also increased considerably comparing with those beneath the natural ground at the same depths. For example, at the depth of 10 m, the permafrost temperatures under CC and NG were 0.14°C and −0.45°C, respectively, with a difference up to 0.6°C. With increase in the operation time, the thaw bulb beneath the channel expanded rapidly in the vertical direction. In 10th year after the channel excavation, the maximum seasonal thawing depth at CC and ST reached to 16.7 and 14.4 m, respectively. Till the 30th year, the permafrost beneath the channel degraded almost totally within the depth of 30 m. After that, the later thermal erosion of the channel became considerable, and the permafrost beneath the channel slope degraded with the depth of 30 m.Figure 4
Thermal regimes of subgrade beneath the channel on October 15 after 5 (a), 10 (b), 20 (c), 30 (d), 40 (e), and 50 (f) years of the channel excavation (unit:°C).
(a)(b)(c)(d)(e)(f)The above results showed that, during the first 30 years of the channel operation, the permafrost beneath the channel experienced a rapid downward degradation because of the significant thermal effect of the flowing water. Along with the rapid downward degradation of permafrost beneath the channel, a large lateral thermal gradient developed beneath the slope of the channel. This significant thermal gradient would cause uneven deformations of the slope and even slumps or collapse because of the temperature dependence of mechanical properties of frozen soils. Thus, how to ensure the slope stability of the channel is a great challenge in permafrost environment.


### 5.2. Long-Term Permafrost Thermal Regimes beneath the Channel in Cold Seasons

In mid-April, the active layer on the QTP generally refreezes completely. In order to investigate the thermal regimes of subgrade beneath the channel in cold seasons, ground temperature distributions within the subgrade on April 15 within the 50 years after the channel excavation were analyzed in this section.Figure5 shows the thermal regimes of subgrade beneath the channel on April 15 within the 50 years after the channel excavation. The shallow ground far away from the slope shoulder of the channel refrozen totally in cold seasons. While, beneath the channel, the thaw bulb beneath the channel did not freeze due to the flowing water. Compared with the last warm season in the same year (Figure 4), the thaw bulb still expanded in the following cold seasons. In the 5th year, the thaw bulb mainly exited vertically beneath the bottom of the channel and the water body within the channel. With the time went on, the thaw bulb developed rapidly both in vertical and lateral directions. Special attention should be paid to ground thermal regimes beneath the slope of the channel. The shallow ground on the slope was frozen but the deep ground was thawed in the cold season due to the lateral thermal erosion of the flowing water. Till the 50th year, the thaw bulb in the cold season reached to the natural ground about 10 m far away from the slope shoulder of the channel. This meant that, all the deep ground beneath the slope of the channel was thawed throughout the year.Figure 5
Thermal regimes of subgrade beneath the channel on April 15, after 5 (a), 10 (b), 20 (c), 30 (d), 40 (e), and 50 (f) years of the channel excavation (unit: °C).
(a)(b)(c)(d)(e)(f)


### 5.3. Impacts of Channel Width on Long-Term Thermal Regimes of Permafrost Subgrade

The width of the channel is not only related to the slope stability of the channel but also the permafrost degradation beneath the channel. To investigate the impacts of the channel width, a series of numerical simulations with different widths of the channel (15, 25, and 35 m) were carried out in this study. The flow of water was determined as a constant, and then the water depths were different for the channels with different widths. With the channel widths of 15, 25, and 35 m, the water depths were set as 1.53, 1.0, and 0.74 m, respectively.Figure6 shows ground temperature profiles at CC, ST, mid-slope (MS), and SS of the channels with the widths of 15, 25, and 35 m on October 15 in the 15th year after the channel excavation. It can be seen that, with different widths of the channel, the long-term thermal regimes of the subgrade differed, but the magnitudes of the difference were different at the four locations. At CC (Figure 6(a)), the ground thermal regimes for the channels with different widths were close to each other. The difference in the maximum thaw depths for the three cases was about 1 m. This means that, with different widths, the thermal erosion from the channels did not differed considerably at this location. While, at ST, the maximum thaw depths for the channels with different widths differed considerably. The smaller the width of the channel was, the grater the maximum thaw depth was at ST (Figure 6(b)). This considerable difference was related to different water depths within the channels with different widths. When the width of the channel increased from 15 to 35 m, the maximum thaw depth at ST decreased from 18.2 to 15.5 m. At MS (Figure 6(c)), the difference in maximum thaw depths for the three cases was also considerable. The wider the channel was, the smaller the maximum thaw depths were. At SS (Figure 6(d)), the ground thermal regimes were also very close to each other with three widths of the channel.Figure 6
Ground temperature profiles at CC (a), ST (b), MS (c), and SS (d) of the channels with widths of 15, 25, and 35 m on October 15 in the 10th year after the channel excavation.
(a)(b)(c)(d)The above results showed that the thermal regimes of subgrade differed considerably for the channels with different widths. The difference mainly existed at slope toe and middle of the slope. While at the centerline and slope shoulder, the ground thermal regimes were very close to each other. For slope stability, the greater thaw depth is generally harmful. Thus, it can be concluded that a channel built on permafrost should be wide-and-shallow rather than narrow-and-deep if conditions permit.


## 5.1. Long-Term Permafrost Thermal Regimes beneath the Channel in Warm Seasons

In the permafrost regions on the QTP, the maximum seasonal thaw depth generally occurs in mid-October. In the following analysis, the time point was chose to investigate the long-term permafrost thermal regimes beneath the channel. For brevity, we only took the channel with the width of 15 m as an example.Figure4 shows the thermal regimes of subgrade beneath the channel on October 15 within the 50 years after the channel excavation. It can be seen that, due to combined thermal effects of the flowing water and climate warming, the permafrost beneath the channel degraded rapidly with time went on. While in the natural ground far away from the channel, the permafrost degraded slowly due to the climate warming effect. In the 5th year, the maximum seasonal thaw depth at NG, ST, and centerline of the channel (CC) was 2.4, 9.2, and 10.5 m, respectively. A thaw bulb developed beneath the channel due to the thermal effect of flow water. The permafrost temperatures beneath the channel also increased considerably comparing with those beneath the natural ground at the same depths. For example, at the depth of 10 m, the permafrost temperatures under CC and NG were 0.14°C and −0.45°C, respectively, with a difference up to 0.6°C. With increase in the operation time, the thaw bulb beneath the channel expanded rapidly in the vertical direction. In 10th year after the channel excavation, the maximum seasonal thawing depth at CC and ST reached to 16.7 and 14.4 m, respectively. Till the 30th year, the permafrost beneath the channel degraded almost totally within the depth of 30 m. After that, the later thermal erosion of the channel became considerable, and the permafrost beneath the channel slope degraded with the depth of 30 m.Figure 4
Thermal regimes of subgrade beneath the channel on October 15 after 5 (a), 10 (b), 20 (c), 30 (d), 40 (e), and 50 (f) years of the channel excavation (unit:°C).
(a)(b)(c)(d)(e)(f)The above results showed that, during the first 30 years of the channel operation, the permafrost beneath the channel experienced a rapid downward degradation because of the significant thermal effect of the flowing water. Along with the rapid downward degradation of permafrost beneath the channel, a large lateral thermal gradient developed beneath the slope of the channel. This significant thermal gradient would cause uneven deformations of the slope and even slumps or collapse because of the temperature dependence of mechanical properties of frozen soils. Thus, how to ensure the slope stability of the channel is a great challenge in permafrost environment.


## 5.2. Long-Term Permafrost Thermal Regimes beneath the Channel in Cold Seasons

In mid-April, the active layer on the QTP generally refreezes completely. In order to investigate the thermal regimes of subgrade beneath the channel in cold seasons, ground temperature distributions within the subgrade on April 15 within the 50 years after the channel excavation were analyzed in this section.Figure5 shows the thermal regimes of subgrade beneath the channel on April 15 within the 50 years after the channel excavation. The shallow ground far away from the slope shoulder of the channel refrozen totally in cold seasons. While, beneath the channel, the thaw bulb beneath the channel did not freeze due to the flowing water. Compared with the last warm season in the same year (Figure 4), the thaw bulb still expanded in the following cold seasons. In the 5th year, the thaw bulb mainly exited vertically beneath the bottom of the channel and the water body within the channel. With the time went on, the thaw bulb developed rapidly both in vertical and lateral directions. Special attention should be paid to ground thermal regimes beneath the slope of the channel. The shallow ground on the slope was frozen but the deep ground was thawed in the cold season due to the lateral thermal erosion of the flowing water. Till the 50th year, the thaw bulb in the cold season reached to the natural ground about 10 m far away from the slope shoulder of the channel. This meant that, all the deep ground beneath the slope of the channel was thawed throughout the year.Figure 5
Thermal regimes of subgrade beneath the channel on April 15, after 5 (a), 10 (b), 20 (c), 30 (d), 40 (e), and 50 (f) years of the channel excavation (unit: °C).
(a)(b)(c)(d)(e)(f)


## 5.3. Impacts of Channel Width on Long-Term Thermal Regimes of Permafrost Subgrade

The width of the channel is not only related to the slope stability of the channel but also the permafrost degradation beneath the channel. To investigate the impacts of the channel width, a series of numerical simulations with different widths of the channel (15, 25, and 35 m) were carried out in this study. The flow of water was determined as a constant, and then the water depths were different for the channels with different widths. With the channel widths of 15, 25, and 35 m, the water depths were set as 1.53, 1.0, and 0.74 m, respectively.Figure6 shows ground temperature profiles at CC, ST, mid-slope (MS), and SS of the channels with the widths of 15, 25, and 35 m on October 15 in the 15th year after the channel excavation. It can be seen that, with different widths of the channel, the long-term thermal regimes of the subgrade differed, but the magnitudes of the difference were different at the four locations. At CC (Figure 6(a)), the ground thermal regimes for the channels with different widths were close to each other. The difference in the maximum thaw depths for the three cases was about 1 m. This means that, with different widths, the thermal erosion from the channels did not differed considerably at this location. While, at ST, the maximum thaw depths for the channels with different widths differed considerably. The smaller the width of the channel was, the grater the maximum thaw depth was at ST (Figure 6(b)). This considerable difference was related to different water depths within the channels with different widths. When the width of the channel increased from 15 to 35 m, the maximum thaw depth at ST decreased from 18.2 to 15.5 m. At MS (Figure 6(c)), the difference in maximum thaw depths for the three cases was also considerable. The wider the channel was, the smaller the maximum thaw depths were. At SS (Figure 6(d)), the ground thermal regimes were also very close to each other with three widths of the channel.Figure 6
Ground temperature profiles at CC (a), ST (b), MS (c), and SS (d) of the channels with widths of 15, 25, and 35 m on October 15 in the 10th year after the channel excavation.
(a)(b)(c)(d)The above results showed that the thermal regimes of subgrade differed considerably for the channels with different widths. The difference mainly existed at slope toe and middle of the slope. While at the centerline and slope shoulder, the ground thermal regimes were very close to each other. For slope stability, the greater thaw depth is generally harmful. Thus, it can be concluded that a channel built on permafrost should be wide-and-shallow rather than narrow-and-deep if conditions permit.


## 6. Conclusions

In permafrost regions, construction of a channel involves a large amount of excavation activities and changes to surface water body, which can exert great impacts on thermal regimes of underlying permafrost. In this paper, a coupled mathematical model of heat and moisture transfer for freeze-thaw soil was constructed to investigate the long-term thermal regimes of subgrade beneath a drainage channel built for an expanding lake on the Qinghai-Tibet Plateau. Using numerical simulations, thermal regimes of the subgrade both in warm and cold seasons were analyzed during a 50-year period, as well as the impact of the channel width. The conclusions were obtained as follows:(1)
In permafrost environment, excavation of the channel and the flowing water within could lead to a significant permafrost degradation. During the first 30 years of the channel operation, the permafrost beneath the channel mainly experienced a rapid downward degradation due to the thermal effect of the flowing water. After that, the lateral thermal erosion of the flowing water caused a rapid permafrost degradation beneath the slope of the channel. With a width of 15 m, the permafrost beneath the channel bottom and slope would degrade totally within a depth of 30 m in a 50-year period.(2)
The ground beneath the channel would not refreeze in cold seasons due to the flowing water within the channel, and the thaw bulb developed throughout a year. During the first 10 years, the thaw bulb mainly existed vertically beneath the channel. After that, the thaw bulb expanded quickly in lateral direction. Till the 50th year after the excavation of the channel, the thaw bulb reached to the natural ground about 25∼30 m far away from the centerline of the channel with a width of 15 m.(3)
In permafrost environment, the width of the channel is an important factor. With different widths, the long-term thermal regimes of the subgrade beneath the channels differed considerably. The maximum difference was at the slope toe of the channel. The narrower the channel was, the larger the maximum thaw depth was at the slope toe. Thus, it is recommended that the channels in permafrost environment should be designed as wide-and-shallow rather than narrow-and-deep if conditions permit.



---

*Source: 6613114-2021-02-12.xml*