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

Abstract

Permafrost is a key element of the cryosphere and sensitive to climate change. High-resolution permafrost map is important to environmental assessment, climate modeling, and engineering application. In this study, to estimate high-resolution Xing’an permafrost map (up to 1 km²), 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 temperature 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⁵ km² and 6.81 × 10⁵ km², respectively. Ground-surface characteristics were introduced into the permafrost models to improve the performance of models. The results provided a basic reference for permafrost distribution research at the regional scale.

1. 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 (CO₂, CH₄ and N₂O). 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 permafrost 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². Furthermore, we considered the impact of surficial condition on the permafrost distribution.

2. Materials and Methods

2.1. 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

$$\mathrm{thawing} \mathrm{index}={\displaystyle \sum \left|T_i\right|} \left(T_i>0 °\mathrm{C}\right)$$

(1)

$$\mathrm{freezing} \mathrm{index}={\displaystyle \sum \left|T_i\right|} \left(T_i<0 °\mathrm{C}\right)$$

(2)

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³/cm³) and the lowest in the west part (0.1 cm³/cm³) [44].

2.2. 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

$$SFN=\frac{\sqrt{n_f\left|FDD\right|}}{\sqrt{n_f\left|FDD\right|}+E\sqrt{n_{t}TDD}}$$

(3)

where SFN is surface frozen number. T_i is temperature (°C). FDD and TDD represent temperature freezing index (°C days) and temperature thawing index (°C days), respectively. n_f and n_t 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

$$MAGT=\left\{\begin{array}{l}\frac{1}{p}\left(n_{f}FDD+E{n}_{t}TDD\right)\\ for \left(n_{f}FDD+E{n}_{t}TDD\right)\le 0\\ \frac{1}{p}\left(\frac{1}{E}{n}_{f}FDD+n_{t}TDD\right)\\ for \left(n_{f}FDD+E{n}_{t}TDD\right)\ge 0\end{array}\right.$$

(4)

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 n_fFDD is greater than En_tTDD (MAGT > 0 °C), seasonally frozen ground will exist. If n_fFDD is less than En_tTDD (MAGT ≤ 0 °C), permafrost will exist.

2.3. FDD and TDD

The Moderate Resolution Imaging Spectroradiometers (MODIS) has provided LST measurements at the spatial resolution of 1 km² 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² resolution. Second, for each grid (1 km²), 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²), the surface temperature was determined. Finally, the surface temperatures were inserted into the corresponding gaps in MODIS LST time series.

2.4. n (n_f and n_t) 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. n_t factor mainly depends on the land cover type (LCT) during the thawing seasons. As shown in

Table 1. List of n_t factors related to land cover type group.

Land Cover Type Group	nt Factor	IGBP Land Cover Type
Forests	0.8	1, 2, 3, 4
Shrubs	0.85	6, 7, 8, 9
Grasslands and croplands	0.95	10, 12, 14
Wetlands	0.9	11
Urban	1.2	13
Barren	1.1	16

, n_t factor value used in this study was determined based on LCT. The values of n_t 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². On the other hand, n_f 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 n_f 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. List of E factors related to land cover type.

Land Cover Type Group	E Factor	IGBP Land Cover Type
Forests	0.9	1, 2, 3, 4
Shrubs	0.8	6, 7, 8, 9
Grasslands and croplands	0.75	10, 12, 14
Wetlands	0.55	11
Urban	0.7	13
Barren	0.95	16

). LCT was also obtained from MCD12Q1.

2.5. Modeling of Subpixel Heterogeneity

It is well known that permafrost properties can change considerably within the 1 km² pixel due to the heterogeneous vegetation, slope, and aspect. To simulate this variability, we first calculated SFN within the 500 m² subpixel (SFN_ij) by using Equation (1). Temperature within subpixel was equal to gaps-filling LST within the 1 km² pixel. Unique values of n_t and E were determined based on LCT within subpixel. Then, for per 1 km² pixel, the SFN value (SFN_i) was calculated as follows:

$$SF{N}_i=\frac{1}{4}{\displaystyle \sum _{j=1}^{4}SF{N}_{ij}}$$

(5)

If SFN_i is less than 0.5, this pixel is classified into seasonal frozen ground. If SFN_i 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.

$$MAG{T}_i=\frac{1}{4}{\displaystyle \sum _{j=1}^{4}MAG{T}_{ij}}$$

(6)

where MAGT_i is MAGT within the 1 km² pixel. If MAGT_i is more than 0 °C, this pixel is classified into seasonal frozen ground. If MAGT_i is less than 0 °C, this pixel is classified into permafrost.

2.6. 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 temperature 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²) were used for accuracy evaluation [52].

$${\mathrm{R}}^2=\frac{{\displaystyle \sum _{i=1}^S{\left(P_i-\overline{O}\right)}^2}}{{\displaystyle \sum _{i=1}^S{\left(O_i-\overline{O}\right)}^2}}$$

(7)

$$RPD=\frac{SDP}{{\left({\displaystyle \sum _{i=1}^S{\left(P_i-O_i\right)}^2}-\left\{{\left[{\displaystyle \sum _{i=1}^S\left(P_i-O_i\right)}\right]}^2/S\right\}/\left(S-1\right)\right)}^{1/2}}$$

(8)

in which

$$SDP={\left\{{\displaystyle \sum _{i=1}^S{P}_i^2}-\left[\left({\displaystyle \sum _{i=1}^S{P}_i}\right)/S\right]/\left(S-1\right)\right\}}^{1/2}$$

(9)

where S is the number of samples, P_i is the MODIS LST, surface temperature derived from ERA5 climate reanalysis data, or gaps-filling MODIS LST, O_i is the air temperature from NMIC/CMA, and $\overline{O}$ 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.

$$OA=\frac{s}{n}$$

(10)

$$Ka=\frac{s/n-\left(a_1{b}_1+a_0{b}_0\right)/n^2}{1-\left(a_1{b}_1+a_0{b}_0\right)/n^2}$$

(11)

where n is the total number of pixels, s is the total number of correctly classified pixels, a₀ is the total number of observed pixels with permafrost, a₁ is the total number of simulated result with permafrost, b₀ is the total pixels of seasonally frozen ground in the observed result, b₁ 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.

3. Results

3.1. 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⁵ km² and 6.81 × 10⁵ km², respectively.

Spatial distribution of simulated SFN (≥0.5) in Northeast China is shown in Figure 3. SFN value was the highest in the north part of Da Xing’anling Mountains and the lowest in Xiao Xing’anling Mountains. Spatial distribution of simulated MAGT (<0 °C) in Northeast China is shown in Figure 4. MAGT value was the highest in Xiao Xing’anling Mountains and the lowest in the north part of Da Xing’anling Mountains. In summary, permafrost area with a higher SFN value had a relatively lower MAGT value.

3.2. 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² and RPD, these data sets agreed with daily air temperatures. Furthermore, the values of R² 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⁵ km², 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.

4. 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⁵ km² in this study, which is larger than the result (4.23 × 10⁵ km²) 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⁵ km²) 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. n_t value of unity was used by Obu, et al. [48], while n_t value was also constrained based on LCT. As suggested by Way and Lewkowicz [49], it is more reasonable that n_t 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.

5. 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⁵ km² and 6.81 × 10⁵ km², 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.