1. Introduction
Permafrost is mainly found at high latitudes and high altitudes, such as Siberia, the Daxingan Mountains, the Tibetan Plateau, Alaska, etc. The area of permafrost in China is second only to those of Russia and Canada [
1,
2]. The total area of permafrost in the Qinghai–Tibetan Plateau is equivalent to about 70% of the total area of permafrost in China [
2,
3]. Compared with the relatively narrow railways and expressways built in the permafrost region, the wide expressways have essential differences in terms of heat sources, heat transfer processes, and heat transfer intensity [
4]. Doubling the width of the subgrade was shown to result in a 60% increase in the average annual heat flux of the subgrade bottom, and the heat was mainly concentrated in the permafrost under the road center [
5]. The subgrade width strongly influences the residual thawed layer of the permafrost subgrade. To achieve the goal where the underlying permafrost has no residual thawed layer in the 20th year after subgrade construction, the subgrade width must be less than 6 m [
6]. The construction of an expressway may accelerate the degradation of underlying permafrost. The protection of permafrost subgrades essentially involves using certain engineering measures to reduce or eliminate the heat transfer through the permafrost subgrade into the underlying permafrost, thereby achieving permafrost subgrade energy balance [
7]. Therefore, against the background of global warming, it is necessary to take measures to protect wide subgrade expressways and the permafrost underneath the subgrade [
8,
9,
10,
11,
12,
13].
Researchers have proposed a series of prevention and control measures for engineering issues in permafrost [
14,
15,
16]. For example, adding an awning to road slopes can significantly uplift the artificial permafrost table and reduce the temperature of the embankment body [
17,
18]. However, due to the harsh meteorological conditions on the Qinghai–Tibet Plateau, awning embankments have been found to be damaged by strong winds [
19]. A two-phase closed thermosyphon (TPCT) is considered an effective measure with which to protect permafrost [
20,
21,
22]. However, due to a non-uniform linear heat source, the strong water accumulation, and the ice formation effects of the TPCT, road surface cracks often appear on TPCT embankments [
23,
24]. A crushed rock embankment (CRE) can reduce the ground temperature via forced convection and stacking effects through the packing gravel layer [
25,
26,
27,
28], which significantly reduces the temperature of the embankment and ensures the long-term thermal stability of the underlying permafrost. However, it has also been found that the CRE is highly influenced by the environment, and the CRE makes little contribution to the permafrost stability in a warm permafrost region [
29,
30,
31]. Furthermore, a ventilation structure is also an effective means of disposal for permafrost subgrade. With the aid of forced convection in a ventilation pipe, a duct-ventilated embankment can effectively stabilize the underlying permafrost, especially for an embankment with a ventilation pipe installed in a lower position [
32,
33,
34,
35].
The above studies regarding the disposal measures of permafrost subgrades focus on the effects and parameter optimization of the protection measures, and there is currently no research that has combined the heat transfer process of permafrost subgrade with its heat control parameters. In this work, through an analysis of the spatial–temporal temperature distribution characteristics and attenuation law of the energy flow transfer process of permafrost subgrade, an interface heat control method—which involves importing supposed cold energy at the interface of the permafrost subgrade to achieve permafrost subgrade energy balance—is proposed. The application time, application position, and required cold energy quantity for the interface heat control are given, the effect of the interface heat control is analyzed, and the parameters determined via the interface heat control are applied to a mechanical ventilation subgrade.
2. Earth–Atmosphere-Coupled Calculation Model for Permafrost Subgrade
2.1. Physical Model
An Earth–atmosphere-coupled calculation model of permafrost subgrade, one that considers the upper air environment, was built (
Figure 1). The calculation model included the air, subgrade, and underlying permafrost foundation regions. According to the Qinghai–Tibet expressway design scheme and geological drilling results in the Wudaoliang area, the subgrade was divided into a 0.12 m thick, asphalt layer; a 0.53 m thick, cement-stabilized gravel layer; and a 2.35 m thick, subgrade filling layer. The underlying permafrost foundation region was divided into a 2.0 m thick, silty clay layer; a 5.0 m thick, gravel clay layer; and a 23.0 m thick, weathered mudstone layer. The cross-section of the subgrade adopted an integral subgrade structure, and its overall width was 26.0 m with a 0.75 m earth shoulder on both sides. Furthermore, in order to more realistically reflect the convection heat exchange between the distant air flow and the subgrade, the height of the air zone was set to 40.0 m, and the width of the inlet and outlet air regions was set to 60.0 m [
36].
2.2. Numerical Method
An unsteady numerical model was adopted for the calculation of the complex heat transfer processes of permafrost subgrade. The implicit solver and SIMPLE algorithm were utilized in the numerical model. The standard
model was used for the air region forced convection heat transfer simulation. The heat conduction of the solid region was solved using the equivalent heat capacity method. The governing equations for the fluid zone and the solid zone are given in our prior published work [
36].
The natural ground, asphalt pavement, subgrade slopes, and shoulders were set as the internal heat source wall boundary conditions, and these included the consideration of solar radiation absorption, long-wave radiation dissipation, and the water evaporation of the subgrade surface. The source wall setting method is described in detail in the published work of [
36]. The bottom of the calculation model was defined as the heat flux boundary condition, which was set as 0.06 W/m
2 to reflect the deep geothermal energy [
37]. The left and right sides of the air region were defined as the velocity inlet and outflow boundary conditions, respectively. Considering the wind direction variation of the actual climatic environment [
38], the above boundary conditions were swapped at a specific time as follows: the wind direction was defined as roughly east–west (left to right in
Figure 1) during the period of March to August, and then west–east at other times. The wind velocity at different heights was calculated as follows:
where v
ref is the wind speed at the velocity reference point and
x is the height from the ground. The variation in v
ref is shown in
Table 1.
The left and right sides of the solid region were set as adiabatic boundary conditions. The upper boundary of the calculation model was set as the temperature boundary condition, and temperature was obtained as follows [
38]:
where
T is the annual average air temperature for the calculation condition,
t is the year of the calculation time, and
x is the calculation time (20 August was set as the starting time of the calculation).
In order to analyze the major permafrost types, the above annual average air temperature (
T in Equation (2)) was set to −3.0 °C, −3.5 °C, −4.5 °C, and −5.5 °C in the present work, and the four average annual air temperature settings corresponded to the four work cases (which are specified in
Table 2 [
39]).
Furthermore, the interface heat control method was realized by importing the ideal homogeneously distributed cooling energy into a specific interface. Using the calculation, the cumulative heat fluxes of some of the key interfaces of the permafrost subgrade were analyzed (such as the 0 m, 0.3 m, 0.5 m, 1.0 m, 1.5 m, 2.0 m, and 2.35 m interfaces, as shown in
Figure 1). The import calculations were carried out using the ANSYS Fluent 15.0 software package, and the calculation residuals were set as follows: the convergence criterion of the residuals of the energy equation calculation was set to 10
−6; and the convergence criteria of the residuals of the momentum and continuity equation calculations were both set to 10
−5 (which enabled a better calculation accuracy). The above complex boundary conditions of the numerical model were imported into Fluent via the UDF program. The calculation time step was defined as 4 h, and the total computation time was set as 20 years. The detailed numerical method is described in our previous published work [
36].
3. Determination of the Interface Heat Control Parameters for Permafrost Subgrade
3.1. Determining the Application Period of the Interface Heat Control
Determining the application period of the interface heat control must take into account the heat control efficiency and not cause new engineering issues. The climate of the permafrost region is unique: the warm season is shorter, i.e., from May to August; and the rest of the year represents the cold season.
Figure 2 shows the variations in the thawing depth of the underlying permafrost in the 1st, 5th, 10th, and 20th years after the construction of the subgrade. As shown in the figure, the thawing depth of the underlying permafrost exhibits a trend of increasing and then decreasing during the year. During the warm season, the permafrost subgrade absorbed heat and transferred heat into the underlying permafrost. The permafrost warmed and thawed, and the thawing depth increased. In the early cold season, the heat absorbed by the permafrost subgrade had not yet been fully transferred to the underlying permafrost, which was still in a heat-absorbing state: thus, the thawing depth increased. In the late cold season, the underlying permafrost released heat into the air, the permafrost began to freeze, and the thawing depth decreased. For the four annual average air temperature cases, the thawing depth of the underlying permafrost reached its maximum in December, December, December, and November, respectively.
The refreezing depth of the permafrost was less than the thawing depth. Long-term repeated action will cause a sustained increase in the artificial permafrost table.
Figure 3 shows a comparison of the artificial permafrost table underneath the subgrade with the natural ground table in the 1st, 5th, 10th, and 20th years after the construction of the subgrade. Natural ground refers to the natural ground at a distance of 40 m from the subgrade, and the comparison of the natural ground and the underlying permafrost of the subgrade served as an important reference for the effect of interface heat control. If no measures were taken for the permafrost subgrade for the four annual average temperature cases, the differences between the artificial permafrost table and the natural ground table in the 20th year would be 7.33 m, 6.34 m, 4.98 m, and 2.22 m, respectively. Therefore, it is crucial to apply heat control to the permafrost subgrade.
If cold energy is imported into the subgrade after the permafrost thawing depth reaches its maximum, the mechanism for importing cold energy into the subgrade involves increasing the cold energy reserve against the incoming heat in the next year. However, this could lead to new heave issues. If cold energy is imported into the subgrade before the permafrost thawing depth reaches its maximum for the year, the process of heat transfer to the underlying permafrost can be interrupted, and the heat absorbed by the permafrost subgrade during the warm season will not be fully transferred to the underlying permafrost, thereby protecting the permafrost. Combined with the above analysis of the variation in the thawing depth, the period before the permafrost thawing depth reaches its maximum was divided into two phases, i.e., the warm season and the early cold season. The early cold season is the period after the warm season when the thawing depth is still increasing, and it ends when the thawing depth reaches its maximum for the year. In this work, this time was defined as the lower thaw season, which covered September–December, September–December, September–December, and September–November for the four cases.
Figure 4 shows the temperature difference between the subgrade interface temperature and the air temperature during the warm, cold, and lower thaw seasons for the four cases. A positive value indicates that the interface temperature was higher than the air temperature, and a negative value indicates that the interface temperature was lower than the air temperature. As shown in the figure, the temperature at all subgrade interfaces was higher than the air temperature in the cold season and the lower thaw season, and the temperature difference in the lower thaw season was larger than that in the cold season. During the warm season, the temperatures at the 0 m, 0.3 m, and 0.5 m interfaces were lower than the air temperature, and the temperatures at the other interfaces were higher than the air temperature. However, the temperature differences during the warm season were smaller than those in the cold season and the lower thaw season. Therefore, for optimal interface heat control efficiency, the lower thaw season should be selected as the application period of the interface heat control.
In conclusion, selecting the lower thaw season as the application period for interface heat control can interfere with the process of heat transfer into the underlying permafrost and could uplift the artificial permafrost table; moreover, the lower thaw season is more efficient for interface heat control. Therefore, the lower thaw season was chosen as the application period of the interface heat control.
3.2. Determining the Application Position of the Interface Heat Control
Selection of the appropriate position from which to import cold energy involves offsetting the heat through this interface into the underlying permafrost, and the efficiency of the imported cold energy is influenced by the selection of the application position. Therefore, the interface with the lower heat flux should be selected, and the zone of 0.05 m above and below this interface was used as the application position of the interface heat control.
Figure 5 shows the average heat flux during the lower thaw season at different subgrade interfaces for the four annual average air temperature conditions. As shown in the figure, the interfaces with lower heat flux values were the 0 m interface, 0.3 m interface, 0.5 m interface, and 1.5 m interface. Based on construction factors, the 0.5 m and 1.5 m interfaces were initially selected as the application positions of the interface heat control.
The −4.5 °C case was used as an example for imported cold energy, and the cold energy quantity was based on the annual average heat flux at the interface.
Table 3 shows the annual average heat flux at the 0.5 m interface and the 1.5 m interface for the four cases.
Figure 6 shows the variations in the artificial permafrost table after importing different amounts of cold energy at the 0.5 m and the 1.5 m interfaces, where q
1 is the annual average heat flux at the 0.5 m interface and q
2 is the annual average heat flux at the 1.5 m interface. As shown in the figure, within 20 years after the construction of the subgrade, when the quantity of cold energy has the same multiple as the interface heat flux, the artificial permafrost table after importing cold energy at the 0.5 m interface was found to be lower than the artificial permafrost table after importing cold energy at the 1.5 m interface (and q
1 is lower than q
2). Therefore, after importing the same quantity of cold energy at the 0.5 m interface and the 1.5 interface, the magnitude of the underlying permafrost table uplift for the 0.5 m interface case was found to be larger than the magnitude of the underlying permafrost table uplift for the 1.5 m interface case. To uplift the artificial permafrost table, the quantity of cold energy consumed by importing the cold energy at the 0.5 m interface was less than that at the 1.5 m interface. Therefore, it was more cost-effective to import the cold energy at the 0.5 m interface than at the 1.5 m interface. As a result, the 0.5 m interface was determined as the application position of the interface heat control.
3.3. Determining the Required Cold Energy Quantity for the Interface Heat Control
3.3.1. Initial Determination of the Required Cold Energy Quantity
The quantity of the imported cold energy determines the uplift height of the artificial permafrost table of the underlying permafrost. Insufficient cold energy imported into the subgrade may lead to thawing disease. Conversely, excessive cold energy may lead to heave disease. Preliminarily, the interface heat control target was for the artificial permafrost table to not be higher than the natural ground table within 20 years after the construction of the subgrade. The −4.5 °C case was used as an example to explore the appropriate cold energy quantity to be imported into the subgrade at the 0.5 m interface during the lower thaw season. During the lower thaw season, which covers September–December for a total of 4 months, the cold energy required to offset the heat absorbed by the permafrost subgrade over a year was based on three times the average heat flux of the 0.5 m interface of normal subgrade (3q), and the cold energy imported into the subgrade was gradually increased until the interface heat control target was reached.
Table 2 provides the value of q, which was taken as 0.65 W/m
2.
Figure 7 shows the variation in the artificial permafrost table and the natural ground table after importing different cold energy quantities at the 0.5 m interface during the lower thaw season. As shown in the figure, the artificial permafrost table tended to lift in the 1st–2nd years, the quantity of cold energy was larger, and the lifting tendency was more obvious. This was due to the thermal resistance effect of the subgrade, which prevented the cold energy imported into the permafrost subgrade in the 1st year from being fully transferred to the underlying permafrost. The artificial permafrost table of the imported cold energy subgrade was found to be significantly higher than that of the normal subgrade; as such, the effect of the imported cold energy to uplift the artificial permafrost table was evident. When the cold energy quantity imported into the subgrade reaches 7q, it meets the interface heat control target, thus ensuring that the artificial permafrost table does not exceed the natural ground table within 20 years after the construction of the subgrade.
3.3.2. Optimization of the Required Cold Energy Quantity
Within 20 years after the construction of the subgrade, 7q cold energy from the 0.5 m interface, which was imported during the lower thaw season, can achieve the interface heat control preliminary target. However, during the early part of the 20 years, the 7q cold energy quantity exceeds the cold energy needed to offset the heat absorption of the permafrost subgrade. This results in an artificial permafrost table that is significantly lower than the natural ground table, which may lead to new heave issues. Therefore, it is not appropriate to import the same cold energy quantity during the 20 years, i.e., the cold energy quantity needs to vary within the service life of the subgrade. As shown in
Figure 7, the cold energy quantity required for the interface heat control in the early part of the 20-year period was found to be less than that required in the later part, and the reduction in the cold energy quantity in the 2nd year may have prevented or attenuated the uplift of the artificial permafrost table.
Table 4 shows the artificial permafrost table, the natural ground table, and the difference between them in the 1st year (positive values indicate that the artificial permafrost table was larger than the natural ground table, whereas negative values indicate that the artificial permafrost table was smaller than the natural ground table). As shown in the table, the difference between the natural ground table and the artificial permafrost table was 0 m in the 1st year after importing 4q cold energy; as such, it was determined that the imported cold energy in the 1st year was 4q.
Table 5 shows the artificial permafrost table, the natural ground table, and the difference between them in the 2nd year. As shown in the table, in the 2nd year, when the cold energy was 4q, the artificial permafrost table was higher than the natural ground table; as such, the cold energy must have been reduced based on 4q. When the cold energy was 0, the difference between the artificial permafrost table and the natural ground table was small, i.e., −0.033 m. Therefore, it was determined that the imported cold energy in the 2nd year was 0.
Starting from the 3rd year, the cold energy quantity imported into the subgrade was gradually increased based on 3q. The results of importing different cold energy quantities into the subgrade are shown in
Figure 8, which includes the artificial permafrost table, the natural ground table, and the differences between the two. As shown in the figure, when the cold energy imported into the subgrade was 3q, 4q, and 5q, the difference between the artificial permafrost table and the natural ground table was too large in the 3rd year. When the cold energy quantity was 6q, the maximum difference between the artificial permafrost table and the natural ground table was −0.081 m during the 3rd–10th years. Thus, it was determined that the imported cold energy during the 3rd–10th years was 6q. In the 11th year, the difference between the artificial permafrost table and the natural ground table was 0.111 m; this difference was too large, so the imported cold energy needed to be adjusted.
Starting from the 11th year, the cold energy quantity imported into the subgrade was gradually increased based on 7q. The results of importing different cold energy quantities into the subgrade are shown in
Figure 9, which includes the artificial permafrost table, the natural ground table, and the difference between the two. During the 11th–20th years, when the imported cold energy was 7q, the maximum difference between the artificial permafrost table and the natural ground table was −0.097 m. Therefore, it was determined that the imported cold energy during the 11th–20th years was 7q.
In conclusion, for the −4.5 °C case, the heat control scheme involved importing cold energy into the 0.5 m interface during the lower thaw season, and the cold energy quantities imported into the subgrade were 2.60 W/m3 in the 1st year, no cold energy in the 2nd year, 3.90 W/m3 during the 3rd–10th years, and 4.55 W/m3 during the 11th–20th years.
The two cold energy schemes will be referred to as Scheme 1 (fixed cold energy quantity) and Scheme 2 (varying cold energy quantity), and a comparison will be made between the two schemes. Within the 20 years after the construction of the subgrade, Scheme 2 consumed significantly less cold energy compared to Scheme 1. Specifically, Scheme 1 consumed 9.526 × 104 MJ cold energy, while Scheme 2 consumed 8.301 × 104 MJ cold energy, which is 12.85% less than Scheme 1. Scheme 2 avoided the large uplift of the artificial permafrost table in Scheme 1, and Scheme 2 reduced the maximum difference between the artificial permafrost table and the natural ground table to 0.097 m (compared to 0.697 m achieved in Scheme 1).
Table 6 shows the required cold energy quantity for the interface heat control for the other cases.
4. Analysis of the Cooling Effect of the Interface Heat Control for Permafrost Subgrade
Through interface heat control for the permafrost subgrade, during the 20-year period, the maximum difference between the artificial permafrost table and the natural ground table was 0.097 m. To further investigate the effect of the interface heat control for the permafrost subgrade, the −4.5 °C case was analyzed as an example.
Figure 10 shows the variations in the ground temperature below the road center in the 20-year period. As shown in the figure, in the 20th year, the permafrost under the center of the normal subgrade had completely thawed. Due to the low air temperature, the shallow permafrost had frozen in the cold season, with the maximum freezing depth occurring at the end of May with a long delay. Due to the interface heat control, the underlying permafrost showed a complete freeze–thaw cycle process, i.e., the permafrost had begun to thaw at the end of June and then reached the maximum thawing depth at the end of November. Next, it started to refreeze and had then frozen completely by the beginning of January. Following this, it begun to thaw again at the end of May (meaning a freeze duration of 5 months).
Figure 11 shows the annual average ground temperature variation in the shallow permafrost (3.5 m below the ground surface) and deep permafrost (9.5 m below the ground surface) during the 20-year period. As shown in the figure, under normal subgrade, the shallow permafrost temperature increased by 1.68 °C, i.e., rising from −1.13 °C to 0.55 °C, and the deep permafrost temperature also increased by 0.76 °C, i.e., rising from −1.86 °C to −1.10 °C. In the 8th year, the average annual temperature of the shallow permafrost exceeded 0 °C. Underneath the subgrade, due to interface heat control, the shallow permafrost temperature increased by 0.32 °C, i.e., rising from −1.13 °C to −0.81 °C, and the deep permafrost temperature increased by 0.37 °C, i.e., rising from −1.86 °C to −1.49 °C; moreover, the warming amplitude of the shallow permafrost and the deep permafrost decreased by 80.95% and 51.32%, respectively.
Figure 12a shows the temperature field of the underlying permafrost when the underlying permafrost of the normal subgrade reached the maximum refreezing depth in the 20th year after the construction of the subgrade.
Figure 12b shows the temperature field of the underlying permafrost of the subgrade due to interface heat control at the same moment. As shown in the figure, before the heat control for the permafrost subgrade, the high-temperature thawing nucleus existed in the underlying permafrost when the permafrost reached the maximum refreezing depth. The high-temperature thawing nucleus under the normal subgrade was not conducive to the thermal stability of the permafrost. The strong heat absorption of the permafrost subgrade caused the underlying permafrost to warm up and form a high-temperature thawing nucleus, and the thawing nucleus was driven by the temperature difference to emit heat outward—upward to offset part of the cold during the cold season and downward to transfer heat to the underlying unstable permafrost. However, via the heat control for the subgrade, there was no high-temperature thawing nucleus in the underlying permafrost.
Figure 13 shows the maximum thawing depth of the underlying permafrost at each location with and without interface heat control. As shown in the figure, the maximum thawing depth of the underlying permafrost at each location was significantly reduced by interface heat control, in which the maximum thawing depth under the road center was reduced by the largest magnitude, i.e., 5.10 m. With the interface heat control, the reduction in the maximum thawing depth of the underlying permafrost of the slope and shoulder was found to be smaller compared to that of the road center.
5. Analysis of the Mechanical Ventilation Subgrade Employing Interface Heat Control Parameters
By laying ventilation pipes in the permafrost subgrade and using mechanical ventilation, due to the temperature difference between the subgrade temperature and air temperature, cold air can be used to take away the heat inside the permafrost subgrade to realize the importation of the supposed cold energy. The application position of the interface heat control can determine the location of the ventilation pipes, the application period of the interface heat control can determine the time required to ventilate the subgrade, and the required cold energy quantity for the interface heat control can be determined by calculating the wind velocity at the inlet of the ventilation pipe. Therefore, the ventilation pipe was laid 0.5 m above the ground, and the ventilation time was set from September to December.
Due to the existence of a separation zone in the middle of the wide subgrade, the ventilation pipe can be selected as a T-pipe ventilation pipe with the outlet placed in the separation zone. A T-pipe is more conducive to the stability of the permafrost subgrade than a straight pipe [
40]. The diameter of the ventilation pipe was 0.1 m, and the chosen ventilation pipe was a T-pipe. Based on the Earth–atmosphere-coupled calculation model, a mechanical ventilation calculation model for the permafrost subgrade was established, and the ventilation pipe subgrade is shown in
Figure 14.
The air mass flow rate per second required for heat exchange is as follows:
The relationship between the air mass flow rate and the air velocity in the ventilation pipe is as follows:
The relationship between the air velocity in the ventilation pipe and the imported cold energy is as follows:
where
is the cold energy quantity imported into the subgrade,
is the specific heat of air,
is the temperature difference between the 0.5 m interface and the air temperature during the ventilation time,
is the density of air, and
Apipe is the area of the ventilation pipe inlet.
Based on the predetermined cold energy quantity, the corresponding inlet air speeds for cold energy values of 2.60 W/m3; 0; 3.90 W/m3; and 4.55 W/m3 were 0.42 m/s, 0 m/s, 0.62 m/s, and 0.72 m/s, respectively. Therefore, the inlet air velocities of the ventilation pipe for ventilation pipe subgrades were 0.42 m/s during the 1st year, 0 m/s during the 2nd year, 0.62 m/s during the 3rd–10th years, and 0.72 m/s during the 11th–20th years.
Using the above-calculated inlet air speeds to ventilate the permafrost subgrade,
Figure 15 shows the variations in the underlying permafrost table of the imported cold energy subgrade, mechanical ventilation subgrade, normal subgrade, and natural ground table. As shown in the figure, after the parameters of the mechanical ventilation subgrade determined via the interface heat control parameters were applied, the artificial permafrost table of the mechanical ventilation subgrade was found to be significantly higher than that of the normal subgrade as it fluctuated near the natural ground table during the 20-year period. In addition, the maximum difference between the artificial permafrost table and the natural ground table was found to be 0.15 m. Therefore, using the parameters determined via the interface heat control, the ventilation pipe subgrade also achieved a protective effect for the permafrost. The artificial permafrost table of the mechanical ventilation subgrade was higher than that of the imported cold energy subgrade, which was attributed to a gradual temperature increase in the cold air inside the pipe; as such, it was necessary to amplify the inlet air speed to a certain extent.
6. Discussion
Existing protection measures for permafrost subgrades can be divided into passive insulation measures and active cooling measures. The passive insulation measures comprise awning embankments and crushed rock embankments (during the warm season). The active cooling measures comprise duct-ventilated embankments, two-phase closed thermosyphon (TPCT) embankments, and crushed rock embankments (during the cold season). TPCT embankments remove heat from the permafrost and absorb cold in the environment through the conversion of the liquid nitrogen and nitrogen in TPCT embankments during the cold season [
41]. However, due to the installation location of TPCT embankments, the lateral ground temperature difference in the underlying permafrost was found to be larger, which led to longitudinal cracks in the road [
42]. Awning embankments reduce the impact of direct sunlight on a subgrade via the awning laid on the roadbed slope; however, due to limitations in the awning application location, the underlying permafrost was still affected by the uneven ground temperature, and—with an increase in the service life of the subgrade—the awning can easily be damaged by the wind [
43]. Crushed rock embankments (CREs) and duct-ventilated embankments can ameliorate the problem of uneven ground temperatures that exist in the underlying permafrost. GRE comprise a crushed rock layer with gaps. During the warm season, the crushed rock layer can reduce the heat entering the subgrade; during cold season, the heat in the gaps in the crushed rock layer can be removed under the action of cold air [
43]. The artificial permafrost table of GREs can be uplifted too much, and the protection effect of GREs is too reliant on the thermal properties of the crushed rock layer, as well as on the environment temperature [
44]. The mechanism of the duct-ventilated embankment is to rely on the forced convection heat transfer between cold air and the subgrade to remove the heat in the subgrade [
45]. However, the ventilation pipe can easily be blocked by snow and sand, and the heat removed by cold air from the subgrade is too dependent on the ambient temperature, as well as the wind speed.
Compared with existing protective measures for permafrost subgrades, interface heat control has several advantages. Firstly, existing protective measures for permafrost subgrades mostly focus on the analysis of effects and parameter optimization, but few studies have determined the quantity of cold energy needed for a permafrost subgrade to reach a stable state. This work provides the quantity of the cold energy imported into the subgrade to reach the stable state of the underlying permafrost through the study of the heat transfer process of the permafrost subgrade. Secondly, the protection effectiveness of the permafrost subgrade in existing research relies on the properties and physical dimensions of the materials selected, and there is no way to control the quantity of cold energy imported into the subgrade or the quantity of heat taken away from the subgrade. Interface heat control for permafrost subgrades can quantitatively control the cold energy imported into the subgrade, such that the difference between the artificial permafrost table and the nearby natural permafrost table is small, thereby avoiding the occurrence of the secondary issues caused by excessively high cold energy imported into the roadbed, as well as the issues caused by excessively low cold energy.
In addition, the investigation of interface heat control for permafrost roadbeds can be informative for other permafrost subgrade protection measures. The application location determined via interface heat control can inform the crushed rock layer locations of CREs and the ventilation pipe locations of duct-ventilated embankments, and the application time determined via interface heat control can inform the ventilation time of a duct-ventilated embankment. Furthermore, the quantity of cold energy required, as determined via interface heat control, can be a reference for the material optimization of a CRE, as well as for the calculation of inlet air velocity and the selection of the application temperature for duct-ventilated embankments.
7. Conclusions
In this study, an Earth–atmosphere-coupled calculation model was established by quantitatively analyzing the complex heat transfer process of a 26 m wide permafrost subgrade. The interface heat control parameters and cooling effects of the interface heat control were analyzed, and the thermal performance of a mechanical ventilation permafrost subgrade was calculated. The conclusions were as follows:
- (1)
Interface heat control for permafrost subgrades can ensure the stability of the underlying artificial permafrost. The recommended heat control interface position is 0.5 m, the application period is the lower thaw season, and the imported cold energy should vary during the subgrade service life.
- (2)
By adopting interface heat control measures, the maximum difference between the artificial permafrost table of the permafrost subgrade and the nearby natural ground table is only 0.097 m, and the ground temperature of the underlying permafrost and the area of the thaw bowl decreases significantly. The smoothness of the bottom of the thaw bowl also increases.
- (3)
Due to the cold air inside the ventilation pipe gradually being heated, the inlet air velocity of the ventilation pipe determined via the interface heat control method needs to be partially enlarged in practical applications.
Author Contributions
Conceptualization, Z.L. and F.C.; methodology, Z.L., F.C., H.X. and J.C.; writing—original draft, Z.L., H.X. and B.D.; writing—review and editing, F.C., J.C. and J.L.; investigation, J.C. and Z.L.; resources, J.C., Z.L. and F.C.; data curation, Z.L. and J.C.; formal analysis, J.L., H.X. and J.C.; validation, H.X. and F.C. All authors have read and agreed to the published version of the manuscript.
Funding
This work was supported by the National Science Foundation of China (grant nos. 42371149 and 42230712), the Natural Science Foundation of Shaanxi Province (grant no. 2022JM-143), and the open fund of the National Key Laboratory of Green and Long-Life Road Engineering in Extreme Environment (grant no. YGY2021KFKT03).
Institutional Review Board Statement
Not applicable.
Informed Consent Statement
Not applicable.
Data Availability Statement
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to this project is still in the research phase.
Acknowledgments
We are grateful to the editors and anonymous reviewers for their constructive comments that have improved this manuscript.
Conflicts of Interest
The authors declare no conflicts of interest. Jianbing Chen is employee of National Key Laboratory of Green and Long-Life Road Engineering in Extreme Environment, CCCC First Highway Consultants Co., Ltd. The paper reflects the views of the scientists and not the company.
References
- Phukan, A. Frozen Ground Engineering; John Weily and Sons Inc.: New York, NY, USA, 1985. [Google Scholar]
- Cheng, G.D.; Jin, H.J. Permafrost and groundwater on the Qinghai-Tibet Plateau and in northeast China. Hydrogeol. J. 2013, 21, 5. [Google Scholar] [CrossRef]
- Jin, H.J.; Zhao, L.; Wang, S.; Jin, R. Thermal regimes and degradation modes of permafrost along the Qinghai-Tibet Highway. Sci. China Ser. D Earth Sci. 2006, 49, 1170–1183. [Google Scholar] [CrossRef]
- Yu, Q.H.; Jiang, Z.Q.; Qian, J.; Fan, K.; Gao, L. Analysis of Key Issues in the Construction of Qinghai-Tibet High Grade Highway in China. In Proceedings of the Sixth High-Level Forum on Highway Science and Technology Innovation in China, Beijing, China, 17–18 April 2013. [Google Scholar]
- Yu, Q.H.; Cheng, G.D.; He, N.W.; Pan, X.C. Heat transfer process of permafrost subgrade under different pavement and width conditions. Prog. Nat. Sci. 2006, 16, 5. [Google Scholar]
- Ma, Q.G.; Lai, Y.M.; Wu, D.Y. Research on temperature field of high grade highway subgrade in permafrost region. J. Cent. South Univ. (Nat. Sci. Ed.) 2016, 47, 2415–2423. [Google Scholar]
- Wang, S.J.; Chen, J.B.; Jin, L.; Dong, Y.H.; Chen, D.G. Design theory of highway in permafrost regions based on energy balance. J. Glaciol. Geocryol. 2014, 36, 782–789. [Google Scholar]
- Mu, Y.H.; Ma, W.; Niu, F.J.; Liu, G.; Zhang, Q.L. Study on Geotechnical Hazards to Roadway Engineering in Permafrost Regions. J. Disaster Prev. Mitig. Eng. 2014, 34, 259–267. [Google Scholar]
- Yu, Q.H.; Fan, K.; Qian, J.; Guo, L.; You, Y.H. Key issues of highway construction in permafrost regions in China. Sci. Sin. Technol. 2014, 44, 425–432. [Google Scholar]
- Li, H.D. Simulation analysis of deformation and stability of wide highway roadbed for permafrost region. Road Eng. 2019, 44, 173–177. [Google Scholar]
- Quan, L.; Tian, B.; Niu, K.M.; Li, S.L. Measured Thermal Effect Analysis of Wide Asphalt Pavement-Roadbed System in High-Altitude Permafrost Areas. J. Civ. Eng. 2019, 52, 111–119. [Google Scholar]
- Biskaborn, B.K.; Smith, S.L.; Noetzli, J.; Matthes, H.; Vieira, G.; Streletskiy, D.A.; Schoeneich, P.; Romanovsky, V.E.; Lewkowicz, A.G.; Abramov, A.; et al. Permafrost is warming at a global scale. Nat. Commun. 2019, 10, 264. [Google Scholar] [CrossRef]
- Hjort, J.; Karjalainen, O.; Aalto, J.; Westermann, S.; Romanovsky, V.E.; Nelson, F.E.; Etzelmüller, B.; Luoto, M. Degrading permafrost puts Arctic infrastructure at risk by mid-century. Nat. Commun. 2018, 9, 5147. [Google Scholar] [CrossRef]
- Cheng, G.D.; Wu, Q.B.; Ma, W. Engineering effect of active cooling roadbed of Qinghai-Tibet Railway. Sci. China (E) Eng. Mater. Sci. 2009, 39, 16–22. [Google Scholar]
- Lai, Y.M.; Xu, X.T.; Dong, Y.H.; Li, S.Y. Present situation and prospect of mechanical research on frozen soils in China—ScienceDirect. Cold Reg. Sci. Technol. 2013, 87, 6–18. [Google Scholar] [CrossRef]
- Goering, D.J. Passively Cooled Railway Embankments for Use in Permafrost Areas. J. Cold Reg. Eng. 2003, 17, 119–133. [Google Scholar] [CrossRef]
- Feng, W.J.; Ma, W.; Li, D.; Zhang, L. Application investigation of awning to roadway engineering on the Qinghai-Tibet Plateau. Cold Reg. Sci. Technol. 2006, 45, 51–58. [Google Scholar]
- Bjella, K. Thule Air Base Airfield White Painting and Permafrost Investigation. Phases 1–4; United States Army Corps of Engineers: Washington, DC, USA, 2013. [Google Scholar]
- Shi, L.; Li, N.; Li, G.; Bi, G. Stability Analysis of the Awning in Road Engineering in Permafrost Regions. J. Glaciol. Geocryol. 2007, 29, 986–991. [Google Scholar]
- Pei, W.; Zhang, M.; Yan, Z.; Li, S.; Lai, Y. Numerical evaluation of the cooling performance of a composite L-shaped two-phase closed thermosyphon (LTPCT) technique in permafrost regions. Sol. Energy 2019, 177, 22–31. [Google Scholar] [CrossRef]
- Song, Y.; Jin, L.; Zhang, J. In-situ study on cooling characteristics of two-phase closed thermosyphon embankment of Qinghai–Tibet Highway in permafrost regions. Cold Reg. Sci. Technol. 2013, 93, 12–19. [Google Scholar] [CrossRef]
- Parametthanuwat, T.; Pipatpaiboon, N.; Bhuwakietkumjohn, N.; Sichamnan, S. Heat transfer characteristics of closed-end thermosyphon (CE-TPCT). Eng. Sci. Technol. Int. J. 2021, 27, 101020. [Google Scholar] [CrossRef]
- Wang, S.J.; Long, J.; Hui, P.; Mu, K. Damage analysis of the characteristics and development process of thermosyphon embankment along the Qinghai-Tibet Highway. Cold Reg. Sci. Technol. 2017, 142, 118–131. [Google Scholar] [CrossRef]
- Yu, F.; Zhang, M.; Lai, Y.; Liu, Y.; Qi, J.; Yao, X. Crack formation of a highway embankment installed with two-phase closed thermosyphons in permafrost regions: Field experiment and geothermal modelling. Appl. Therm. Eng. 2017, 115, 670–681. [Google Scholar] [CrossRef]
- Li, D.; Zhang, K.; Tong, G.; Feng, M.; Xing, H. Analysis on Cooling Effect of Crushed-Rocks Embankment of Qinghai-Tibet High-Grade Road. Model. Simul. Eng. 2015, 264, 1–8. [Google Scholar] [CrossRef]
- Liu, M.H.; Niu, F.J.; Fang, J.H.; Lin, Z.J.; Luo, J.; Yin, G.A. In-situ testing study on convection and temperature characteristics of a new crushed-rock slope embankment design in a permafrost region. Sci. Cold Arid Reg. 2014, 6, 110–119. [Google Scholar]
- Wu, Q.B.; Lu, Z.; Zhang, T.; Wei, M.; Liu, Y. Analysis of cooling effect of crushed rock-based embankment of the Qinghai-Xizang Railway. Cold Reg. Sci. Technol. 2008, 53, 271–282. [Google Scholar] [CrossRef]
- Hernandez, M.F.B.; Lemieux, C.; Dore, G. Long-Term Monitoring of Mitigation Techniques of Permafrost Thaw Effects at Tasiujaq Airport. In Proceedings of the 18th International Conference on Cold Regions Engineering and 8th Canadian Permafrost Conference, Nunavik, QC, Canada, 18–22 August 2019. [Google Scholar]
- Wang, Q.Z.; Fang, J.H.; Chao, G. Analysis of cooling effect of block-stone expressway embankment in warm temperature permafrost region. Rock Soil Mech. 2020, 41, 305–314. [Google Scholar]
- Darrow, M.M.; Jensen, D.D. Modeling the performance of an air convection embankment (ACE) with thermal berm over ice-rich permafrost, Lost Chicken Creek, Alaska. Cold Reg. Sci. Technol. 2016, 130, 43–58. [Google Scholar] [CrossRef]
- Malenfant-Lepage, J.; Doré, G.; Fortier, D.; Murchison, P. Thermal Performance of the Permafrost Protection Techniques at Beaver Creek Experimental Road Site. In Proceedings of the International Conference on Permafrost, Yukon, YT, Canada, 25–29 June 2012. [Google Scholar]
- Yuan, K.; Zhang, J.Z.; Fan, K.; Zhu, D.P. Effectiveness and numerical simulation of hollow block ventilated roadbed application in perennial permafrost region. Highw. Transp. Sci. Technol. 2013, 30, 56–62. [Google Scholar]
- Zhang, Z.; Yu, Q.; Wang, J.; Wang, X.; Luo, X. Bidirectional convection mechanism and cooling performance of road embankment with a new duct-ventilated slope in permafrost regions. Cold Reg. Sci. Technol. 2021, 191, 103360. [Google Scholar] [CrossRef]
- Qian, J.; Yu, Q.H.; Wu, Q.B.; You, Y.H.; Guo, L. Analysis of asymmetric temperature fields for the duct-ventilated embankment of highway in permafrost regions. Cold Reg. Sci. Technol. 2016, 132, 1–6. [Google Scholar] [CrossRef]
- Jorgensen, A.S.; Dore, G.; Voyer, E.; Chataigner, Y.; Gosselin, L. Assessment of the effectiveness of two heat removal techniques for permafrost protection. Cold Reg. Sci. Technol. 2008, 53, 179–192. [Google Scholar] [CrossRef]
- Wang, S.J.; Cui, F.Q.; Chen, J.B.; Liu, Z.Y.; Jin, L. Numerical simulations of temperature field for wide subgrade in permafrost regions under earth-atmosphere coupled system. China J. Highw. Transp. 2016, 29, 169–178. [Google Scholar]
- Zhang, M.; Lai, Y.; Liu, Z.; Gao, Z. Nonlinear analysis for the cooling effect of Qinghai-Tibetan railway embankment with different structures in permafrost regions. Cold Reg. Sci. Technol. 2005, 42, 237–249. [Google Scholar] [CrossRef]
- Liu, Z.Y.; Chen, J.B.; Jin, L.; Zhang, Y.J.; Lei, C. Roadbed temperature study based on earth-atmosphere coupled system in permafrost regions of the Qinghai-Tibet Plateau. Cold Reg. Sci. Technol. 2013, 86, 167–176. [Google Scholar] [CrossRef]
- Cheng, G.D. Discussion on the zonation law of high-altitude permafrost in China. J. Geogr. China 1984, 39, 185–193. [Google Scholar]
- Wang, S.W. Optimization of Mechanical Ventilation System for Monolithic Permafrost Subgrade of Qinghai-Tibet Expressway. Master’s Thesis, Chang’an University, Xi’an, China, 2023. [Google Scholar]
- Yu, F.; Qi, J.; Zhang, M.; Lai, Y.; Yao, X.; Liu, Y.; Wu, G. Cooling performance of two-phase closed thermosyphons installed at a highway embankment in permafrost regions. Appl. Therm. Eng. 2016, 98, 220–227. [Google Scholar] [CrossRef]
- Jin, M.; Shang, K.; Yu, Q.; Chen, K.; Guo, L.; You, Y. Study on working performance and cooling effect of a novel horizontal thermosyphon applied to expressway embankment in permafrost regions. Cold Reg. Sci. Technol. 2024, 104147. [Google Scholar] [CrossRef]
- Chen, L.; Voss, C.I.; Fortier, D.; McKenzie, J.M. Surface energy balance of sub-Arctic roads with varying snow regimes and properties in permafrost regions. Permafr. Periglac. Process. 2021, 32, 681–701. [Google Scholar] [CrossRef]
- Wang, J.; Zhao, Y.; Chen, L.; Mu, Y.; Zhang, K.; Li, G. Observation and prediction of long-term thermal condition of permafrost subgrade in high-temperature permafrost region. China Railw. Sci. 2023, 44, 34–45. [Google Scholar]
- Qin, Y.; Zhang, J. A review on the cooling effect of duct-ventilated embankments in China. Cold Reg. Sci. Technol. 2013, 95, 1–10. [Google Scholar] [CrossRef]
Figure 1.
Schematic diagram of the permafrost subgrade physical model.
Figure 1.
Schematic diagram of the permafrost subgrade physical model.
Figure 2.
The thaw curves of the underlying permafrost in the 1st, 5th, 10th, and the 20th years. (a) Tair = −3.0 °C; (b) Tair = −3.5 °C; (c) Tair = −4.5 °C; (d) Tair = −5.5 °C.
Figure 2.
The thaw curves of the underlying permafrost in the 1st, 5th, 10th, and the 20th years. (a) Tair = −3.0 °C; (b) Tair = −3.5 °C; (c) Tair = −4.5 °C; (d) Tair = −5.5 °C.
Figure 3.
Comparison of the artificial permafrost table underneath the subgrade and the natural ground table.
Figure 3.
Comparison of the artificial permafrost table underneath the subgrade and the natural ground table.
Figure 4.
Temperature difference between subgrade interface temperature and air temperature. (a) Tair = −3.0 °C; (b) Tair = −3.5 °C; (c) Tair = −4.5 °C; (d) Tair = −5.5 °C.
Figure 4.
Temperature difference between subgrade interface temperature and air temperature. (a) Tair = −3.0 °C; (b) Tair = −3.5 °C; (c) Tair = −4.5 °C; (d) Tair = −5.5 °C.
Figure 5.
The average heat flux through different interfaces during the lower thaw season.
Figure 5.
The average heat flux through different interfaces during the lower thaw season.
Figure 6.
Variations in the artificial permafrost table after importing different cold energy quantities at the 0.5 m and the 1.5 m interfaces.
Figure 6.
Variations in the artificial permafrost table after importing different cold energy quantities at the 0.5 m and the 1.5 m interfaces.
Figure 7.
Variations in the natural ground table and the artificial permafrost table based on different imported cold energy quantities during the 20-year period.
Figure 7.
Variations in the natural ground table and the artificial permafrost table based on different imported cold energy quantities during the 20-year period.
Figure 8.
Variations in the natural ground table, the artificial permafrost table, and the difference between the two after importing different cold energy quantities during the 3rd–11th years.
Figure 8.
Variations in the natural ground table, the artificial permafrost table, and the difference between the two after importing different cold energy quantities during the 3rd–11th years.
Figure 9.
Variations in the natural ground table, the artificial permafrost table, and the differences between the two after importing different cold energy quantities during the 11th–20th years.
Figure 9.
Variations in the natural ground table, the artificial permafrost table, and the differences between the two after importing different cold energy quantities during the 11th–20th years.
Figure 10.
Variations in the underlying permafrost temperature in the middle of the road in the 20th year. (a) Normal subgrade; (b) Subgrade with interface heat control.
Figure 10.
Variations in the underlying permafrost temperature in the middle of the road in the 20th year. (a) Normal subgrade; (b) Subgrade with interface heat control.
Figure 11.
Variations in the annual average ground temperatures of the underlying permafrost during the 20-year period.
Figure 11.
Variations in the annual average ground temperatures of the underlying permafrost during the 20-year period.
Figure 12.
Temperature field of the underlying permafrost at the moment of the underlying permafrost of the normal subgrade reaching a maximum refreezing depth in the 20th year. (a) Normal subgrade; (b) Subgrade with interface heat control.
Figure 12.
Temperature field of the underlying permafrost at the moment of the underlying permafrost of the normal subgrade reaching a maximum refreezing depth in the 20th year. (a) Normal subgrade; (b) Subgrade with interface heat control.
Figure 13.
Comparison of the maximum thawing depths at different sites in the 20th year.
Figure 13.
Comparison of the maximum thawing depths at different sites in the 20th year.
Figure 14.
Schematic diagram of the ventilation pipe subgrade.
Figure 14.
Schematic diagram of the ventilation pipe subgrade.
Figure 15.
Variation in the permafrost table during the 20-year period.
Figure 15.
Variation in the permafrost table during the 20-year period.
Table 1.
Wind speed at the velocity reference point by month.
Table 1.
Wind speed at the velocity reference point by month.
Month | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
---|
vref/(m·s−1) | 5.69 | 5.83 | 5.92 | 4.79 | 4.37 | 4.11 | 3.72 | 3.37 | 3.38 | 3.58 | 4.29 | 5.23 |
Table 2.
Specifics of four calculation cases.
Table 2.
Specifics of four calculation cases.
Case Indicator | Average Annual Air Temperature | Average Annual Ground Temperature | Active Layer | Permafrost Type |
---|
Case 1 | −3.0 °C | −0.5 °C | 2.67 m | Unstable permafrost |
Case 2 | −3.5 °C | −0.9 °C | 2.36 m | Transitional permafrost |
Case 3 | −4.5 °C | −1.9 °C | 2.03 m | Sub-stabilized permafrost |
Case 4 | −5.5 °C | −2.5 °C | 1.56 m | Sub-stabilized permafrost |
Table 3.
The average annual heat flux at the 0.5 m interface and the 1.5 m interface.
Table 3.
The average annual heat flux at the 0.5 m interface and the 1.5 m interface.
Calculation Condition | Average Annual Heat Flux through the 0.5 m Interface | Average Annual Heat Flux through the 0.5 m Interface |
---|
−3.0 °C | 0.75 W/m2 | 0.76 W/m2 |
−3.5 °C | 0.72 W/m2 | 0.74 W/m2 |
−4.5 °C | 0.65 W/m2 | 0.68 W/m2 |
−5.5 °C | 0.59 W/m2 | 0.61 W/m2 |
Table 4.
The permafrost table and the differences between the natural ground table and artificial permafrost table in the 1st and 2nd years.
Table 4.
The permafrost table and the differences between the natural ground table and artificial permafrost table in the 1st and 2nd years.
Imported Cold Energy Quantity | 3q | 4q | 5q |
---|
Artificial permafrost table | 2.067 m | 2.025 m | 2.011 m |
Difference between artificial table and natural table | 0.042 m | 0 m | 0.014 m |
Table 5.
The permafrost table and differences between the natural ground table and artificial permafrost table in the 2nd year.
Table 5.
The permafrost table and differences between the natural ground table and artificial permafrost table in the 2nd year.
Imported Cold Energy Quantity | 0q | 1q | 2q | 3q | 4q |
---|
Artificial permafrost table | 2.066 m | 1.932 m | 1.865 m | 1.797 m | 1.769 m |
Difference between artificial table and natural table | 0.033 m | 0.101 m | 0.168 m | 0.236 m | 0.264 m |
Table 6.
Required cold energy quantities for the interface heat control for the −3.0 °C case, −3.5 °C case, and −5.5 °C case.
Table 6.
Required cold energy quantities for the interface heat control for the −3.0 °C case, −3.5 °C case, and −5.5 °C case.
Calculation Case | Required Cold Energy Quantity |
---|
−3.0 °C | 3.00 W/m3 for 1st year, 0 for 2nd year, 6.01 W/m3 for 3rd–10th years, 7.50 W/m3 for 11th–20th years |
−3.5 °C | 2.88 W/m3 for 1st year, 0 for 2nd year, 5.76 W/m3 for 3rd–10th years, 6.48 W/m3 for 11th–20th years |
−5.5 °C | 1.77 W/m3 for 1st year, 0 for 2nd year, 2.36 W/m3 for 3rd–10th years, 2.66 W/m3 for 11th–20th years |
| Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).