Abstract
Permafrost degradation induced by climate warming is widely observed in the Northern Hemisphere. However, changes in permafrost sensitivity to climate warming (PSCW) in the future remains unclear. This study examined the changes in permafrost distribution in the Northern Hemisphere under global warming of 1.5 °C and 2 °C, and then characterized the spatial and temporal characteristics of PSCW. Global warming of 1.5 °C and 2 °C would result in 17.8 ± 5.3% and 28.3 ± 7.2% degradation of permafrost area under the climate scenario of Representative Concentration Pathway (RCP) 4.5, respectively, and 18.7 ± 4.6% and 28.1 ± 7.2% under the RCP 8.5, respectively. Permafrost tends to be more sensitive to climate change under the RCP 8.5 than RCP 4.5. PSCW shows small temporal variations in the 21st century under both RCPs, indicating a relatively stable sensitivity to warming on a hemisphere scale. However, PSCW varies greatly among regions, with high values at low latitudes and low values towards high latitudes. Air temperature is a major cause for the spatial heterogeneity of PSCW, explaining 66% of its variations. Permafrost under a warmer climate scenario tends to be more sensitive to the warming. Reducing snow depth and rising air temperature collectively enhances the permafrost sensitivity. Increasing in soil water content, by contrast, reduces the effect of warming. Permafrost in the south of the Northern Hemisphere is most vulnerable to climate warming. Our study highlights that permafrost in the region will respond differently under different warming scenarios across space (e.g. north vs south) and time (e.g. summer vs winter) in this century.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Climate warming is more pronounced in cold regions of high altitudes (Pepin et al
2015) and high latitudes (Cohen et al
2014, Huang et al
2018), leading to widespread permafrost degradation (Romanovsky et al
2010, Schuur et al
2015, Biskaborn et al
2019, Meredith et al
2019). It is estimated (Chadburn et al
2017) that permafrost degradation due to global warming could be at 
million km2 °C−1 (1σ confidence). Permafrost thaw has far-reaching influences on environments and human society in cold regions. The downward moving of permafrost table and the drainage of thawed water can lead to the subside of ground surface, undermining the stability of human infrastructure, and changing landscape and hydrologic processes of surface and subsurface (Liljedahl et al
2016, O'Neill et al
2019). As a result, plant and soil microbial activities and consequent ecosystem processes can also be affected significantly (Wrona et al
2016, Pelletier et al
2019). The interaction of these physical, chemical, and biogeochemical processes in turn affect permafrost thermal regime (Loranty et al
2018). It remains largely uncertain whether these interactions will accelerate or decelerate permafrost degradation caused by climate warming.
Increasing air temperature (Ta) contributes to permafrost degradation dominantly. During 1960–2009, changes in Ta contributed to 84% of the change in permafrost area (Mcguire et al 2016). Many model studies have examined the permafrost sensitivity to climate warming (PSCW). For instance, Koven et al (2013) calculated permafrost degradation in the 21st century using monthly predicted soil temperature with 18 CMIP5 models. They found there is a wide range of permafrost degradation rate, varying from 0.2 to 3.5 million km2 °C−1. Mcguire et al (2016) also compared results of 15 models and found that permafrost area decreasing rate shows great differences among these models from 0.2 to 58.8 × 103 km2 yr−1. The difference in complexity of the physical processes involved in these models and the lack of some critical processes such as thermokarst and ground ice dynamics, are the main reasons for these inconsistencies among models.
Field observations also showed that the influence of warming on permafrost state is uneven around the world. Observations from the circumpolar active layer monitoring (CALM) program showed there is a great spatial heterogeneity of the changing rate of active layer thickness (ALT) (Luo et al 2016a). ALT decreased slightly at five sites, while the other 12 sites experienced distinct increasing trends, varying from 0.05 cm a−1 at Site U1 (Barrow) to 8.4 cm a−1 at Site K0 (Kazakhstan). Even under the same climate in a small basin, Sun et al (2019) concluded that permafrost is more sensitive to warming at low elevations and sunny slopes. The different responses of permafrost (or ALT) to warming indicate that the PSCW varies from sites to sites and regions to regions.
The difficulty in PSCW prediction and the variation of PSCW in observation is largely attributed to that climate warming does not affect permafrost thermal state directly. Heat from the atmosphere must be buffered by surface snow, plant, and soils, before affecting permafrost. Snow aggravates permafrost thaw by preventing cold air from penetrating into soils in winter and therefore keeping permafrost relatively warmer (Zhang 2005). Vegetation does not affect permafrost temperature directly, but it can intercept snow and expose underlying soil to cold winter air. Organic soil layers and humus made of plant litters impact soil thermal properties considerably. Frost shade also has an influence on ground surface temperature (Loranty et al 2018). The effects of water on permafrost temperature are complex. Wet soil has a higher thermal conductivity than dry soil. Therefore, heat conduction in wet soils is higher in summer. However, because of the high heat capacity of water (4200 J kg−1 °C−1), much more heat is needed to increase temperature in wet soils than in dry soils.
Although many studies have examined permafrost dynamics in the past and future, considering the influence of snow, vegetation, and soil properties, few studies have focused on PSCW directly, especially on the global scale. Our study adopted a semi-empirical model to analyze the spatial and temporal variations of permafrost responses to unit temperature change during the 21st century. This study shall increase our understanding of the dynamic responses of permafrost due to climate warming and help identify the most vulnerable permafrost zones to temperature increasing in the future.
2. Data and methods
2.1. Methods
2.1.1. An overview of the GIPL model
We used the GIPL (Geophysical Institute Permafrost Laboratory) model (Sazonova and Romanovsky 2003) to calculate ALT and mean annual ground temperature (MAGT). The core of the GIPL model is a modified Kudryavtsev's approach (Romanovsky and Osterkamp 1997), which has several advantages over the classical Stefan equation (Shiklomanov and Nelson 1999) and allows calculation of ALT and Tps under a wide variety of climate conditions (Anisimov et al 1997). The GIPL model treats air, snow cover, surface vegetation, and active layer as separate layers with different thermal effects on the heat flow from the atmosphere to permafrost. This model has been evidenced to allow estimation of permafrost temperature with high accuracy in comparison with observations (within 0.5 °C) (Sazonova et al 2004), and has accurately estimated ALT in a broad range of permafrost regions, like Alaska (e.g. Shiklomanov and Nelson 1999), Siberia (e.g. Romanovsky et al 2007), and Qinghai Tibet Plateau (e.g. Luo et al 2014). The model is driven by monthly Ta, precipitation, soil water content (SWC), and thermal properties of vegetation, snow cover, and soil layers (refer to table S1 (available online at stacks.iop.org/ERL/16/034038/mmedia) for all model input and output variables). MAGT is the same as mean annual temperature at the top of permafrost (Tps) for regions where permafrost underneath, or mean annual temperature at the bottom of seasonal frozen layer. Tps in GIPL model is given by equation (1):
and ALT (Z) is calculated as equation (2):
where Tgs and Ags are annual mean ground surface temperature and annual amplitude, respectively (°C), accounting for the thermal influence of snow cover and surface vegetation. The calculation of Tgs and Ags are documented by Sazonova and Romanovsky (2003) and have a brief description in the supplementary information. Kt and Kf are soil thermal conductivity in thawed and frozen state, respectively (W m−1 °C−1). C is volumetric heat capacity of thawed ground (J m−3°C−1, refer to Anisimov et al (1997) for calculation), P is the period of the temperature cycle (1 year in second) and Qph is the latent heat of phase change (J kg−1).
The 0 °C threshold of MAGT is critical for permafrost. When MAGT changes from negative to positive, the winter frozen front cannot reach permafrost table by the end of cold season, thus taliks begin to develop and permafrost starts to degrade. It should be noted that the GIPL model focuses on calculation of ALT and MAGT, but it does not calculate the unfrozen depth from ground surface, which includes ALT and the thickness of taliks. For the projection of permafrost area in this study, we excluded the regions where permafrost is already experiencing degradation (MAGT > 0 °C).
There are some shortcomings of the GIPL model. First, it does not take into account unfrozen water and heat flows deeper into permafrost, which can play a critical role in permafrost thermal dynamics (Nicolsky and Romanovsky 2018). Second, ground ice is another important factor that was missing but has significantly impacts on permafrost thermal state (Jorgenson and Osterkamp 2005, Jorgenson et al 2010). Besides, the GIPL model is an equilibrium model based on the assumption of a periodical, quasi-steady-state temperature regime (Kudryavtsev et al 1974). The MAGTs calculated from the GIPL model fluctuate more greatly than measurements, and the difference between observed and calculated MAGTs can be as large as 1.5 °C–2 °C (Sazonova et al 2003). However, compared with transient models which have clear mechanisms but also require much more parameters, initial and boundary conditions, and model inputs that are difficult to acquire on a hemisphere scale, the equilibrium GIPL model demands much fewer model inputs and parameters and fit the needs of simulation on regional to continental and global scales (Anisimov et al 1997, Anisimov and Reneva 2006, Luo et al 2014). Furthermore, when applied to decadal and longer time scales, the GIPL model shows an acceptable degree of accuracy of ±0.2 °C–0.4 °C for the MAGTs and ±0.1–0.3 m for the ALT calculations. In addition, ALT and MAGTs calculated from the GIPL model show a strong correlation (r equal to 0.7–0.9, and p-value close to 0.0) with that calculated from transient models (Sazonova and Romanovsky 2003), demonstrating that the GIPL model is as efficient as transient models in detecting permafrost long-term changes.
In this study, we do not intend to predict the exact permafrost dynamics year by year. Instead, we pay more attention on long-term permafrost stability. The PSCW of each permafrost grid is calculated from all the Ta and ALT from 2006 to 2099 in that grid. Therefore, it actually represents the stability of the permafrost grids on a century scale. The projection of permafrost area under 1.5 °C and 2 °C of global warming is also based on decade average. Taken together, the equilibrium GIPL model shall be suitable for our research objectives.
2.1.2. Soil thermal conductivity
There are two key parameters in the GIPL, including soil thermal conductivity of thawed (Kt) and frozen (Kf) ground. Among many soil thermal conductivity models, the generalized thermal conductivity model developed by Côté and Konrad (2005) provided the best fit between estimation and experimental data (Barry-Macaulay et al 2015). The generalized thermal conductivity model integrates the effects of porosity, degree of saturation, mineral content, grain-size distribution, and particle shape on the thermal conductivity of thawed and frozen soils. Soil thermal conductivity (K) of thawed and frozen ground is given below:
where Ksat and Kdry are thermal conductivity of saturated and dry soils (in thawed or frozen state), respectively (W m−1 °C−1); Kr and Sr are the normalized thermal conductivity and the degree of saturation, respectively; κ is an empirical parameter used to account for the different soil types in the unfrozen and the frozen states (refer to Cote and Konrad (2005) for κ value); ρd and ρw are the density of dry soil and water (kg m−3), respectively; n is soil porosity. Kt and Kf were calculated using these equations separately, based on Ksat, Kdry, and SWC in thawed and frozen state, respectively.
2.1.3. Statistical analysis methods
The geodetector method (Wang et al 2010; www.geodetector.org/) was applied to calculate the contribution of different factors to the spatial distribution of PSCW. The method is based on the concept that the observations can be divided into strata (or categories and zones), within which the values are homogeneous but not between them. The stratified heterogeneity is related to the ratio between the variance within the strata and the pooled variance of the entire study area (Wang et al 2016). q-statistic is used in this method to detect spatial stratified heterogeneity and to measure the association between independent factor (X) and dependent factor (Y), both linearly and nonlinearly. q-statistic is given as:
where L is the number of stratum that Y composed of; N and Nh are the number of data points of Y and the stratum h, respectively; σ2 and σh 2 are variance of the total Y and the stratum h, respectively. The strata of Y are a partition of itself, divided by an explanatory variable X, which should be stratified if it is a numerical variable. In this study, we stratified X (e.g. Ta) into 10 categories (L = 10) through the K-means method recommend by Wang and Xu (2017). The q-statistic can vary from 0 to 1, which means X explains q*100% of Y; q = 0 indicates that there is no coupling between Y and X; and q = 1 indicates that Y is completely determined by X. With definite physical meaning and no linear hypothesis, this method has already been applied in many fields of natural and social sciences (Wang et al 2010, Luo et al 2016b, Zhang et al 2019, Yang et al 2020).
2.1.4. Calculation of permafrost sensitivity to climate warming
PSCW in this study was defined as the deepening of ALT upon 1 °C of warming. When plotting ALT against Ta, the slop of the curve is PSCW. To explore the temporal variation of PSCW in the 21st century, we plotted ALT against Ta in the three periods of 2010–2019, 2050–2059 and 2090–2099 (represent the beginning, middle and end of the 21st century), respectively. The temporal variation of PSCW was then deduced from the differences among the relation between ALT and Ta during these three periods. Spatially, we obtained PSCW for each permafrost grid by constructing a linear fitting equation for yearly ALT and Ta for 2006–2099 in that grid, and the slope of the fitted line is defined as PSCW for that grid. During the simulation, if Tps for a model grid exceeds 0 °C, which means the surface permafrost is already thawed and a talik starts to form (Sazonova and Romanovsky 2003), the calculated ALT for the grid at this time was set to null value (because the change of ALT in this case does not directly relate to the thermal dynamics of underneath permafrost below taliks). PSCW for the grid is then calculated from the remain ALT and Ta pairs. For model grids where there are few ALT and Ta pairs remaini (less than 20 pairs), we did not calculate PSCW for them because small number of data pairs may cause great uncertainty during curve fitting and thus cannot get a robust PSCW. We also calculated the p-value (probability that the regression coefficient is zero) for the regression equation in each permafrost grid. If p-value in a model grid is greater than 0.05, we excluded PSCW for that grid.
2.2. Data
2.2.1. Climate data for model simulations
Climate data of monthly near surface Ta, snowfall and soil moisture content were derived from four global circulation models (GCMs) in the second simulation round of the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP 2b). Known issues of previous round of ISI-MIP have been solved for the ISI-MIP 2b through a series of adjustments, and atmospheric GCM data provided by the ISI-MIP 2b have also been bias-adjusted to a new reference dataset of EWEMBI (Frieler et al 2017). The four GCMs (IPSL-CM5A-LR, GFDL-ESM2M, MIROC5, and HadGEM2-ES) provided a similar fractional range coverage to that of randomly chosen four-member sets of CMIP5 GCMs (Frieler et al 2017). All climate variables have the same spatial resolution at 0.5° × 0.5°, and the daily datasets were averaged to monthly. The temporal coverage of the datasets was 1986–2099. We used the historical (1986–2005) reconstruction climate data to drive our model and then validated the model results of permafrost distribution and ALT against previous studies and observations. For future (2006–2099) projections, Representative Concentration Pathway 4.5 (RCP 4.5) and RCP 8.5, which correspond to radiative forcing levels of 4.5 and 8.5 W m−2 by 2100, respectively (Moss et al 2010), were selected to evaluate the effects of climate change on PSCW. Data for the Southern Hemisphere were excluded.
It has been estimated that 1.5 °C and 2 °C warming threshold would be reached by about 2030 and 2050 under the RCP 4.5, respectively, and about 2025 and 2040 under the RCP 8.5, respectively (Donnelly et al 2017, Karmalkar and Bradley 2017, Li et al 2019). To account for the uncertainties among different GCMs that were used to estimate the time threshold of 1.5 °C and 2 °C global warming, we adopted a 11 year span around each threshold year, that was 2025–2035 and 2045–2055 for 1.5 °C and 2 °C warming under the RCP 4.5, respectively, and 2020–2030 and 2035–2045 for 1.5 °C and 2 °C warming under the RCP 8.5, respectively.
2.2.2. Snow depth, vegetation, and soil data for model simulation
Snow and vegetation property data are needed for the calculation of Tgs and Ags. Snow depth was calculated by the equation given by Nelson and Outcalt (1987) from monthly precipitation of the ISI-MIP 2b dataset. Snow thermal properties in this study were assumed constant over time and throughout the study area. ρsn, Csn, and λsn were set equal to 300 kg m−3, 0.32 W m−1 °C and 5.07 × 10−7 m2 s−1, respectively (Shiklomanov and Nelson 1999). Although these snow parameters are realistic for the whole Northern Hemisphere, they have much smaller effect on the model's overall accuracy than snow depth (Shiklomanov and Nelson 1999). Because of the lack of vegetation thermal property data at the hemispheric scale, we assumed a 10 cm layer of moss across the study area with thermal diffusivity equal to 1.39 × 10−6 and 5.56 × 10−8 m2 s−1 in frozen and thawed states, respectively (Kudryavtsev et al 1974). The 10 cm thick layer of moss has been proved to achieve the best fit against data from tundra regions by Anisimov et al (1997).
Thermal conductivity of saturated and dry soils (Kwet and Kdry, respectively), and the density (ρd) and heat capacity (Cdry) of dry soil are needed during the calculation of Kf, Kt and C. The typical value for Kwet, Kdry, ρd, and Cdry of sand, silt, clay, and peat soil are given by Anisimov et al (1997) and Shiklomanov and Nelson (1999), showing good model performance when compared with observations (Anisimov et al 1997, Shiklomanov and Nelson 1999). For a given model grid, we first obtained the soil texture data for the grid from the Harmonized World Soil Database (v 1.2; FAO/IIASA/ISRIC/ISSCAS/JRC, 2012). Based on the soil texture data, we obtained the proportion of sand, silt, clay, and peat for that grid. By multiplying the proportion and corresponding thermal property and then summing up the results, we obtained the thermal property for that grid.
The lack of spatially explicit snow and vegetation thermal parameters on the hemisphere scale induces uncertainties in our model estimation. Additionally, this study did not consider the impacts of vegetation dynamics on permafrost degradation, which can be significant (Jorgenson et al 2010, Loranty et al 2018). Although enduring these uncertainties and the model shortcomings, by regressing ALT to Ta on a century scale, some fluctuations caused by these factors might have been filtered out, resulting in a robust estimation of the long-term PSCW trend.
2.2.3. Observed ALT data for model validation
Observed ALT from CALM network sites (www2.gwu.edu/~calm/data/north.htm) were used to validate the model results during the historical period (1986–2005). Three primary methods are employed at CALM sites. Some sites use spatially oriented mechanical probing at regular intervals across a grid (like 100 m × 100 m) and/or transect(s). Some sites employed thaw-tube measurements and some other sites infer thaw depth from ground temperature measurements. ALT is measured at the end of the thawing season each year. For spatially oriented measurements, ALT of the site is calculated by averaging all the sampling points at the site, while for point observations (e.g. thaw-tube and ground temperature measurements), the single measurement is ALT for this site. In this study, we mainly chose the sites that employ the spatially oriented mechanical probing method to get a more general representation of the local value. The earliest ALT observation at CALM network sites can date back to 1990 and continue into present. But for the purpose of model validation of ALT during the historical period, we selected sites that have sufficient (no less than 5) observations during 1990–2005. Sixty-five CALM sites (black dot in figure 1(a)) around the Northern Hemisphere permafrost regions were finally selected. For the model validation, we first averaged the observed and calculated ALT from 1990 to 2005 to get the mean value for the historical period. Then located these 65 CALM sites into model grids and validated the model by comparing the historical mean ALT of the observations and model results in corresponding model grids.
Figure 1. (a) Simulated permafrost distribution in the Northern Hemisphere under recent past climate (1986–2005) by GIPL model; the red lines denote the boundary of discontinuous and sporadic permafrost regions Brown's (2002) map. The black dots represent the location of the 65 observations from CALM, with averaged measured versus calculated ALT during 1995–2005 for each site shown in (b). The horizontal and vertical error bars in panel (b) give the ±1 standard deviation for measured and calculated ALT during 1990–2005, respectively.
Download figure:
Standard image High-resolution image3. Results
3.1. Model validation
Based on ISI-MIP 2b historical climate data, permafrost distribution in the Northern Hemisphere was simulated for the recent past (1986–2005, figure 1(a)). Permafrost extent in this study only includes regions where MAGTs are negative (Tps < 0). According to Obu et al (2019), the 0 °C of MAGTs roughly corresponds to the boundary of the discontinuous and sporadic permafrost regions. It can be seen from figure 1(a) that the simulated permafrost regions well match the continuous and discontinuous permafrost distribution (Brown et al 2002). Permafrost area was estimated to be 14.9 ± 0.2 × 106 km2 (figure 1(a)), which is within the range of 12.2–17.0 × 106 km2 given by Zhang (2000), but a little larger than 13.9 × 106 km2 estimated by Obu et al (2019) using the same indicator of negative MAGTs at the top of permafrost. But because the study period of Obu et al (2019; from 2000 to 2016) is later than the historical period (1986–2005), considering the climate warming from 1986 to 2016, the larger extent of our study could be reasonable.
The comparison between measured and calculated ALT was shown as figure 1(b). The modeled ALT is generally higher than observed values, suggesting that the model might have slight warming biases.
3.2. Permafrost degradation in the 21st century
All the four GCMs predicted a consistent warming during the 21st century under the two RCP scenarios (figure 2). The change of permafrost distribution is closely related to Ta changes. Under the RCP 4.5 scenario, when warming slows down after about 2070 for the IPSL-CM5A-LR and MIROC5 models, the decreasing in permafrost area would also decelerate. While under the RCP 8.5, persistent warming would lead to permafrost thaw continually. Among different models under the same RCP (figure 2), stronger warming would normally result in larger permafrost extent decrease, and vice versa.
Figure 2. Changes in the Northern Hemisphere permafrost area and Ta in permafrost regions under the RCP 4.5 and RCP 8.5 scenarios, respectively.
Download figure:
Standard image High-resolution imageWe also estimated permafrost distribution when global mean Ta rises by 1.5 °C and 2 °C relative to pre-industrial levels (1850–1900), and permafrost distribution at the end of the 21st century (2090s) under the RCP 4.5 and RCP 8.5 scenarios, respectively (figure 3 and table 1). When Ta reaches the 1.5 °C warming threshold, 2.66 ± 0.82 and 2.80 ± 0.71 millionkm2 of permafrost would thaw under the RCP 4.5 and RCP 8.5, respectively. Another 1.57 ± 0.30 and 1.40 ± 0.42 million km2 of permafrost would thaw under the RCP 4.5 and RCP 8.5, respectively, when 2 °C warming threshold reached. Under both global warming threshold of 1.5 °C and 2 °C, permafrost degradation rate is greater under the RCP 8.5 than RCP 4.5 (table 1), indicating a higher sensitivity of permafrost extent to climate warming under the RCP 8.5.
Figure 3. Distribution of the Northern Hemisphere permafrost when global mean near surface Ta rise by 1.5 °C and 2 °C above pre-industrial levels (1850–1900) and at the end of 21st century under the RCP 4.5 and RCP 8.5 scenarios respectively. Also shown in thick red lines are the boundary of the permafrost map provided by Brown et al (2002).
Download figure:
Standard image High-resolution imageTable 1. Projected permafrost (PF) extent in the Northern Hemisphere when global mean near surface Ta rise by 1.5 °C and 2 °C above pre-industrial levels (1850–1900), and at the end of the 21st century.
| 1.5 °C | 2 °C | End of the 21st century | ||||
|---|---|---|---|---|---|---|
| RCP 4.5 | RCP 8.5 | RCP 4.5 | RCP 8.5 | RCP 4.5 | RCP 8.5 | |
| PF degradation (×106 km2) | 2.66 ± 0.82 | 2.80 ± 0.71 | 4.23 ± 1.13 | 4.20 ± 1.13 | 5.73 ± 1.85 | 11.33 ± 2.39 |
| Remain fraction | 0.82 ± 0.05 | 0.81 ± 0.05 | 0.71 ± 0.07 | 0.72 ± 0.07 | 0.61 ± 0.12 | 0.24 ± 0.15 |
| Warming (°C) a | 1.85 ± 0.52 | 1.89 ± 0.54 | 2.91 ± 0.81 | 2.83 ± 0.84 | 3.95 ± 1.37 | 8.05 ± 1.92 |
| PF degradation per degree (×106 km2 °C−1) | 1.43 ± 0.13 | 1.48 ± 0.09 | 1.45 ± 0.07 | 1.48 ± 0.07 | 1.47 ± 0.08 | 1.42 ± 0.05 |
Note. The second number in each cell is the standard deviation (std) for the four GCMs. a Denotes warming in permafrost regions relative to the regional mean Ta of 1986–2005.
By the end of the 21st century, Ta in permafrost regions would rise by 3.95 ± 1.37 °C relative to present level (1986–2005) under the RCP 4.5, and the temperature sensitivity of permafrost distribution was projected at 1.47 ± 0.08 million km2 °C−1 during the 21st century. As a result, 39.3 ± 11.9% (5.73 ± 1.85 million km2) of the present permafrost extent (14.9 × 106 km2) would experience degradation. Under the RCP 8.5 scenario, only 24.1 ± 15.2% of the current permafrost area would persist at the end of the 21st century, due to a much greater warming of 8.05 ± 1.92 °C above present level. These permafrost regions would mainly distribute in northern Canada and in the east of Siberia inside the Arctic Circle (figure 3).
3.3. Temporal and spatial variation of PSCW
To analyze the temporal and spatial variation of PSCW, we plotted ALT against mean annual, summer (May–October), and winter (November–April in the next year) Ta for all the Northern Hemisphere permafrost grids in three periods of 2010–2019, 2050–2059, and 2090–2099, respectively (figure 4). We used these three periods to represent the begin, middle, and end of the 21st century. Because of the degradation of permafrost, there are less data points for the period of 2090–2099 (as shown in figure 4), especially for the RCP 8.5. On the hemisphere scale, ALT increases exponentially with mean summer Ta, leading PSCW, the slope of the fitted curves, increasing with mean summer Ta. However, the relationship between ALT and Ta is not clear in winter, and therefore resulting in a weak relationship for all year around. These results indicate that ALT is sensitive to temperature increase in summer but resistant to winter warming.
Figure 4. Scatter plot of ALT against mean annual, summer (May–October), and winter (November–April in the next year) Ta of all Northern Hemisphere permafrost grids during the three periods of 2010–2019, 2050–2059, and 2090–2099 under the RCP 4.5 and RCP 8.5, respectively.
Download figure:
Standard image High-resolution imageThe differences in the slopes of the fitted curves at low and high summer Ta show that PSCW has great spatial heterogeneity in the Northern Hemisphere. As for temporal variation, there are little differences between the fitted curves of the three periods under the RCP 4.5. Although under the RCP 8.5, PSCW increases slightly during 2050–2059, compared to spatial heterogeneity, the differences of PSCW in the three periods are relatively small, indicating that the temporal variation of PSCW is relatively stable on the hemisphere scale.
The spatial distributions of PSCW based on mean annual, summer and winter Ta under the RCP 4.5 and RCP 8.5 are shown in figure 5, respectively. For all the distributions, PSCW are greatest in low latitudes and decrease with increasing latitudes. Furthermore, PSCW based on summer temperature is generally higher than that based on mean annual temperature, while PSCW based on winter temperature is the lowest. This phenomenon demonstrates again that ALT is more sensitive to summer warming. It should also be noted that although PSCW based on winter temperature is generally small, temperature increasing in winter still deepens ALT, especially for low latitude permafrost regions under the RCP 8.5 scenario.
Figure 5. Distribution and histogram of PSCW (m °C−1) based on mean annual, summer, and winter Ta in the 21st century under the RCP 4.5 and RCP 8.5 scenarios.
Download figure:
Standard image High-resolution imageAs for the difference in PSCW under the two RCP scenarios, ALT tends to be more sensitive to warming (on matter it is annual, summer, or winter warming) under the RCP 8.5 (the histograms of figure 5), since more permafrost grids have higher PSWC under the RCP 8.5.
4. Discussion
4.1. Sensitivity of ALT and PSCW
We used a simplified method to examine the sensitivities of ALT and PSCW to Ta, seasonal air temperature amplitude (Aa), SWC, and snow depth (Zsn), respectively. Take Ta as an example, first, we run the model with the original data of the RCP 4.5 scenario. Then we kept variables of Aa, SWC and Zsn unchanged but increased Ta in all the model grids during the whole study period by 0.5 °C at each step. The responses of ALT and PSCW to increasing temperature then can be plotted out (as shown in figure 6). The sensitivities of ALT and PSCW to Aa, SWC, and Zsn were also tested respectively (figure 6) using the same method (but with different steps: 0.5 °C for Aa, 0.02 for SWC, and 0.04 m for Zsn). Because in some regions, permafrost thaws quickly before the sensitivity study finishing (especially for increasing Ta), making the calculation of ALT and PSCW based on different permafrost extent, thus the comparison between different steps is not ideal. Therefore, we selected eastern Siberia (60°–75° N, 120°–150° E), which persists during all the simulation periods, as a representative region for this sensitivity study.
Figure 6. Sensitivity study of ALT (m) and PSCW (m °C−1) for a representative permafrost region in Siberia (60°–75° N, 120°–150° E). The first dot (from left to right) in each panel is the calculated ALT or PSCW from original RCP 4.5 data. The dots (from left to right) after the first one was calculated from gradually increased Ta (°C), Aa (°C), SWC (m3 m−3), or Zsn (m), respectively. Ta and Aa increase by 0.5 °C at each step; SWC increase by 0.02 at each step; and Zsn increase by 0.04 m at each step.
Download figure:
Standard image High-resolution imageBoth ALT and PSCW increase quickly with Ta (figure 6). The influence of Ta on ALT is direct, through deepening ALT in summer and increasing soil temperature and unfrozen SWC in winter (Li et al 2018). The increasing of PSCW with Ta means that ALT would deepen greater under warmer climate, suggesting that ALT is more sensitive to Ta in warmer climate. This can result from several processes. First, soil temperature would be relatively higher in a warm winter (Hansen et al 2014), which can make the frozen ground easier to thaw during the coming spring and summer. Second, a warmer summer would transfer more heat into soil and compel permafrost table to move to a lower level. Moreover, freezing front moves shallower underground in a warmer winter. When freezing front cannot reach the permafrost table, talik will form between them and help thaw the surface frozen ground form bottom to up in next summer (Walvoord and Kurylyk 2016).
ALT increased gradually with Aa even when Ta remained unchanged (figure 6). However, the influence of Aa on PSCW is limited. The increasing of Aa with Ta remaining unchanged means a warmer summer but a colder winter. A warmer summer under greater Aa would lead to the deepening of ALT, but a cold winter would decrease the soil temperature and unfrozen water content, which may improve the resistance of the frozen ground to temperature increasing in the next spring and summer and somewhat offset the warming effects on PSCW in summer, resulting in a limited influence on annual scale.
ALT and PSCW decrease with SWC noticeably (figure 6). Because of the high volumetric heat capacity of water, soil volumetric heat capacity will increase with SWC significantly. Therefore, it needs more heat to increase 1 °C of temperature in wet soils than dry ones, which makes wet soils more resilient to warming.
Zsn has important effects on permafrost thermal regime (Zhang 2005, Park et al 2015, Bisht et al 2018). Because a large fraction of snow layer is filled with air, snowpack has extremely low thermal conductivity, insulating thermal transfer between the atmosphere and substrate ground. With the increasing of snow depth, insulation effects of snowpack will be stronger than albedo effects (Zhang 2005), warming the ground and therefore increasing both ALT and PSCW (figure 6).
4.2. Temporal stability of PSCW
Figure 4 shows that although PSCW would change over time, their differences among the three periods (i.e. the three fitted lines) are relatively small. When plotting mean ALT against Ta in the Northern Hemisphere permafrost regions during 2006–2035, 2036–2065 and 2066–2095 (figure 7), we find that PSCW (the slop of ALT against Ta) are similar during these three periods, no matter based on mean annual Ta, or mean summer and winter Ta.
Figure 7. Relationship between mean ALT and mean annual, summer and winter Ta in three 30 year periods under the RCP4.5 and RCP 8.5, respectively.
Download figure:
Standard image High-resolution imageThe dynamics of PSCW can be controlled by heat balance between the atmosphere and ground, and local factors such as vegetation, snow covers, soil thermal properties, surface water and groundwater. As for the temporal variation of a specific region, factors like topography and soil texture can be treated as constant values in the 21st century. And since vegetation dynamics in the 21st century is out of the scope of this study and may introduce additional uncertainty, we fixed vegetation properties and focused mainly on the dynamics of Ta, SWC, and Zsn.
During the 21st century, Aa of the Northern Hemisphere permafrost region would decrease by 1.33 ± 0.94 °C and 3.29 ± 1.74 °C under the RCP 4.5 and RCP 8.5, respectively (figure S1). SWC and Zsn would increase by 0.02 ± 0.0 and 0.01 ± 0.00 m, respectively, under the RCP 4.5, and 0.03 ± 0.02 and 0.02 ± 0.0 m, respectively, under the RCP 8.5. The sensitivity study (figure 6) shows that the magnitude of the changes in SWC and Zsn under both RCPs would have limited influence (less than 0.002 m °C−1) on PSCW. PSCW is sensitive to Ta (figure 6). Figure 7 shows that PSCW tends to increase with summer warming under the RCP 8.5. However, since climate warming in permafrost regions would be more significant in winter (increasing by 4.42 ± 2.36 °C and 11.25 ± 3.89 °C under the RCP 4.5 and RCP 8.5, respectively) than in summer (increasing by 2.17 ± 1.21 °C and 6.09 ± 1.66 °C under the RCP 4.5 and RCP 8.5, respectively), and since PSCW is resistant to winter warming as described above, the relationship between ALT and Ta during the three periods on annual scale shows little variation (figure 7). Apart from this, although the decrease in Aa and the increase in SWC have limited influence on PSCW, both of them contribute to decreasing PSCW. As a result, PSCW would keep relatively stable on the hemisphere scale in the 21st century under both RCPs.
4.3. Spatial heterogeneity of PSCW
Figure 5 shows that PSCW would decrease significantly from low latitudes to high ones. We plotted the main impact factors of PSCW along latitude in figure 8. It can be seen that Ta shares the same spatial pattern of PSCW, while Aa and Zsn have opposite trends. SWC peaks at 55°–60° N and then decreases slightly. Geodetector method results (table S2) show that Ta could explain more than 50% of the spatial variation of PSCW, which is much greater than that of Aa, SWC and Zsn (about 10%, 14% and 14%, respectively). The variation of PSCW explained by Ta in this study is smaller than that in Mcguire et al (2016), which concluded that changes in Ta explain 84% of the change in permafrost area during 1960–2009. The difference may result from different study periods (historical vs future), but the predominant role of Ta on permafrost dynamics is clear.
Figure 8. Variation of annual mean Ta, Aa, SWC and Zsn along latitude in Northern Hemisphere (a); PSCW along latitude (b). 'All' indicates Ta, Aa, SWC and Zsn all varied along altitude as RCP 4.5 scenario; 'Ta' indicates only Ta varied along altitude as RCP 4.5, Aa, SWC and Zsn were set as their mean value of the latitude between 45° and 50° N; 'Ta & (Aa, SWC or Zsn)' indicates Ta and one of the three variables of (Aa, SWC or Zsn) changed along altitude while the other two variables were equal to their mean value of the latitude between 45° and 50° N.
Download figure:
Standard image High-resolution image
We further tested the sensitivity of PSCW to the spatial variation of Ta, Aa, SWC and Zsn. First, we set values of Aa, SWC and Zsn in all the pixels equal to their mean values between latitudes 45° and 50° N, respectively, but kept Ta the same as the original RCP 4.5 data. We run the model with these modified data and found that the calculated PSCW followed the original results (when Ta, Aa, SWC and Zsn all changes along latitude) (figure 8). This indicates again that Ta is the main reason for the spatial heterogeneity of PSCW, which is consistent with our geodetector results.
In order to find out whether the spatial variation of Aa, SWC and Zsn would strengthen or weaken the effects of Ta on PSCW, we let one of the three variables (Aa, SWC and Zsn) change with Ta along altitude, but keep the other two equal to their mean values between latitudes 45° and 50° N. Figure 8 shows that when Ta and Zsn change together, PSCW is noticeably larger than that when Ta changes alone. This result demonstrates that Zsn would promote the influence of Ta on PSCW. However, SWC tends to have an opposite effect. Increased SWC in high latitudes tends to weaken ALT temperature sensitivity. Aa increases slightly with latitude and has weak regulation on the influence of Ta on PSCW.
5. Conclusions
PSCW analysis is critical for understanding the dynamic responses of permafrost to climate warming. The spatial and temporal variation of the simulated PSCW in the Northern Hemisphere during 2006–2100 were analyzed. We found that, the 21st century would witness a wide permafrost degradation, starting from the south boundary towards high latitudes. Permafrost degradation would follow climate warming closely. Global warming of 1.5 °C and 2 °C relative to 1860–1900 would lead to 3.26 ± 0.89 and 5.01 ± 1.27 million km2 of permafrost degradation, respectively, under the RCP 4.5, and 3.63 ± 0.74 and 6.40 ± 1.40 million km2 of permafrost degradation, respectively, under the RCP 8.5. Permafrost is more sensitive to climate change under the RCP 8.5 than RCP 4.5.
PSCW shows small temporal variation on the hemisphere scale during the 21st century under both RCPs. The relatively small change in Aa, SWC, and Zsn and the mainly warming in winter (which has little influence on PSCW) may result in the temporal stability of PSCW. However, PSCW shows a great spatial heterogeneity, peaking at low latitudes and reducing towards high latitudes. Ta is the main factor affecting the spatial heterogeneity of PSCW, explaining more than 50% of its spatial variation. Permafrost in warmer climate tends to be more sensitivity to climate warming. Snow depth would interact with Ta to promote PSCW. But the regulation of Ta on PSCW tends to be weakened towards wet permafrost regions. The influence of Aa on the regulation of Ta on PSCW is relatively small.
Acknowledgments
This work was supported by projects of the 'Strategic Priority Research Program' of the Chinese Academy of Sciences (Grant No. XDA20020202) and the National Natural Science Foundation of China (Grant Nos. 41571193 and 41530749). The research is also logistically supported by the Department of Earth, Atmospheric, and Planetary Sciences at Purdue University. We gratefully acknowledge financial support from China Scholarship Council. We thank three anonymous reviewers for their critical comments and suggestions. We also thank the ISIMIP cross sectoral science team for their roles in producing, coordinating, and making available the ISI-MIP 2b data (available at: https://esg.pik-potsdam.de/search/isimip/). The authors declare no conflict of interest.
Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.