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

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. +e 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 15m, the thaw bulb beneath the channel could expand laterally to the natural ground about 10m 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.


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]. e 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][7][8].
e 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][12][13][14]. e 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. e 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][19][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 tailend-lake [21][22][23][24]. From 2011 to 2019, the surface area of the Zonag Lake decreased from 269 to 150 km 2 , while the surface area of Yanhu Lake increased by almost 3 times from 73 to 209 km 2 [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][27][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][31][32][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]. us, 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][36][37][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. rough 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.
en, a numerical model was constructed based on the engineering geotechnical investigation conducted in the field and the dimensions of the drainage channel. e numerical model was validated with field observations carried out at the drainage channel. en, 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.

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: where θ is the equivalent volume of water content (m 3 ·m −3 ); θ u is the volume of unfrozen water content (m 3 ·m −3 ); θ i is the volume of ice content (m 3 ·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 q lh is the liquid water flux density which can represent liquid flows due to a pressure head gradient (m·s −1 ). e 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]. en, it is mainly influenced by the factors including water potential gradient and hydraulic conductivity of soils. e 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]: where y is the vertical coordinate (m); h is the pressure head (m); and K lh is the water conductivity coefficient of liquid water under the action of soil-water potential gradient (m·s −1 ). en, 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]: e relationship between the amount of water and energy in the soil can be reflected by soil-water characteristic curve (SWCC). en, 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]. e hydraulic properties can be written as follows: where S e is the effective saturation; K s 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 (m 3 ·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].

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][56][57]. en, the governing equations of heat transfer in freeze-thaw soil can be written as follows [49]: where C m 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 C m and equivalent thermal conductivity λ m in freeze-thaw soils can be written as follows [55,58]: where T f ± ∆T is the temperature range of phase change; C u and λ u are the volume heat capacity and thermal conductivity of unfrozen soil; C f and λ f are the volume heat capacity and thermal conductivity of frozen soil; and L s 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]: Advances in Materials Science and Engineering 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 × 10 5 J·kg −1 ). After soil freezing, not all liquid water can be transferred to solid ice. e rest of liquid water is called as unfrozen water. en, the unfrozen water content is related to factors including the temperature, pressure, water salinity, mineralogy, soil specific surface area, and soil surface chemistry [59][60][61][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]: where a and b are parameters related to the soil properties. e volume liquid water content can be determined by temperature and equivalent volume [50,52]: where T f 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].

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. e meteorological data were collected from the Wu Daoliang weather station, which is about 70 km far away from the channel. e 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. Figure 1 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: T o � −3.02 + 11.6 sin 2πt 365 where t is the time (d).

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 (Figure 2). e 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. e 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.

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. e depth of the channel was 4.7 m, and the bottom width was set as 25 m. e 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. e 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. e thermal and hydraulic parameters of the three soil layers are listed in Table 1, which were gained based on the borehole drillings and related laboratory tests [34,49,50,52,65].

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. en, based on the adherent layer theory [62], the upper boundaries of GF and FE were set as temperature boundary: T � −0.502 + 11.6 sin 2πt 365 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. e relationship between the water temperature in the lake and the air temperature can be linearly fitted as follows: where T w (t) and T a (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). Parameters

Advances in Materials Science and Engineering
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: T a1 � −4.02 + 11.1 sin 2πt where t is the time (d).
Without considering the impacts of artificial heat inputs, ground water inputs, stream shading, and wind sheltering, the values of A and B in equation (14) in the same area have very slight differences [66]. en, the water temperature in the channel can be written as equation (17).
us, the boundary BH and BG were set as follows: where t is the time (d). e boundary condition at CD was set as heat flux with the value of 0.03 W/m 2 , which was gained from the geothermal gradient within the ground temperature boreholes. e 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. e spatial and temporal discretization of governing equations was carried out by using the finite method. e 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.

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 ( Figure 3). 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. e 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.

Long-Term Permafrost
ermal 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. Figure 4 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 5 th 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. e 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 10 th 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 30 th 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.
e 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. is 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. us, how to ensure the slope stability of the channel is a great challenge in permafrost environment.

Long-Term Permafrost
ermal 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. Figure 5 shows the thermal regimes of subgrade beneath the channel on April 15 within the 50 years after the channel excavation.
e 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 5 th year, the thaw bulb mainly exited vertically beneath the 6 Advances in Materials Science and Engineering  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. e 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 50 th 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. is meant that, all the deep ground beneath the slope of the channel was thawed throughout the year.

Impacts of Channel Width on Long-Term ermal Regimes of Permafrost
Subgrade. e 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. e 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. Figure 6 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 15 th 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. e difference in the maximum thaw depths for the three cases was about 1 m. is 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. e smaller the width of the channel was, the grater the maximum thaw depth was at ST (Figure 6(b)). is 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. e 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.
e above results showed that the thermal regimes of subgrade differed considerably for the channels with different widths. e 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  Advances in Materials Science and Engineering generally harmful. us, it can be concluded that a channel built on permafrost should be wide-and-shallow rather than narrow-and-deep if conditions permit.

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. e 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) e 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. e maximum difference was at the slope toe of the channel. e narrower the channel was, the larger the maximum thaw depth was at the slope toe. us, it is recommended that the channels in permafrost environment should be designed as wide-and-shallow rather than narrow-and-deep if conditions permit.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest related to this manuscript.