The Mongolian Plateau is located in the interior of Asia. It covers approximately 2.72 × 10⁶ km², including Mongolia, southern Russia, the northern part of the Inner Mongolia Autonomous Region in China, and parts of the Xinjiang Uygur Autonomous Region (Li et al., 2020) (Figure 1). The elevation of the Mongolian Plateau gradually decreases from northwest to southeast, with an average elevation of 1,580 m above sea level (a.s.l.). The western and northwestern regions are mainly mountainous with elevations between 2000 and 4,000 m; the central and eastern regions are mainly plains and hills, with elevations of approximately 800–1,500 m; and the southern and southeastern regions are the Gobi Desert, with elevations of approximately 1,000 m (Chen et al., 2020).

The Mongolian Plateau is located in the center of the Asian high-pressure system, which is one of the origins of cold waves in the winter monsoon climate zone (Chen et al., 2020). The climate type is dominated by arid and semiarid climate characteristics with large intra-annual variability in temperature and precipitation (Zhou et al., 2012). Due to the influence of water vapor from the Arctic Ocean to the north and the Pacific Ocean to the east, precipitation on the Mongolian Plateau decreases spatially from southeast to northwest. The mean annual precipitation is approximately 267 mm, and the Selenge River basin in the north can reach an annual precipitation of 500 mm. The southwestern region has an annual precipitation of only approximately 100 mm (Wang et al., 2008). In addition, evapotranspiration in Mongolia showed a gradual increasing trend from 2001 to 2015. The maximum evapotranspiration values were found in the continuous and discontinuous permafrost regions (Yu, 2017).

The permafrost distribution data in the Northern Hemisphere (Figure 1A) are derived from Obu et al. (2019) via the data source https://apgc.awi.de/dataset/pex.

The monthly reanalysis data of air temperature were used to analyze the spatiotemporal variations in the freezing and thawing index on the Mongolian Plateau. These data were obtained from the CRU time-series Version 4.05 (CRU TS v4.05) datasets (https://crudata.uea.ac.uk/cru/data/hrg/), and the time series is from 1901 to 2019 with a spatial resolution of 0.5 × 0.5. The CRU dataset combines centennial-scale climate change grid data covering the global scale. The dataset is integrated from several well-known datasets through interpolation and other methods, which are widely used as credible proxy data after comparison (Wen et al., 2006; Chen et al., 2021). In addition, the CRU dataset is continuous in space and time and can reflect long-term trends (Shi et al., 2019). The observational records in most of the meteorological sites are from the post-1980s period, thus limiting the study of early climate and freezing and thawing index changes. The long time series provided by these data can effectively support the research in this study.

The applicability of the CRU data on the Mongolian Plateau was examined and evaluated using observations from 106 meteorological stations in the study area obtained from the National Oceanic and Atmospheric Administration (NOAA https://www.climate.gov/maps-data). Since the time series recorded at the selected meteorological stations were not consistent, we selected temperature data from January 1975 to December 2019 for the quality check of the CRU data.

For the accuracy assessment of the CRU reanalysis dataset, four indicators, namely, the mean deviation error (MBE), correlation coefficient (R), normalized standard error (NSSE), and root mean square error (RMSE) are used in this study and are calculated as follows (Qin et al., 2021): M B E = ∑ n = 1 M ( C R U n − O B n ) M (1) R M S E = ∑ n = 1 M ( C R U n − O B n ) 2 M (2) N S S E = ∑ n = 1 M ( C R U n − O B n ) 2 ∑ n = 1 M ( O B n ) 2 (3) R M S E = ∑ n = 1 M ( C R U n − C R U ¯ ) ( O B n − O B ¯ ) ∑ n = 1 M ( C R U n − C R U ¯ ) 2 ∑ n = 1 M ( O B n − O B ¯ ) 2 (4) Where M is the total number of samples and C R U n and O B n (n = 1, 2, 3, M) represent the monthly air temperature corresponding to the CRU and meteorological stations. C R U ¯ and O B ¯ represent the monthly mean values corresponding to the years of the CRU and meteorological stations, respectively.

The air freezing (thawing) index is the cumulative sum of the air temperatures corresponding to the number of days when the mean air temperature is equal to or below (above) 0°C during the freezing (thawing) period (Eqn. 5). The freezing period is usually defined as July 1 of each year to June 30 of the following year, while the thawing period is defined as January 1 to December 31 of each year (Wu et al., 2011; Luo et al., 2014; Shi et al., 2019). T D D = ∫ t 0 t 1 | T i | d t , T i < 0°C F D D = ∫ t 2 t 3 T j d t , T j > 0°C (5) T D D = ∑ i = 1 M T T j , T j > 0°C F D D = ∑ i = 1 M F | T i | , T i < 0°C (6) Where T i and T j refer to the freezing and thawing periods when the air temperature is greater than 0°C and less than 0°C, respectively. t 0 and t 1 represent the beginning and end times of the freezing period, and t 2 and t 3 represent the beginning and end times of the thawing period, respectively. Eqn. 5 can be simplified to Eqn. 6.

However, the acquisition of daily air temperatures on the Mongolian Plateau during 1901–2019 is limited by many factors, and thus the CRU TS V4.05 monthly dataset was used to determine the freezing and thawing index in this study, so the equation can be further written as Eqn. 7 (Shi et al., 2019). T D D = ∑ j = 1 M T T j ¯ · D j , T j ¯ > 0°C F D D = ∑ i = 1 M F T i ¯ · D i , T i ¯ < 0°C (7) Where T i ¯ and T j ¯ represent the monthly mean values of the freezing and thawing periods, respectively. D i and D j represent the number of days in the months corresponding to the freezing and thawing periods, respectively.

The M–K method (Mann, 1945; Kendall, 1948) is a nonparametric statistical test that has been widely used in trend analysis of time series (You et al., 2010; Wu et al., 2011). In addition, the algorithm can obtain the year in which the trend changes when the value of intersection points is within the confidence intervals (Shi et al., 2019). In the M–K trend analysis, the statistic S of the test is calculated based on the relationship between the time series n and the corresponding continuous values x j and x i , and is defined as follows: S = ∑ i = 1 n − 1 ∑ j = i + 1 n sgn ( x j − x i ) sgn ( x j − x i ) = { 1 , x j − x i > 0 0 , x j − x i = 0 − 1 , x j − x i < 0 (8)

When sample number n > 10, the standard normal system variables are given by Eqn. 9, and the variance in the sample, V ar ( S ) is calculated by Eqn. 10: Z = { S − 1 V ar ( S ) S > 0 0 S = 0 S + 1 V ar ( S ) S < 0 (9) V ar ( S ) = n ( n − 1 ) ( 2 n + 5 ) 18 (10)

In the bilateral trend test, the trend of the time series data can be judged by the statistic Z at a given confidence level: it is considered an upward trend when Z >0 and a downward trend when Z <0. The magnitude of the Z value is used to determine whether the trend is significant.

In the M–K mutation test, the statistic is defined under the assumption that the time series are independent as follows: S k = ∑ i = 1 k R i ( k = 2,3 , ... , n ) U F k = S k − E ( S k ) V a r ( S k ) ( k = 1,2 , ... , n ) U B k = − U F k k = n , n − 1 , ... , 1 (11) Where R i denotes the cumulative number of samples in which x i > x j (1 ≤ j ≤ i). When k = 1, U F 1 = 0, and U B 1 = 0, E ( S k ) and V a r ( S k ) represent the mean and variance, respectively, which can be calculated as follows: E ( S k ) = n ( n + 1 ) / 4 V a r ( S k ) = n ( n − 1 ) ( 2 n + 5 ) / 72 (12)

If the two curves U F k and U B k intersect between the critical lines, the change year is considered to correspond to the intersection point.

Figure 2 shows the annual and monthly fitting results of the meteorological station data and CRU reanalysis data from 1901 to 2019. According to the relevant impacts of different months (Figure 2B), the CRU reanalysis data, compared with meteorological observation data, underestimate the actual air temperature on the Mongolian Plateau, especially in winter (December–February) and spring (March–May). In contrast, the results for both summer and autumn are similar, except for a slight underestimation in July and September. The annual fitting result from 1901 to 2019 (Figure 2A) shows a significant positive correlation, with a Pearson’s correlation coefficient of 0.96 and an explained variance of 0.92, indicating that the CRU dataset is credible for studying air temperature correlation on the Mongolian Plateau.

The indicators for evaluating the applicability of monthly CRU reanalysis data on the Mongolian Plateau were spatially visualized via ArcMap spatial analysis (Figure 3). The monthly CRU data underestimated the air temperature for 70% of stations, with the mean MBE value being −0.77 and the variation ranging between −7.38 and 2.47. The variation in MBE on the Mongolian Plateau was mostly concentrated between −1 and 0°C. The spatial distribution of the correlation between the CRU data and the meteorological observations at site locations is greater than 0.95 for most sites, with a range of 0.90–0.99. The mean value of NSSE is 0.68, and the value is less than 0.6 for approximately 68% of the sites. The range of RMSE was 0.13–6.78°C, with a mean value of 0.80°C, and approximately 82% of the sites had RMSE values less than 2°C. Based on the statistical results and spatial variations in the four evaluation indices, the R values are high, and the MBE, RMSE, and NSSE values are low. In terms of spatial distribution, most stations show a positive relationship, but a few stations are in the western and northern border areas. Therefore, the CRU reanalysis dataset has good applicability to the Mongolian Plateau region and can be used as a proxy for a long time series for the Mongolian Plateau region.

Based on the multivariate statistical relationship (Eqn. 8) between the grid data of the mean freezing index calculated from 1901 to 2019 and the longitude, latitude, and elevation of the corresponding points, it can be seen that the spatial distribution of the freezing index on the Mongolian Plateau is influenced by latitude and altitude. As shown in Figure 4, the spatial distribution of the freezing index on the Mongolian Plateau exhibits a significant latitudinal gradient. It increases with latitude, with an average value ranging from 1,200 to 3,600°C·d. On the Mongolian Plateau, the highest freezing index value occurs in northern regions, with a value of nearly 5,600°C·d, whereas in the southern Gobi Altai and surrounding regions, the values are less than 600°C·d. The spatial distribution of the freezing index also shows altitude differences, and the index decreases faster in high-elevation regions, such as the northeast and northwest on the Mongolian Plateau, than in low-elevation regions. F D D = 0.18 l o n + 0.94 l a t + 0.49 e l e − 0.31 ( R 2 = 0.93 ) (13)

The overall freezing index of the Mongolian Plateau from 1901 to 2019 shows a downward trend of 4.10°C·d/yr. According to long-term changes (Figure 5A), the freezing index shows a decreasing trend and changes at different rates. Various characteristics of the freezing index are shown in Figure 5C. A decreasing trend appeared in 1902–1903, and a brief increasing trend appears in 1905–1924 and 1929–1930. A decreasing trend with a slow rate of 0.08°C·d/yr was observed from 1930 to 2019.

According to the spatial interpolation map of the thawing index (Figure 6) and multivariate relationship (Eqn. 9), we found that latitude is the main factor and that the thawing index increases gradually with decreasing latitude on the Mongolian Plateau. The thawing index in the northern region varies from 450 to 2,658°C·d, with a mean value of 1,596°C·d, while in the southern region, it varies from 2,000 to 4,000°C·d, with a mean value of 3,216°C·d. High values occur in the Gobi Altai region and the northern part of the Inner Mongolia Autonomous Region of China. During the study period, the region with high values exceeding 3,500°C·d gradually expanded to the south and east. The lowest values appeared in the northern part of the Khuvsgul Mountains, Khangai Mountains, and Khentii Mountains in Mongolia, the Greater Khingan Range in China, and the southern part of Russia. The thawing index varies greatly in the western part. In the northern part of the Inner Mongolia Autonomous Region, the high value region gradually expanded to the east over time. The lowest value gradually increased with time, and the thawing index showed an increasing trend, with significant increases in the Khangai and Khentii Mountain regions located in central Mongolia. T D D = − 0.10 l o n − 0.89 l a t − 0.76 e l e + 1.20 ( R 2 = 0.89 ) (14)

The change in the thawing index on the Mongolian Plateau from 1901 to 2019 is opposite that of the freezing index, with an overall upward trend and a positive slope of 2.25°C·d/yr. According to the M–K mutation test (Figure 5D), the thawing index began to increase abruptly after 1995 (based on the 0.05 and 0.10 significance levels) (Figure 5B). The increasing trend from 1901 to 1995 had a flat slope of 0.44°C·d/yr, whereas the increasing trend from 1995 to 2019 had a significantly steeper slope of 6.19°C·d/yr.

The permafrost on the Mongolian Plateau is mainly distributed in the northern regions, especially in the northern alpine regions (Figure 1A). Sen’s slope can reflect the interannual variation in the freezing and thawing index. The slope change in the air freezing index (Figure 7) shows that the largest freezing index occurs in the northern alpine regions, and the permafrost temperature increases rapidly in the Khuvsgul and Khangai Mountain regions (Kynicky et al., 2009; Zhao et al., 2010; Ishikawa et al., 2018; Munkhjargal et al., 2020; Dashtseren et al., 2021).

The climate determines the air thawing index, and most regions on the Mongolian Plateau exhibit trends consistent with global warming (Munkhjargal et al., 2020; Peng et al., 2020b). Although the change rates of regions listed in Table 1 are based on different time periods and data sources, the thawing index values of the regions show increasing trends (Wu et al., 2011; Luo et al., 2014; Wang et al., 2019; Liao et al., 2021). The influence of latitude on the thawing index in this region is smaller than that in the Arctic Circle, which has simple topographic conditions and is more sensitive to climate warming (Johannessen et al., 2016) than the Mongolian Plateau (Shi et al., 2019). The influence of latitude of the thawing index on the Mongolian Plateau is greater than that on the Qinghai-Tibetan Plateau, which has high elevations and complex topographic conditions (Shi et al., 2019).

The spatial variation in the thawing index of the Mongolian Plateau is also affected by both latitude and altitude (Dashtseren et al., 2021), and the change in the thawing index in the low-latitude region is greater than that in the high-latitude region (Figure 8). The factors that may cause variations in the freezing and thawing index are elevation and solar radiation. Air temperature and solar radiation both gradually decrease with increasing elevation. In addition, the content of water vapor in the atmosphere may also be an influencing factor. For example, in the low-latitude regions with low water vapor content and weak solar radiation, the freezing index changes only slightly. The air freezing and thawing index changes reflect climate warming, which may accelerate permafrost degradation (Dashtseren et al., 2021).

Climate change studies usually require long time series of meteorological data, but the low number and scattered nature of the meteorological stations and boreholes on the Mongolian Plateau make research difficult. Thus, we selected the CRU reanalysis dataset as the data source for analyzing the air freezing and thawing index. However, we found that compared to the observed data, the CRU reanalysis data underestimate the air temperature. This may be caused by the fact that, first, the elevation of a grid cell in the CRU reanalysis is higher than that of the meteorological station. This may result in a difference because of the lapse rate of air temperature (Wang et al., 2018; Shi et al., 2019). Generally, the air temperature decreases with elevation at a rate of 0.65°C/100 m, but there will be differences among regions (Frauenfled et al., 2007; Wang et al., 2018). Second, the CRU reanalysis data are produced by integrating data from meteorological stations worldwide. The original uncertainty is present in the new reanalysis dataset (Harris et al., 2014). Third, the CRU reanalysis dataset still has shortcomings at a spatial resolution of 0.5 × 0.5, as each grid corresponds to actual spatial coverage of approximately 3,000 km² (Frauenfled et al., 2007; Chen et al., 2021), and the spatial representativeness is low. We checked the CRU quality based on a single station, and the deviations were inevitable.

There are also uncertainties in calculating the air freezing and thawing index (Frauenfeld et al., 2007; Chen et al., 2021). Since the CRU reanalysis data based on calculations in this paper provide monthly mean temperature data, the modification of the calculation formula will further increase the uncertainty of the freezing and thawing index results. Moreover, we use temperature data from January 1975 to December 2019 to perform a quality check of the CRU data using a dependence test due to a lack of comparable data, and thus the CRU reanalysis has uncertainty. The underestimation by the reanalysis data undoubtedly affects the air freezing and thawing index results. The air freezing and thawing index are very important for assessing the permafrost status on the Mongolian Plateau (Wu et al., 2011; Dashtseren et al., 2021; Gao et al., 2022). Therefore, improving data accuracy and reducing the uncertainty of the analysis results are of great importance for studying climate change and the permafrost environment (Frauenfled et al., 2007; Harris et al., 2014; Ran et al., 2015).