Prediction Model of Thermal Thawing Sensibility and Thaw Depth for Permafrost Embankment along the Qinghai-Tibet Engineering Corridor Using MODIS Data

The aim of this paper was to reveal the distribution law of permafrost thermal thawing sensibility and thaw depth caused by road construction in Qinghai-Tibet engineering corridor (QTEC). The prediction models of permafrost thermal thawing sensibility and thaw depth have been developed by incorporating the MODIS and in situ soil temperature observation data. The comprehensive earth-atmosphere-coupled numerical models of different embankment structures have been utilized to calculate the thaw depth of the underlying permafrost foundation. Finally, using the given data and above developed prediction models, the distribution maps of permafrost thermal thawing sensibility and thaw depth in QTEC are obtained by grid calculation. The results show the following: (1) Insensitive permafrost of QTEC mainly distributes in the large-scale mountain and high latitude area, and highly sensitive permafrost is located in the perennial river bed, flood plain, and terrace regions. (2) Road construction has a strong thermal disturbance to underlying permafrost, and the proportion of large thaw depth area of separate embankment is obviously smaller than that of 26m full-width embankment. (3) Increase of subgrade interval reduces the proportion of large thaw depth areas, and the application of separate embankment structure is an effective engineering means for the Qinghai-Tibet expressway.


Introduction
With the construction of the Qinghai-Tibet Railway, a series of communication pipelines, oil and gas pipelines, and power grid transmission projects have plan to build in QTEC [1]. The construction of the Qinghai-Tibet expressway has also been put on the agenda, which is included in the "National Highway Network Planning" and "National Highway Network Planning (2013~2030)." Therefore, ensuring the safety of major projects, avoiding engineering instability and reducing the occurrence of geological disasters will be the major issues that need be solved in the engineering construction of the Qinghai-Tibet plateau permafrost regions [2][3][4]. Per-mafrost is a special natural and ecological product of the Qinghai-Tibet Plateau, and its spatial distribution characteristics directly affect the safety and stability of engineering constructions [5][6][7]. Influenced by human activities or engineering structures, the frozen soil at the permafrost table will melt. The different ground temperature, surface characteristics, and ice-containing conditions will cause different thermal respond sensibility to external thermal disturbance and make the difference of thaw depth. The thermal thawing sensibility (Se) can be used as the indicator of frozen soil thermal responds speed [8], and it is also one of the key influencing factors of permafrost characteristics. Thus, the research on the permafrost thermal thawing sensibility and thaw depth distribution in QTEC is not only a prerequisite for environmental monitoring, management, and protection but also a necessary requirement for the engineering structures design and construction [9][10][11][12].
The coupled heat transfer processes between permafrost and the external environment are quite complex and are affected by many factors such as latitude, longitude, elevation, vegetation coverage, soil characteristics, and climate environment [13,14]. So distribution characteristics of thermal thawing sensibility are very regional and special [6,8,15]. Wu et al. [8] proposed that the permafrost thermal thawing sensitivity has a strong correlation with frozen soil temperature and is also related to the seasonal thaw depth (active layer thickness) of frozen soil. It is essential to obtain the distribution of above two influence factors for thermal thawing sensibility determination. Earlier studies on the regionalization of permafrost temperature mainly used artificial drilling and monitoring methods [16]. But the vast area of Qinghai-Tibet Plateau and harsh natural climatic conditions prevented the researchers from obtaining complete, accurate, and timely geothermal observation data, which can usually only be determined the permafrost temperature distribution in a certain area. In recent years, with the development of remote sensing technology and the access to related data, the quantitative studies on the permafrost Z h a M a Q u G a E r Q u (G a R i Q u ) D a n g Q u D a n g Q u Ground temperature measurement point 94°0´0˝E Figure 1: Soil temperature monitoring sections along the QTEC.  [17,18]. Meanwhile, an automatic permafrost temperature monitoring system based on the integrated Global System for Mobile Communications-Railway (GSM-R) and General Packet Radio Service technology (GPRS) has also been applied continuously, thereby realizing permafrost temperature's realtime monitoring and analyzing [16,19,20]. For the permafrost thaw depth research, Nelson and Shiklomanov and Nelson [21,22] used the Stefan method to map the high-precision distribution of active layer in the Kuparuk watershed of northern Alaska. Klene et al. [23] considered the thermal performance and topography of the melt index, soil, water, and surface coverage characteristics of the study area, used the Stefan formula to calculate the active layer thickness of the urbanized area of arctic, and developed a regional probability stochastic model of active layer thickness. Pang et al. [24] considered the influence of vegetation and soil properties on the active layer thickness and obtained active layer thickness distribution map of Qinghai-Tibet Plateau using the Kudryavtsev formula. Wu et al. [25] analyzed the ground temperature data of 10 field observation sites in the Beiluhe region from 2002 to 2012, and the results showed that the average increase of active layer thickness was 4.26 cm/year, and the main reason was the increased rainfall in summer. Chaves et al. [26] monitored the soil temperature of the Keller peninsula in Antarctica, obtained the characteristics of the interannual variability of active layer thickness and near-surface soil heat flow under the climate warming condition, and proved that the thermal regime of patterned     3 Journal of Sensors ground soils was sensitive to air temperatures and atmospheric variations.
With the climate warming, the degradation of the underlying permafrost and the increasing thaw depth caused by engineering construction in QTEC are becoming more serious, and the accompanying engineering risks are also increasing. Therefore, the distribution law of thermal thawing sensibility and thaw depth are the important research topics for permafrost engineering. However, in most of recent studies on the thermal thawing sensibility and thaw depth, the research works are based on the long-term soil temperature monitoring and required a relatively complex mathematical calculation process to obtain the permafrost thermal thawing sensibility [8]. Furthermore, restricted by environmental and economic factors, it cannot arrange a sufficient soil temperature monitoring site along the QTEC, which makes it impossible to obtain the distribution map of permafrost thermal thawing sensibility and thaw depth in the entire corridor. In the present work, to reveal the distribution law of permafrost thermal thawing sensibility and thaw depth caused by road construction in QTEC, the corresponding prediction models have been developed by incorporating the MODIS and in situ soil temperature observation data. The comprehensive earth-atmosphere-   Journal of Sensors coupled numerical models of 26 m full-width and 13 m separate embankment are used to calculate the thaw depth of the underlying permafrost foundation since road construction. Finally, using the given data and above developed prediction models, the distribution maps of permafrost thermal thawing sensibility and thaw depth in QTEC are obtained.
The research results of this paper can contribute to the thermal design and maintenance of engineering structure in permafrost regions.  [27,28]. The daily LST is calculated by the average value of the 11:30 am and 11:30 pm MODIS LST data product. As the cloud, snow, or sensors sometimes induce outliers, NASA adopts the mask processing or outlier identification for the daily LST data during product distribution. The average LST is obtained by the band calculation method using the effective value of daily LST.

Methodology
The vegetation index (VI) reflects the difference between the reflection of vegetation in visible light, near-infrared band, and soil background, which is widely utilized for indicating the status of vegetation growth [19]. In our study, we calculated the annual NDVI values to analyze the trends in vegetation change. The NDVI is expressed on a scale of -1 to +1. The spatial distribution characteristics of surface vegetation and nonvegetation can be determined by reasonable selection of thresholds. In the present study, the average vegetation index, N, is calculated by the arithmetic mean of NDVI which is greater than 0.
Slope is an important index to describe terrain information, and it mainly indicates the inclination of the ground surface. Using the SRTM3-DEM data, surface analysis has been carried out in the ARCGIS 10.2 Software to generate the slope and aspect distribution of the study area. And the equivalent latitude is calculated and applied for determining the sunny or shady of slope direction [29].
where φ is the latitude, k is the slope of location, φ′ is the equivalent latitude, and l is the aspect of location. When the equivalent latitude is less than or equal to the latitude, it is the sunny slope, and vice is shady.

In Situ Observation Data.
To investigate the long-term changes of permafrost along the QTEC, a series of soil temperatures monitoring sections are built along the Qinghai-Tibet railway and highway, which is of 127 mean annual ground temperature (MAGT) monitoring points and 371 active layer depth (ALD) monitoring points. As shown in Figure 1, the monitoring sections are located from Xidatan to the southern foot of Tanggula Mountain, which is also the main permafrost distribution region in QTEC. Soil temperatures are measured at depths of 0.5 m to 30.0 m using thermal probes with a precision of ±0.02°C, and these measurements are made twice a month. The ALD of a certain year is defined as the maximum depth at which the ground temperature is zero in this year. The MAGT is defined as the soil temperature of certain depth basically does not change with time. Thus, the ALD and MAGT of 2016 are calculated from the ground temperature observations.

Prediction
where w 0 is the constant of fitting formula, ðw 1 , w 2 , ⋯, w n Þ are the quantitative indicators of environment, terrain, and other factors related to the permafrost distribution, and ðα 1 , α 2 , ⋯, α n Þ are the coefficients for corresponding multiple regression equation. To achieve better fitting goodness, the prediction accuracy of different combinations of fitting parameters of sunny and shady slope models has been compared. It is found that the influence of LST is significantly greater than equivalent latitude and NDVI for the sunny slope model, and vice versa. Therefore, the optimal parameter combinations of sunny and shady slope models are given in Table 1. Furthermore, statistical test on parameters in the regression model has also been done. It can be seen that the Wals value of every parameters are large and the Sig value is almost less than 0.05, which means each variable in the equation has a high significance and the regression model is considered to pass the test. According to filed survey data of the ratio of permafrost and talik region in QTEC and comparative analysis of prediction accuracy of different identification probability, a fixed value of 0.7 is determined as the identification probability of permafrost: Table 2 is the comparison results of permafrost logistic identification model and measured data. It can be inferred that for the sunny and shady slopes, the corresponding prediction accuracy of permafrost is 93.8% and 92.3%,   Journal of Sensors respectively, and the prediction accuracy of talik is 88.9% and 87.5%, which validates accuracy of probability model.

Multiple Linear Regression Models of MAGT and ALD.
The prediction of MAGT adopts the stepwise regression analysis strategy. The parameters are incorporated into multiple regression models one by one. The selected parameters are tested for significance and partial regression square value. During the fitting process, the significance level value of F is used as the criterion of the stepwise regression method, and the probability of selecting and excluding independent parameters are set to 0.05 and 0.01, respectively. Taking MAGT as the dependent variable and remote sensing data as the independent variables, the stepwise regression analysis is performed to obtain the multiple linear regression model of MAGT (T g ): T g = 56:871-0:952φ-4:736H-0:067φ ′ -5: where φ is the latitude, H is the elevation, φ ′ is the equivalent latitude, and N is the arithmetic mean of vegetation index greater than 0. The statistical test of MAGT prediction model is shown in Table 3. It shows that the Sig values of all fitting parameters are less than 0.05, which satisfy the significance requirements. From the perspective of collinearity, the parameters' tolerance value of prediction model is relatively large, and the VIF value is small, which can be concluded that the collinearity among various fitting parameters is excluded. Considering the ALD is mainly influenced by MAGT, geographical location, surface characteristics, slope, etc., the above multiple regression method is used to obtain the ALD prediction model either. Using a total of 371 in situ observation data along QTEC, the multiple linear regression model of ALD (h) is calculated as follows: where B is the constant of fitting formula, ðx 1 , x 2 , ⋯, x n Þ are the multivariate fitting parameters, and ðβ 1 , β 2 , ⋯, β n Þ are the coefficients for corresponding multiple regression equation. The specific values of the coefficients of regression equation under different MAGT conditions are given in Table 4.         The partial correlation analysis between the permafrost thermal thawing sensitivity and the MAGT and ALD is shown in Table 5. It can be seen from Table 5 that the permafrost thermal thawing sensitivity is significantly correlated with the MAGT and ALD. Through linear regression, the multiple linear regression model of the permafrost thermal thawing sensitivity is obtained as follows: Se = 0:806 + 0:107T g + 0:036h: ð6Þ Figure 2 is the comparison of thermal thawing sensibility calculation results using the prediction model and the empirical formula. It can be seen that the predicted values agrees well with the calculated results of the Kudryavtsev formula using the soil temperature monitoring data [8,15] (R 2 = 0:935), which testifies reasonableness and reliability of the model.

Numerical Calculation Model of Thaw Depth for
Permafrost Embankment. The numerical method is used to calculate the time variation of thaw depth of the underlying permafrost foundation after road construction. Based on the analysis of the complex heat transfer processes among the surrounding environment, embankment, and underlying permafrost, the earth-atmosphere-coupled open numerical model is developed [30,31], and two corresponding physical models for numerical calculation are built: (i) the 26 m width expressway subgrade model with adjacent natural permafrost and air environment and (ii) the separate 13 m width subgrade model with different intervals and comprehensive earth-atmosphere coupled system. The schematic of the above two numerical models are shown in Figure 3.
Considering the distribution pattern of top layer permafrost and external environment conditions in QTEC, 12 computation cases for the above two models have been built to investigate the thaw depth of different soil types, moisture content, and yearly average air temperature [27] (see in Table 6). A 2D unsteady model has been developed in this work to predict the spatial and temporal thermal exchange processes between the earth-atmosphere-coupled system. The standard κ-ε model is used for the air region turbulence simulation. The complex boundary conditions of air regions and source term of ground surface are imported into the FLUENT by the UDF program. The air temperature, wind speed, wind direction, solar radiation intensity, and other data are taken from the meteorological monitoring data, and the annually temperature ascending speed is defined as 0.022°C/year. The thermal thawing sensitivity of the above cases can be calculated from the temperature field computation results of the nonroadbed model in 100th year.

Fitting Formula of Thaw Depth.
Using the developed numerical model, the underlying permafrost thaw depth of the 26 m full-length and separate embankment in 20th year of different thermal thawing sensitivity conditions is calculated (as shown in Figure 4). It can be seen that as the sensitivity of thermal ablation increases, the thaw depth, Δh (amount of the artificial permafrost table to natural permafrost table), increases as well, and the growth trend follows the exponential law. It infers that the thaw depth of permafrost embankment can be predicted with thermal thawing sensitivity.
For the formula of the fitting thaw depth and thermal thawing sensitivity of 26 m full-length and separate embankments, the results are as follows: where Δh w−26m , Δh s−0m , and Δh s−9m are the thawing depth of 26 m full-width embankment, 0 m interval separate embankment, and 9 m interval separate embankment in 20th year, respectively. The determined coefficients (R 2 ) of three fitting formulas are 0.979, 0.838, and 0.886, which exhibit relative better fitness of the calculation formula.

Permafrost Thermal Thawing Sensibility Distribution.
Importing the LST, NDVI, and SRTM3-DEM data into ARCGIS 10.2 Software, the grid calculation has been done based on the thermal thawing sensibility prediction model (Eqs. (2)-(6)), and the permafrost thermal thawing sensibility distribution map can be obtained (as shown in Figure 5). The criteria for classification the thermal ablation sensitivity are shown in Table 7 [8,15].
The results show that the present insensitive permafrost mainly distributes in the large-scale mountain and high latitude area, such as the Kunlun mountain, Fenghuo mountain, and Tanggula mountain. The highly sensitive permafrost is located in the perennial river bed, flood plain, and terrace regions, such as the Chumaer river, Tuotuohe river, and Tongtian river. Meanwhile, as shown in Table 8, sensitive and extremely sensitive permafrost account for 38% for the entire QTEC region based on the current data, and weakly sensitive and insensitive account for 13.8% and 22.2%, 9 Journal of Sensors respectively. It can be expected that permafrost in the QTEC will become more fragile, and the proportion of sensitivetype permafrost will also upsurge with global warming and intensified human engineering activities.

Thaw Depth Distribution of Wide and Separate
Embankment. Figure 6 is the interannual variation of thaw depth of 26 m full-width embankment with different thermal thawing sensibility within 20 years. It can be seen that the thaw depth of the underlying permafrost foundation has a strong correlation with its thermal thawing sensitivity. With the increase of the thermal thawing sensitivity, the permafrost is more sensitive to the external thermal disturbance. Under the same thermal disturbance (3 m height embankment), permafrost thaw depth increases with thermal thawing sensitivity and time. It also should be noticed that speed of thawing becomes slower in the later period (>10 years), indicating that the external thermal disturbance has gradually penetrated into the interior of permafrost foundation, and the whole heat-exchange process gradually tends to a new dynamic equilibrium.
In order to further clear the thermal impact of separated roadbed and full width embankment on underlying permafrost in QTEC, the thaw depth distribution maps of different embankment structures have been calculated with the above proposed thaw depth prediction model (see in Figure 7). The figure shows that road construction has a strong thermal disturbance to the underlying permafrost. Similar to thermal thawing sensitivity distribution, the relatively small thaw depth areas are mainly distributed in the large-scale mountain and high latitude area, such as the Kunlun mountain, Fenghuo mountain, and Tanggula mountain. Their thaw depth is roughly less than 3 m, and some are even below 1 m. The large thaw depth area is mainly distributed in the perennial river bed, floodplain, and river terrace area of the Chumal River, Tuotuo River, or Tongtian River. The thaw depth is generally between 5 and 10 m. Furthermore, it also can be seen that the large thaw depth area of separate embankment is obviously smaller than that of 26 m full-width embankment, especially in the Chumar River, Beilu River, and Buqu River Valley regions. Therefore, it can be concluded that the heat disturbance of separate embankment to the underlying permafrost is weaker than that of the full-width embankment. Table 9 shows the proportion distribution of thaw depth for different embankment structures in QTEC. It can be seen that the proportion of the talik zone in QTEC is basically unchanged (about 31%) for all three kinds of adopted embankment structure. However, the proportions of thaw depth between 5~7 m and >7 m of separate embankment are much less than that of 26 m full-width embankment, and the range of 0~2 m, vise versa, especially for the separate embankment with 9 m interval (best spacing). For the 26 m full-width embankment, the area ratio of large thaw depth (>5 m) is 22.9%, and the small thaw depth (<3 m) area proportion is 21.3%. Comparatively, large and small area proportion of separate embankment with 9 m interval is only 3.7% and 38.8%, respectively. The increase of subgrade interval also reduces the proportion of large thaw depth areas. In summary, the use of separate embankment structure is an effective engineering means for the Qinghai-Tibet expressway.

Conclusion
The distribution law of permafrost thermal thawing sensibility along QTEC, the degradation of the underlying permafrost, and the increasing thaw depth caused by road construction are the important research topics of road engineering in permafrost regions. In the present work, the prediction models of permafrost thermal thawing sensibility and thaw depth have been developed by incorporating the remote sensing and in situ observation data. The comprehensive earth-atmosphere-coupled numerical models of 26 m full-width and 13 m separate embankment are used to calculate the thaw depth of the underlying permafrost foundation since road construction. Distribution maps of permafrost thermal thawing sensibility and thaw depth in QTEC are obtained by grid calculation in the ARCGIS 10.2 Software. The research results show that (1) The insensitive permafrost mainly distributes in the large-scale mountain and high latitude area, and highly sensitive permafrost is located in the perennial river bed, flood plain, and terrace regions. With climate warming and intensified human engineering activities, the permafrost in the QTEC will become more fragile, and the proportion of sensitive-type permafrost will also increase (2) The road construction has a strong thermal disturbance to the underlying permafrost, and the proportion of a large thaw depth area of separate embankment is obviously smaller than that of 26 m full-width embankment (3) The increase of subgrade interval reduces the proportion of large thaw depth areas. Therefore, the application of separate embankment structure is an effective engineering means for the Qinghai-Tibet expressway.

Data Availability
The raw LST, NDVI, and SRTM3 data used to support the findings of this study were downloaded from the Land Processes Distributed Active Archive Center of NASA and Consortium for Spatial Information of CGIAR, which are freely available at https://lpdaac.usgs.gov/ and http://srtm. csi.cgiar.org/ respectively.

Conflicts of Interest
The authors declare that they have no conflicts of interest. Journal of Sensors Central Universities, CHD (Nos. 300102269207, 300102269303), and the Applied Fundamental Research Project of China Communications Construction Co, Ltd (Nos. 2018-ZJKJ-PTJS03, 2016-ZJKJ-02). We are grateful to the anonymous reviewers for their constructive comments to improve this manuscript.