Attribution of historical near-surface permafrost degradation to anthropogenic greenhouse gas warming

Given the current confirmed permafrost degradation and its considerable impacts on ecosystems, water resources, infrastructure and climate, there is great interest in understanding the causes of permafrost degradation. Using the surface frost index (SFI) model and multimodel data from the fifth phase of the Coupled Model Intercomparison Project (CMIP5), this study, for the first time, investigates external anthropogenic and natural forcing impacts on historical (1921–2005) near-surface permafrost change in the Northern Hemisphere. The results show that anthropogenic greenhouse gas (GHG) forcing produces a significant decrease in the area of near-surface permafrost distribution at a rate of 0.46 × 106 km2 decade−1, similar to observations and the historical simulation (ALL). Anthropogenic aerosol (AA) forcing yields an increase in near-surface permafrost distribution area at a rate of 0.25 × 106 km2 decade−1. Under natural (NAT) forcing, there is a weak trend and distinct decadal variability in near-surface permafrost area. The effects of ALL and GHG forcings are detectable in the observed change in historical near-surface permafrost area, but the effects of NAT and AA forcings are not detected using the optimal fingerprint methods. This indicates that the observed near-surface permafrost degradation can be largely attributed to GHG-induced warming, which has decreased the near-surface permafrost area in the Northern Hemisphere by approximately 0. 21 × 106 km2 decade−1 on average over the study period, according to the attribution analysis.


Introduction
Twenty-two percent of the land surface of the Northern Hemisphere is underlain by permafrost (Obu et al 2019). Approximately 11 000-37 000 km 3 of ground ice  and 1330-1580 Pg of organic carbon (Schuur et al 2015) are stored in permafrost, which is defined as ground that is frozen for at least two consecutive years. Permafrost degradation, accompanied by ground ice melting and organic carbon release, could have adverse effects on ecosystems (Vonk et al 2015), water resources (Liljedahl et al 2016, Pastick et al 2019, infrastructure (Nelson et al 2001, Hjort et al 2018, and climate (Koven et al 2011). In addition, permafrost contains a large pool of the harmful substance mercury (Schuster et al 2018) and various bacteria (D'Costa et al 2011). Their emergences due to permafrost thaw may have serious consequences for human health.
Given the considerable potential impacts of permafrost degradation, growing research communities have focused on estimating the amount of nearsurface permafrost degradation using observations and simulation methods (Romanovsky et al 2010, Lawrence et al 2012, Luo et al 2016, Aalto et al 2017, Chadburn et al 2017, Biskaborn et al 2019. Despite distinct uncertainties due to sparse monitoring stations, as well as the variable climate sensitivities and soil processes among the models, there is good agreements among different studies that the extent of areas with the presence of near-surface (within ∼3-4 m, IPCC 2019) permafrost has shrunk during recent decades and will probably continue to shrink for the next one hundred years Lawrence 2013, Biskaborn et al 2019).
To understand the processes of permafrost degradation, some studies have attempted to explore the factors that induce the degradation (Zhang 2005, Shkolnik et al 2010. Air temperature is considered to be the most dominant control factor for permafrost change through affecting ground temperature (Biskaborn et al 2019). In addition, snow properties and their seasonal evolution have significant control on soil thermal conditions (Zhang 2005). A recent study showed that an increase in summer precipitation can result in the expansion of near-surface permafrost extent (Guo and Wang 2017). These studies have shed light on the local factors that induce the thaw of permafrost. But it is not well understood how permafrost is influenced by large-scale external anthropogenic (greenhouse gases and aerosols) and natural (volcanic aerosols and solar variability) forcing agents, which is important for in-depth understanding of the mechanisms of permafrost change.
Detection and attribution provides a rigorous approach to quantify the relative contributions of internal climate variability and different external forcings to the observed climate changes (Barnett et al 2005, Huber andKnutti 2011). Detection and attribution was initially applied to global temperature patterns (Stott et al 2000). Subsequently, detection and attribution studies have focused on global or regional mean precipitation patterns (Lambert et al 2005, Zhang et al 2007, surface climate extremes (Stott et al 2004, Zhang et al 2013, Sun et al 2014, Fisher and Knutti 2015, Arctic sea ice (Min et al 2008), the hydrological cycle (Barnett et al 2008), and snow cover (Rupp et al 2013). To date, no detection and attribution studies of Northern Hemisphere permafrost have been published.
This study aims to explore the influence of external anthropogenic and natural forcings on historical permafrost change in the Northern Hemisphere using the optimal fingerprint detection and attribution method. The surface frost index (SFI) model was used to diagnose the permafrost from climate models from the fifth phase of the Coupled Model Intercomparison Project (CMIP5).

Data
The historical climate data used to diagnose permafrost were derived from CMIP5 simulations. They include monthly 2 m air temperature, snow depth, and snow mass (used to calculate snow density).
Historical simulation (ALL) and three external forcing experiments were used: greenhouse gas (GHG) forcing, natural (NAT) forcing, and anthropogenic aerosol (AA) forcing. Pre-industrial control (CTL) experiments were used for detection and attribution analysis. Climate models that include all the five of the above experiments were selected. Accordingly, 11 models were selected and their key details are shown in Tables S1 and S2. The models provided monthly 2 m air temperature, snow depth, and snow mass data. However, three models (GFDL-CM3, GFDL-ESM2 M, and IPSL-CM5 A-LR) did not provide snow depth and snow mass data. For these three models, their precipitation data were used to calculate mean snow depth in winter (details can be seen in the following section 2.2.2). Additional details with respect to the CMIP5 simulations can be found in (Taylor et al 2012). The data are available at http://pcmdi9.llnl.gov/esgf-web-fe/, and they have been widely employed for climate change evaluation studies (IPCC 2013).
The Climatic Research Unit (CRU) highresolution monthly 2 m air temperature and precipitation version 4.01 data were used to drive permafrost model to diagnose the observed permafrost, and were obtained from http://www.cru.uea.acuk/. The precipitation data were used to calculate mean snow depth in winter. The data cover the period from 1901 to 2016 and have a resolution of 0.5 • × 0.5 • . The CRU data are developed based on station observations using the climate anomaly approach. Station anomalies are resampled into a high-resolution grid, and then combined with an existing climatology to produce absolute values. For further details of the CRU dataset, refer to (Harris et al 2014). CRU is one of the best-known gridded observation datasets and it has been used worldwide to detect historical climate change (IPCC 2013).
As one of the best available datasets on the distribution of permafrost, the International Permafrost Association (IPA) map was used as the observed permafrost extent to validate the simulations (Brown et al 1997). The data are available at http://nsidc. org/data/docs/fgdc/ggd318_map_circumarctic/index .html at a resolution of 0.5 • × 0.5 • . The data include continuous, discontinuous, isolated, and sporadic permafrost and glaciers in the Northern Hemisphere. Previous studies have indicated that climate model (e.g. Community Climate System Model version 3) generally can identify only continuous and discontinuous permafrost due to its coarse spatial resolution (1.4 • × 1.4 • ) (Burn and Nelson 2006). This judgement is resolution-dependent, which may be not true for higher-resolution model. In this study, a simulation resolution of 0.5 • × 0.5 • was used, which still may be not sufficiently high to identify isolated and sporadic permafrost. Thus, we used only continuous and discontinuous permafrost observations for the validation.

Model (surface frozen index)
The surface frozen index (SFI) model was used to diagnose near-surface permafrost extent. The model was developed by (Nelson and Outcalt 1987) as where DDF is the sum of freezing degree days. The label ' * ' means a consideration of snow insulation effect that results in a reduction in the DDF (Nelson and Outcalt 1987). DDT is the sum of thawing degree days that are calculated in terms of the sinusoidal climate. These two variables can be calculated using specific equations that are described in detail in (Nelson and Outcalt 1987). The model requires four climate variables as input data: 2 m air temperature in the warmest and coldest months, and mean snow depth and snow density in winter. In this study, the 2 m air temperature in the warmest and coldest months were calculated using monthly 2 m air temperature data. (Nelson and Outcalt 1987) indicated that the use of such monthly 2 m air temperature causes an underestimation of freezing or thawing indices. Mean snow depth and snow density in winter were calculated in terms of a weighted mean of monthly snow data with monthly 2 m air temperature below 0 • C, as shown in (Slater and Lawrence 2013). This method intends to obtain more meaningful mean snow data in winter for expressing the snow insulation effect than the arithmetic mean. However, for the three models and CRU data that did not provide snow data, their precipitation data were directly used to calculate mean winter snow depth and a snow density of 250 kg m −3 was used (Slater and Lawrence 2013). The SFI model emphasizes the importance of climate in forming permafrost. In addition, the snow insulation effect on permafrost is also considered. Because climatic stationarity is implicitly assumed in the model, climate variables are averaged over a certain period before they are input into model. Similar to (Slater and Lawrence 2013), this study diagnoses the permafrost in a certain year (e.g. 2000) using the previous 20 year average (i.e. in this case, 1981-2000) climate variables. We examine that diagnosed decadal variability in near-surface permafrost area does not depend on the timescale (e.g. 20 or 10 years) chosen to average climate variables to meet the climatic stationarity, although shorter timescale in general gives more obvious inter-annual variation (not shown). The SFI values range from 0 to 1 and a criterion of SFI > 0.6 is used to determine the continuous and discontinuous permafrost in this study (Guo and Wang 2016). Note that the SFI model actually indicates the sustainability of the upper (near-surface) permafrost layer under a certain stationarity climate condition. The superiority of the SFI model is embodied in the requirement of readily available climate data and easy and rapid calculation, which makes it suitable to use for diagnosing permafrost distribution during a long-term period under various climate forcings (Nelson and Outcalt 1987). A previous study indicated that the SFI model can provide more information for the evaluation of change in permafrost than raw diagnoses using soil temperature from the current climate models (Slater and Lawrence 2013).

Method for calculating mean winter snow depth using precipitation data
For the three models (GFDL-CM3, GFDL-ESM2 M, and IPSL-CM5 A-LR) and CRU data that did not provide snow data, the method introduced by (Nelson and Outcalt 1987) was used to calculate the mean winter snow depth using monthly precipitation data. This method weights snowfall by its duration, i.e. earlier snowfall in winter has greater weighting than later snowfall in spring. Furthermore, the fact that the magnitude and duration of snow thaw changes with altitude is taken into account.
whereZ s is mean winter snow depth, P i is precipitation in the i th month (i = 1, 2,…, k) when the mean temperature is ⩽ 0 • C. ρ r is snow density (250 kg m −3 , (Slater and Lawrence 2013). ∅ is the latitude of the site.

Optimal fingerprint method
The optimal fingerprint method is used to perform detection and attribution analysis of permafrost in this study, which includes the following equations (Ribes et al 2013): where Y is the observations. X i is the permafrost response to the i th external forcing (e.g. GHG, NAT, and AA), β i is the scaling factor that adjusts the fingerprint value to yield the best match to the observations, n is the total number of external forcings, and ε is the internal climate variability, which is calculated using the CTL simulations. Here X i is assumed to be unknown andX i is computed from the ensemble mean of multimodel.X i involves a residual term, ε x i , which is related to internal variability. The corresponding residual consistency test is used to check the reliability of the residual term (Allen and Stott 2003). We divide the CTL simulations into 77 chunks that denote non-overlapping 85 year samples. To avoid spurious detections, we use half of these chunks to calculate the regression coefficients and use the rest to  [1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005] near-surface permafrost areas (green area) to observations (areas outlined in blue, made ∼1960-1990) and (c) differentiation of the simulated present-day near-surface permafrost areas from 11 climate models. In panel (a), the simulated near-surface permafrost area has been diagnosed by the SFI model driven by the CRU data. In panel (b), the simulated near-surface permafrost area has been diagnosed by the SFI model driven by the CRU data climatology plus multimodel ensemble-mean anomalies. In panel (c), the simulated near-surface permafrost area has been diagnosed by the SFI model driven by CRU data climatology plus anomalies of each model. The unit of the color bar is the total number of models that captured the same near-surface permafrost area (e.g. gray indicates the regions where there are 11 models capturing near-surface permafrost distribution). Bias between CRU data-based and simulated and observed near-surface permafrost areas is given in the bottom right corner of the panels (a) and (b).
carry out a residual consistency test. If the uncertainty range (5%-95%) of the scaling factor is greater than zero, the fingerprint of the corresponding external forcing is detectable in the observations (Min et al 2011). In addition, if the uncertainty range also contains '1' , the simulation is considered to be in agreement with observations (Min et al 2011). The attributable trends in permafrost area for a certain external forcing are calculated as the linear trends of permafrost area of this external forcing multiplied by the corresponding 5%-95% scaling factors (Sun et al 2014, Li et al 2017.

Ensemble empirical mode decomposition
The ensemble empirical mode decomposition (EEMD) method (Wu and Huang 2009) is used to separate the original series of permafrost area into variability (interannual and interdecadal) and trend terms. EEMD is an adaptive and temporally local sparse decomposition method and it decomposes the data series, X (t), into a suite of components from high frequency to low frequency, C j , together with the residual term of R n .
The principle of EEMD is that the added white noise is uniformly distributed in the whole timefrequency space that consists of components of different scales. In this study, the output of EEMD includes seven components. The first component is the original series, the sum of the second to sixth components is variability, and the last component is the nonlinear trend.

Other methods
To reduce systematic biases in the climate data from CMIP5 simulation, yearly anomalies of each climate variable (2 m air temperature, snow depth and snow density) were first calculated relative to the 1901-1920 and they were then added to CRU data climatology for 1901-1920. The anomalies of 2 m air temperature (snow depth and snow density) denote a change (proportional change). This method has been used in the previous researches on permafrost change Lawrence 2013, Guo andWang 2016).
This study uses simulations and gridded observations, but they have different horizontal resolutions. To perform homogenous calculations and comparisons, all data are sampled into a common resolution of 0.5 • × 0.5 • . In addition, the linear trend of permafrost area is calculated as the slope of the linear fit based on least squares regression (Guo et al 2018).

Model validation
Multimodel ensemble-mean present-day (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) near-surface permafrost extent fits well with observations, except for slight overestimations in the Labrador Peninsula in northeastern Canada and the eastern Tibetan Plateau as well as underestimations in southern Alaska and the northern Western Siberian Plain ( figure 1(b)). The ensemble-mean near-surface permafrost area is 17.2 × 10 6 km 2 , which is close to observed area of 15.2 × 10 6 km 2 , with a bias of 2.0 × 10 6 km 2 . Differences in permafrost extent between the models are small and they are mostly distributed at the southern edges of the permafrost region ( figure 1(c)). Note that, to date, there are not observed temporal series of change in permafrost area for the Northern Hemisphere. In this study, observed CRU air temperature and precipitation are used to drive the SFI model to obtain a series of near-surface permafrost area from 1921 to 2005. Our validation shows that simulated present-day permafrost extent based on CRU data is in good agreement with the observations from the IPA map (Brown et al 1997), with a bias in area of 1.8 × 10 6 km 2 ( figure 1(a)). Thus, we take this CRU data-based series of permafrost area as the 'observed series of near-surface permafrost area' in this study.
In addition to the IPA map, we also compare the simulated near-surface permafrost distribution to the results in Alaska from (Pastick et al 2015) based on multiple data fusion approach. Main features of the simulated distribution are very comparable to those from (Pastick et al 2015). Such as, all two results show that the near-surface permafrost is mostly present in the northern and eastern parts of Alaska, but absent in the southern part of Alaska. These similarities add the reliability of our simulation. It is unsurprising that some deviations are shown between the two results, such as, (Pastick et al 2015) shows that near-surface permafrost is absent in Brooks Range in Alaska, in contrast to our results. These may be due to that (Pastick et al 2015) refers to the presence of permafrost in the upper 1 m, which differs from the nearsurface (in general the upper ∼3-4 m) in this study.

Near-surface permafrost response to external anthropogenic and natural forcings
From 1921 to 2005, the ensemble-mean near-surface permafrost area profile first decrease until approximately mid-1960 s and then slightly increase until approximately mid-1980 s with a fast decrease until 2005 under the ALL forcing, which is very similar to the CRU data-based profile ( figure 2(a)). The ensemble-mean trend of near-surface permafrost area under the ALL forcing is −0.16 × 10 6 km 2 decade −1 , which is very close to −0.18 × 10 6 km 2 decade −1 of the CRU data-based trend ( Figure S1 (available at stacks.iop.org/ERL/15/084040/mmedia)). Near-surface permafrost area significantly decreased under GHG forcing from 1921 to 2005, with an ensemble-mean trend of −0.46 × 10 6 km 2 decade −1 ( figure 2(b)). Near-surface permafrost area increased by a rate of 0.25 × 10 6 km 2 under the AA forcing (figure 2(d)), which is suppressed by the significant decrease resulting from the anthropogenic GHG for- cing. Under NAT forcing, there is a weak trend and distinct decadal variability in the near-surface permafrost distribution area (figure 2(c)).
Change in permafrost extent mostly occurs at the edge of near-surface permafrost region, particularly in the northern part of the Western Siberian Plain in Russia and the southern edge of the permafrost in Canada, similar to the results from a previous study by (Guo and Wang 2017) based on a numerical land surface model ( figure 3). The spatial pattern of nearsurface permafrost change from the GHG forcing is similar to that from the CRU data and the ALL forcings, but shows an obviously larger magnitude of change. Both GHG and ALL forcings produce an decrease in all permafrost regions, except for a few grids located in the northern and eastern edges of Greenland that show an increase. The AA forcing produces an increase in all permafrost regions, except for a few grids located in the northern and eastern edges of Greenland that shows a decrease. The NAT forcing produces a weak decrease in the northern part of the Western Siberian Plain in Russia, whereas an increase is found in the southern edge of the permafrost in Canada.

Detection and attribution of historical near-surface permafrost degradation
Based on the optimal fingerprint method, the scaling factors of the near-surface permafrost area series under different forcings are calculated and are shown in figure 4. The standard residual consistency test shows that the null hypothesis that the observed near-surface permafrost area series is equivalent to ensemble-mean results from all forcings is established at the 90% confidence level (all the p values larger than 0.1). This indicates that the regression models well fit the data. In the single-signal analysis, the scaling factor of ALL forcing is 1.10 and the associated 90% uncertainty range (0.71-1.43) is greater than zero, indicating that the ALL simulation fits well with the observations and its signal is detectable in permafrost in the Northern Hemisphere. For the other three external forcings, the GHG signal also can be detected, with a scaling factor of 0.46, and the associated 90% uncertainty ranges is 0.28-0.65. The best estimate of the scaling factor is smaller than '1' , indicating an overestimation of near-surface permafrost change. The aforementioned results regarding temporal series and spatial patterns of near-surface permafrost change in section 3.2 support this analysis. However, the NAT-only and AA-only signals are not detected, with their scaling factors smaller than zero. Similar results can be obtained using the threesignal analysis that is considered to be able to separate the individual forcing contributions from the combined effects of the three forcings ( Figure S2). These analyses suggest that changes in near-surface permafrost area in the Northern Hemisphere are related to anthropogenic GHG warming. When this conclusion is suggested, one would want to ask why permafrost apparently loses in the southern edge of permafrost region but the greatest warming associated with GHGs is at higher latitudes (IPCC 2013). This is because permafrost at the southern edge of permafrost region has higher soil temperature (closer to the melting point) (Biskaborn et al 2019) and the associated larger sensitivity to climate warming than permafrost in higher-latitude regions (Guo and Wang 2016). The original series of permafrost area from the observations and ALL, GHG, NAT, and AA forcings are further separated into variability and trend items ( Figure S3). Both ALL and GHG signals can be detected in the observed trend and variability items, while both NAT and AA signals cannot be detected ( Figure  S4). This indicates that anthropogenic GHG warming is as an important factor to influence not only the trend but also the variability of the near-surface permafrost area in the Northern Hemisphere.
The attributable trends of near-surface permafrost change to anthropogenic GHG warming are further analyzed (figure 4). In general, the trend of changes in near-surface permafrost area from observations diagnosed based on CRU data is −0.18 × 10 6 km 2 decade −1 (−0.16 to −0.21 × 10 6 km 2 decade −1 ) from 1921 to 2005. The trend of near-surface permafrost area from the one-signal GHG forcing analysis is −0.21 ×10 6 km 2 decade −1 (−0.13 to −0.30 × 10 6 km 2 decade −1 ), which is consistent and comparable to the observed change. This indicates that anthropogenic GHG warming has contributed to the decrease in near-surface permafrost area during recent decades.

Discussion and conclusions
The GHG forcing produces a significant decrease in near-surface permafrost area, while the AA forcing yields a significant increase in near-surface permafrost area. The GHG forcing signal is detectable in the observed near-surface permafrost change. However, the AA forcing cannot be detected and the negative scaling factors of AA forcing imply an opposite effect from anthropogenic aerosols, which partially offsets the decreasing role of the GHG forcing in permafrost change. The NAT forcing produces a weak trend but distinct decadal variability of near-surface permafrost area and its signal also cannot be detected. The observed near-surface permafrost change can be attributed to anthropogenic GHG warming, which is estimated to have caused a decrease in near-surface permafrost area of approximatel −0.21 × 10 6 km 2 decade −1 (−0.13 to −0.30 × 10 6 km 2 decade −1 ) on average during recent decades.
Further analyses indicate that changes in the 2 m air temperature in the permafrost region correspond well to those in the near-surface permafrost area from different external forcing experiments. Specifically, ALL and GHG forcings produce a significant increasing air temperature, AA forcing produces a significant decreasing air temperature, and NAT forcing produces a weak change in air temperature (Figures S5 and S6). Despite this, snow depth in the permafrost region shows very weak change under all four forcings ( Figures S5 and S6). This indicates that the external forcings control the near-surface permafrost extent through mainly influencing the 2 m air temperature.
Possible uncertainties in this study could be due to the coarse resolution (0.9 • × 1.3 • to 2.8 • × 2.8 • ) of atmospheric forcing data from the climate model data providing less regional climate information, which may prevent the SFI from capturing detailed change in permafrost at the edge of the permafrost area. In addition, the observed near-surface permafrost series is obtained using the CRU observations to drive the SFI model in this study. The series may deviate from the true value to some extent and thus contribute to the possible uncertainties in detection and attribution analyses. Sporadic and isolated permafrost are easy to thaw due to their warm heat condition. However, the coarser resolution (0.5 • × 0.5 • ) of the model could not identify sporadic and isolated permafrost. Therefore, the changes in these two type of permafrost may not be included in the present study.
Historically, research on permafrost has moved from historical and future permafrost changes to the impact of permafrost change on climate, engineering facilities, and ecosystems. The present study, for the first time, shed light on the influences of external anthropogenic and natural forcings on permafrost change in the Northern Hemisphere during recent decades. The results advance the understanding to the mechanisms of permafrost change in the Northern Hemisphere.