Spatial Distribution of Permafrost in the Xing’an Mountains of Northeast China from 2001 to 2018

: Permafrost is a key element of the cryosphere and sensitive to climate change. High ‐ reso ‐ lution permafrost map is important to environmental assessment, climate modeling, and engineer ‐ ing application. In this study, to estimate high ‐ resolution Xing’an permafrost map (up to 1 km 2 ), we employed the surface frost number ( SFN ) model and ground temperature at the top of permafrost (TTOP) model for the 2001–2018 period, driven by remote sensing data sets (land surface tempera ‐ ture and land cover). Based on the comparison of the modeling results, it was found that there was no significant difference between the two models. The performances of the SFN model and TTOP model were evaluated by using a published permafrost map. Based on statistical analysis, both the SFN model and TTOP model efficiently estimated the permafrost distribution in Northeast China. The extent of Xing’an permafrost distribution simulated by the SFN model and TTOP model were 6.88 × 10 5 km 2 and 6.81 × 10 5 km 2 , respectively. Ground ‐ surface characteristics were introduced into the permafrost models to improve the performance of models. The results provided a basic refer ‐ ence for permafrost distribution research at the regional scale.


Introduction
Permafrost accounts for about one quarter of the Northern Hemisphere and thus is an important component of the cryosphere [1]. It is estimated that large amounts of ground ice and organic carbon are stored in permafrost [2][3][4]. During recent decades, permafrost could experience significant degradation due to global warming at the regional and global scales [5,6]. The permafrost degradation could lead to the melting of ground ice and substantial emissions of greenhouse gases (CO2, CH4 and N2O). Such a change will, in turn, have serious impacts on ecosystem [7,8], local hydrological system [9,10], and engineering projects [11,12]. In order to estimate accurately the historical or future change of permafrost, it is necessary to study the accurate contemporary permafrost distribution [13]. On the other hand, the permafrost distribution is mainly controlled by microclimatic factors such as snow, air temperature, moisture, vegetation, and soil texture [14,15]. Due to the obvious spatial difference of microclimatic factors, understanding the permafrost distribution at the regional scale is important for investigating the global permafrost distribution [16].
Xing'an permafrost is mainly located in the southern limit of permafrost and restricted to the Da and Xiao Xing'anling Mountains, Northeast China [17,18]. This perma-frost is the highest latitude and the coldest region in China [19,20]. The continuous permafrost is located in the north part of Da Xing'anling Mountains. Different from Norway or the Qinghai-Tibet Plateau (QTP), most of the area is covered by forest, shrub, and grassland [21,22]. These properties make it a challenging area for permafrost study. Xing'an permafrost is more sensitive to climate change due to its higher temperature (close to 0 °C). Many studies found that air temperature in Northeast China showed a significant increasing trend during recent decades [23][24][25][26][27]. However, few studies analyzed the Xing'an permafrost distribution [28,29]. The present available permafrost maps were overly simple, which were compiled only based on air temperature. Furthermore, these maps were at low resolutions.
Traditional methods for mapping permafrost distribution are field investigations (geophysical survey and borehole) that are extremely costly, time-consuming, and reliable over small-scale areas [16]. During recent decades, the permafrost models are developed to map permafrost distribution at the large scale. The permafrost models are mainly divided into two types: empirical models and physical models. The empirical model is developed by establishing relationships between permafrost probability and many external factors, such as elevation, mean annual ground temperature, and vegetation [30][31][32]. However, the empirical model neglects the impact of soil factors, such as soil moisture and soil texture, on permafrost distribution. These factors take important roles in the hydrothermal process within permafrost [33]. As a result, the empirical model may lead to biases in mapping permafrost distribution. Different from the empirical model, the soil factors can be expressed by parameters in the physical model [34,35]. Furthermore, the physical model can describe more details on the hydrothermal exchange processes between ground and the atmosphere [36]. Therefore, the physical model is theoretically superior to the empirical model. However, many important parameters in the physical models are difficult to obtain at the regional scale, which may make it difficult for the physical models to accurately determine the spatial distribution of permafrost.
The remote sensing technique has provided land surface data sets that can be used to characterize the surfical conditions [37][38][39]. These remote sensing data sets can be applied to determine the spatiotemporal distribution of model parameters at the regional scale [40]. Therefore, integration of remote sensing data sets and the physical models provides an effective way to simulate the permafrost distribution and estimate the permafrost extent [15,41]. Furthermore, the higher spatiotemporal resolution remote sensing data sets can improve the resolution of permafrost map and accuracy of permafrost extent [42,43].
In this study, we used two physical models (surface frost number (SFN) model and ground temperature at the top of permafrost (TTOP) model) to simulate the contemporary permafrost distribution in Northeast China. Remote sensing data sets were used as inputs to physical models. Outputs were compared with a published permafrost map. Different the previous studies, we compiled the Xing'an permafrost map with a resolution of 1 km 2 . Furthermore, we considered the impact of surficial condition on the permafrost distribution.

Study Area
Xing'an permafrost is distributed in Da and Xiao Xing'anling Mountains in the north part of Northeast China. It is characterized by lower elevation in the middle part and higher in the west and east part ( Figure 1). The surface geology of study area consists of platforms, hills, and mountains. Most of the area is covered by moss layer, grasslands, shrubs, and forests. In aspect of climate condition, the mean annual temperature was close to 0 °C. The thawing and freezing index are important for permafrost study. They are expressed as From the lower to the higher elevations, the thawing index ranges from 1887 °C•d to 3192 °C•d, and the freezing index ranged from −1633 °C•d to −3464 °C•d [29]. The precipitation amount gradually decreased from approximate 550 mm in the east to 300 mm in the west. The magnitude of soil moisture was the greatest in the east part (0.4 cm 3 /cm 3 ) and the lowest in the west part (0.1 cm 3 /cm 3 ) [44].

The SFN Model and TTOP Model
The SFN and TTOP models are two simple physical models. These models employ only a small number of input parameters and are successfully used to estimate the permafrost distribution at the regional or global scales. The SFN model determines the permafrost distribution by the comparison between the thaw depth and the freeze depth [45]. The SFN model is expressed as where SFN is surface frozen number. Ti is temperature (°C). FDD and TDD represent temperature freezing index (°C days) and temperature thawing index (°C days), respectively. nf and nt represent n factor of the freezing and thawing seasons. E is the ratio of thermal conductivity of the active layer when soil is thawing and freezing. The SFN value is used to distinguish between permafrost and seasonally frozen ground. The SFN value of ≥0.5 is taken as permafrost. The SFN value of <0.5 is taken as seasonally frozen ground. The TTOP model was developed to determine the occurrence of permafrost by using the MAGT at the top of permafrost or the base of seasonally frozen layer [35,46]. The TTOP model is expressed as where MAGT is the mean annual ground temperature at the top of permafrost or the base of seasonally frozen layer. p is the annual period (365 days). The other parameters in the TTOP model are consistent with the SFN model. If nfFDD is greater than EntTDD (MAGT > 0 °C), seasonally frozen ground will exist. If nfFDD is less than EntTDD (MAGT ≤ 0 °C), permafrost will exist.

FDD and TDD
The Moderate Resolution Imaging Spectroradiometers (MODIS) has provided LST measurements at the spatial resolution of 1 km 2 since 2000. Recently, the MODIS data have been widely used in different models [15,33,40]. Therefore, in this study, we used the MODIS MOD11A1 level 3 product ranged from 2001 to 2018 to calculate FDD and TDD. The MODIS MOD11A1 level 3 product contains up to daytime and nighttime measurements per day. The daily LST was estimated by averaging the two observations.
There are numerous data gaps in the MOD11A1 time series duo to cloud cover or other factors, which could cause calculation errors in FDD and TDD. To prevent this, gaps were filled by surface temperature derived from ERA5 climate reanalysis data (ERA5 hourly data on pressure levels). The ERA5 climate reanalysis data were obtained from climate data store website (//cds.climate.copernicus.eu/, accessed on 20 June 2020). These data provide temperature at different elevation with the temporal resolution of 1 h and spatial resolution of 0.25° × 0.25°. The EAR5 climate reanalysis data below 2000 m were used in this study.
The method of filling gaps of MODIS LST data in this study was consistent with the method proposed by Fiddes and Gruber [47] and Obu, et al. [48]. First, ERA5 climate reanalysis data were accumulated to the 1 day resolution and then downscaled to the 1 km 2 resolution. Second, for each grid (1 km 2 ), the linear interpolation method was used to estimate temperature profile with vertical resolution of 1 m. According to the elevation provided by GMTED2010 data (about 1 km 2 ), the surface temperature was determined. Finally, the surface temperatures were inserted into the corresponding gaps in MODIS LST time series.

n (nf and nt) and E Factors
n factors not only quantify the relationship between near surface temperature and ground surface temperature but also characterize the surface energy balance including net wave radiation, latent heat, sensible heat, and ground surface heat flux. It is an effective way to transform temperature near the ground surface to the ground surface temperature. nt factor mainly depends on the land cover type (LCT) during the thawing seasons. As shown in Table 1, nt factor value used in this study was determined based on LCT. The values of nt factor defined by LCTs were consistent with the study proposed by Way and Lewkowicz [49]. LCT was obtained from remotely sensed gridded products (MODIS MCD12Q1), which ranged from 2001 to 2018. The MODIS LCT product includes 17 types with the spatial resolution of 500 m 2 . On the other hand, nf factor mostly depends on snow depth. Snow depth was generally less than 30 cm in Northeast China [40], which cannot protect the ground surface temperature from the impact of the shortwave radiation [50].
Thus, the value of nf factor was set to unity in this study. The parameter E is mainly associated with soil moisture. Wang, et al. [51] found that soil moisture mainly depends on LCT based on continuous monitor data. Thus, LCT was a suitable proxy for soil moist. Following the study of Obu, et al. [48], E value was also determined according to LCT in this study (Table 2). LCT was also obtained from MCD12Q1.

Modeling of Subpixel Heterogeneity
It is well known that permafrost properties can change considerably within the 1 km 2 pixel due to the heterogeneous vegetation, slope, and aspect. To simulate this variability, we first calculated SFN within the 500 m 2 subpixel (SFNij) by using Equation (1). Temperature within subpixel was equal to gaps-filling LST within the 1 km 2 pixel. Unique values of nt and E were determined based on LCT within subpixel. Then, for per 1 km 2 pixel, the SFN value (SFNi) was calculated as follows: If SFNi is less than 0.5, this pixel is classified into seasonal frozen ground. If SFNi is more than 0.5, this pixel is classified into permafrost. It should be noted that Xing'an permafrost is a binary classification in this study.
For the TTOP model, the procedure of determining permafrost distribution was consistent with the SFN model procedure.
where MAGTi is MAGT within the 1 km 2 pixel. If MAGTi is more than 0 °C, this pixel is classified into seasonal frozen ground. If MAGTi is less than 0 °C, this pixel is classified into permafrost.

Accuracy Evaluation Method
The filling gaps method is important for calculating the model inputs (FDD and TDD). To evaluate the performance of filling gaps method, we employed daily air tem-perature from 47 meteorological stations. Figure 1 shows the distribution of meteorological stations in Northeast China. This dataset is provided by National Meteorological Information Center, China Meteorological Administration (NMIC/CMA). We compared air temperatures with the raw MODIS LSTs, surface temperatures derived from ERA5 climate reanalysis data, and gaps-filling MODIS LSTs, respectively. Ratio of prediction to deviation (RPD) and determination coefficient (R 2 ) were used for accuracy evaluation [52].
where S is the number of samples, Pi is the MODIS LST, surface temperature derived from ERA5 climate reanalysis data, or gaps-filling MODIS LST, Oi is the air temperature from NMIC/CMA, and is the averaged air temperature. In order to evaluate the accuracy of simulated results, we used a published permafrost map [53], the 1:10,000,000 map of permafrost in China (Map08), as verification data. The corresponding study region for Map08 was extracted to compare our results. The Map08 provided by National Tibetan Plateau Data Center was compiled by field survey. The Map08 was widely used to explore the permafrost distribution in China [20,29,30,54]. Thus, Map08 was worked as the true in this paper. The overall accuracy (OA) and kappa coefficient (Ka) were used to evaluate the agreement of simulated results and Map08 [55]. OA is the proportion of the number of correctly classified pixels to the total number of sample pixel. Ka is used to measure the agreement between the simulated map and the observed map.
where n is the total number of pixels, s is the total number of correctly classified pixels, a0 is the total number of observed pixels with permafrost, a1 is the total number of simulated result with permafrost, b0 is the total pixels of seasonally frozen ground in the observed result, b1 is the total pixels of seasonally frozen ground in the simulated result. Shi, et al. [55] suggested that the value of Ka as follow: Ka > 0.8 indicates excellent agreement, 0.6 < Ka < 0.8 indicates substantial agreement, 0.4 < Ka < 0.6 represents moderate agreement, 0.2 < Ka < 0.4 represents fair agreement, and Ka < 0.2 corresponds to lack of agreement.

Permafrost Distribution
Based on the SFN model and TTOP model binary classification, the contemporary Xing'an permafrost distribution in Northeast China is shown in Figure 2. Xing'an permafrost was mainly distributed in Da and Xiao Xing'anling Mountains in which the southern limit of permafrost occurred. The permafrost distribution obtained by the SFN model and TTOP model agreed excellently with each other with Ka of 0.989 and OA of 99.5%. The differences were not significant (data not shown). The extent of the permafrost region inferred from the SFN model and TTOP model were 6.88 × 10 5 km 2 and 6.81 × 10 5 km 2 , respectively.

Filling Gaps Method and Models Validation
First, we evaluated the performance of filling gaps method. Figure 5 shows the comparison between air temperature from 47 meteorological stations and the raw MODIS LST, surface temperature derived from ERA5 climate reanalysis data, and gap-filling MODIS LST, respectively. Based on the values of R 2 and RPD, these data sets agreed with daily air temperatures. Furthermore, the values of R 2 and RPD were close to each other, respectively. These results indicated that the gap of MODIS LST can be filled by the surface temperature derived from ERA5 reanalysis climate data. Second, we evaluated the accuracy of simulated Xing'an permafrost maps. Figure 6 shows the comparison between Map08 and the simulated results. Map08 estimated the extent of permafrost region in Northeast China to be 5.74 × 10 5 km 2 , which are close to our results. The Ka values of the SFN model and TTOP model were 0.765 and 0.773. The OA values of the SFN model and TTOP model were 0.88 and 0.885. All these indices indicated that the SFN model and TTOP model reasonably capture the permafrost distribution in Northeast China and agreed substantially with the published permafrost maps. As shown in Figure 6, the differences mainly appear near the boundary of permafrost and seasonally frozen ground. The Xing'an permafrost distribution was overestimated. The SFN and TTOP models overestimated permafrost area in the east of 130° E and south boundary of Da Xing'anling Mountains permafrost region.

Discussion
The comparison result showed that the map of permafrost distribution obtained by using the SFN model agrees excellently with that obtained by the TTOP model. This may be the reason that both models take into account the change of soil physical properties caused by phase transition. In addition, if value of parameter E is set to unity, the map of permafrost distribution by using the SFN model is consistent with that by using the TTOP model based on Equations (1) and (2).
The extent of Xing'an permafrost was about 6.8 × 10 5 km 2 in this study, which is larger than the result (4.23 × 10 5 km 2 ) from Zhang, et al. [29] using the SFN model. Although the SFN model was used in the two studies, many differences still caused the inconsistency of area. First, E value of unity was applied by Zhang, et al. [29], while E value was determined based on LCT in this study. Obu, et al. [48] suggested that E value determined by LCT is more reasonable. Second, different spatial resolution of forcing data may lead to discrepancies. Recently, Obu, et al. [48] used the TTOP model to map the North Hemisphere permafrost. In order to compare our results, we extracted the corresponding study region. The simulated result from Obu, et al. [48] (5.64 × 10 5 km 2 ) is lower than our result. Although the simulated result from Obu, et al. [48] is closer to Map08, the values of OA and Ka (0.844 and 0.687) are lower than our results. The difference may be due to the different parameterization scheme. nt value of unity was used by Obu, et al. [48], while nt value was also constrained based on LCT. As suggested by Way and Lewkowicz [49], it is more reasonable that nt value varies according to LCT for using the TTOP model.
The SFN and MAGT values were probably related to the effect of elevation and latitude. Both elevation and latitude had cooling effect on temperature. Temperature decreased significantly with elevation, leading to increase of FDD with elevation and decrease of TDD with elevation. Latitude had the same impact on FDD and TDD. FDD increased with increasing latitude, and TDD increased with decreasing latitude. According to Equations (1) and (2), elevation and latitude are positively correlated with SFN and negatively correlated with MAGT, respectively. As shown in Figure 1, elevation in the south part of Da Xing'anling Mountains was similar with the north part of Da Xing'anling Mountains. Therefore, latitude was the dominant factor to cause higher SFN and lower MAGT in the north part of Da Xing'anling Mountains. On the hand, elevation in Da Xing'anling Mountains was higher than that in Xiao Xing'anling Mountains ( Figure 1). Therefore, elevation was the main reason that caused higher SFN and lower MAGT in Da Xing'anling Mountains than those in Xiao Xing'anling Mountains.
Although the SFN and TTOP models used in this study perform well, there are many uncertainties in simulated results. First, the SFN model and TTOP model are equilibrium models, which neither consider the effect of warming climate on the ground temperature nor the history of ground thermal conditions. The air temperature in Northeast China experienced a significant warming during recent decades [56][57][58][59][60][61], which lead to the disequilibrium between surface climate and ground temperature. This may produce some biases when applying the SFN model and TTOP model to obtain the permafrost distribution in Northeast China. In further studies, some numerical models using more complex parameters should be used to further accurately simulate and analyze the dynamic change of permafrost distribution caused by climate change in the Xing'an Mountains. Second, the limited availability of soil surveys data sets may lead to another problem when simulating the permafrost distribution. Third, the cut values for the SFN model and TTOP model are critical parameters for distinguish permafrost and seasonally frozen ground. However, those cut values are only theoretical and need relatively sufficient knowledge. Further investigation for cut value is required. Reasonable cut value could possibly improve the accuracy of models.

Conclusions
In this study, two physical models (SFN and TTOP) were used to explore the extent of permafrost region in Northeast China. According to the values of OA and Ka, the SFN model agrees excellently with the TTOP model. The performances of the two models were compared with a published permafrost map. Model validations showed that the SFN model and TTOP model agreed substantially with the published permafrost maps, respectively. The OA values were 0.88% and 0.885% for the SFN model and the TTOP model, respectively. The Ka values were 0.765 and 0.773 for the SFN model and TTOP model, respectively. It was concluded that the SFN model and TTOP model can be used to simulate the Xing'an permafrost distribution in Northeast China. The extent of permafrost region simulated by the SFN model and TTOP model are 6.88 × 10 5 km 2 and 6.81 × 10 5 km 2 , respectively, in Northeast China from 2001 to 2018. Based on the distribution of the SFN and TTOP values, the value of SFN and TTOP was associated with elevation and latitude. The new Xing'an permafrost map is necessary for future permafrost research.

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 first author.

Conflicts of Interest:
The authors declare no conflict of interest.