Permafrost is warming at a global scale

Permafrost warming has the potential to amplify global climate change, because when frozen sediments thaw it unlocks soil organic carbon. Yet to date, no globally consistent assessment of permafrost temperature change has been compiled. Here we use a global data set of permafrost temperature time series from the Global Terrestrial Network for Permafrost to evaluate temperature change across permafrost regions for the period since the International Polar Year (2007–2009). During the reference decade between 2007 and 2016, ground temperature near the depth of zero annual amplitude in the continuous permafrost zone increased by 0.39 ± 0.15 °C. Over the same period, discontinuous permafrost warmed by 0.20 ± 0.10 °C. Permafrost in mountains warmed by 0.19 ± 0.05 °C and in Antarctica by 0.37 ± 0.10 °C. Globally, permafrost temperature increased by 0.29 ± 0.12 °C. The observed trend follows the Arctic amplification of air temperature increase in the Northern Hemisphere. In the discontinuous zone, however, ground warming occurred due to increased snow thickness while air temperature remained statistically unchanged.

O ne quarter of the Northern Hemisphere and 17% of the Earth's exposed land surface is underlain by permafrost 1 , that is ground with a temperature remaining at or below 0°C for at least two consecutive years. The thermal state of permafrost is sensitive to changing climatic conditions and in particular to rising air temperatures and changing snow regimes [2][3][4][5][6][7] . This is important, because over the past few decades, the atmosphere in polar and high elevation regions has warmed faster than elsewhere 8 . Even if global air temperature increased by no more than 2°C by 2100, permafrost may still degrade over a significant area 9 . Such a change would have serious consequences for ecosystems, hydrological systems, and infrastructure integrity [10][11][12] . Carbon release resulting from permafrost degradation will potentially impact the Earth's climate system because large amounts of carbon previously locked in frozen organic matter will decompose into carbon dioxide and methane [13][14][15] . This process is expected to augment global warming by 0.13-0.27°C by 2100 and by up to 0.42°C by 2300 15 . Despite this, permafrost change is not yet adequately represented in most of the Earth System Models 14 that are used for the IPCC projections for decision makers. One major reason for this was the absence of a standardized global data set of permafrost temperature observations for model validation.
Prior to the International Polar Year (IPY, 2007(IPY, -2009, ground temperatures were measured in boreholes scattered across permafrost regions. However, a globally organized permafrost data network and a standard reference period against which temperature change could be measured did not exist. One key outcome of the IPY was strenghtening the Global Terrestrial Network for Permafrost (GTN-P) 16,4 . This initiative established a temperature reference baseline for permafrost and led to an increase in the number of accessible boreholes used for temperature monitoring.
To analyze the thermal change of permafrost we assembled a global permafrost-temperature data set that includes time series of data attributed to the IPY reference boreholes. We compiled a time series for the decade from 2007 to 2016 that comprises mean annual ground temperatures T, determined from temperatures measured in boreholes within the continuous and discontinuous permafrost zones in the Arctic (including the Subarctic), Antarctica and at high elevations outside the polar regions. The measurements were made at, or as close as possible to the depth of zero annual amplitude Z * , where seasonal changes in ground temperature are negligible (≤0.1°C). Rates of permafrost temperature change calculated for the 2007-2016 decade were indexed in each borehole to suppress near-surface and deep geothermal changes. Regional and global change rates were calculated as area-weighted means. To compare single borehole sites, due to the higher availability of full-year records after 2007, we ranked the temperature difference between the biennial means of 2008-2009 and 2015-2016. We used linear regression on T between 2007-2016 to estimate decadal change rates. To calculate annual departures, we compared consecutive years to the reference mean of 2008-2009. We concluded, that ground temperature near the depth of zero annual amplitude increased in all permafrost zones on Earth, that is continuous and discontinuous permafrost in the Northern Hemisphere, as well as permafrost in the mountains and in Antarctica. The observed trend followed increased air temperature and snow thickness, each in varying degrees depending on the region.

Results
Permafrost temperature changes. Measurements from borehole sites established prior to the IPY generally indicated warming driven by higher air temperatures ( Fig. 1)  Similar to warming of the Arctic continuous permafrost zone, the Antarctic permafrost warmed by 0.37 ± 0.10°C dec À1 Ref . However, the remoteness of the continent and its limited accessibility resulted in far fewer boreholes drilled to Z * compared to the Northern Hemisphere. Consequently, permafrost temperature departures and trends were statistically not significant and had large uncertainty bands (Fig. 3d).
Air temperature changes. The relation between air and soil temperature development in permafrost regions is not straightforward due to highly variable buffer layers such as vegetation, active layer soils, or snow cover. To compare permafrost temperature changes to those in the atmosphere, we applied the same calculation method for each borehole site using mean annual air temperatures (T) at 2-m height above ground level (Fig. 4a, d), spatially interpolated from the ERA Interim gridded reanalysis data set 23 . We calculated general snow thickness changes for Arctic sites in Fig. 4a, b. However, there is not, as yet, a reliable consistent data set on snow thickness applicable for high elevation regions or Antarctica.
The propagation of temperature change in the atmosphere downward to the depth of Z * can take up to several years, but the time varies depending on the surface characteristics, the subsurface ice content, and the soil thermal diffusivitiy 24,25 . We took this lag into account by averaging over the previous 4 years for each year considered, but there was no significant correlation at an annual resolution between permafrost temperature departures at Z * depth and 2-m air temperature anomalies derived from ERA Interim data alone (Fig. 4). This lack of correlation can be attributed to the discrepancy between the scale at which borehole observations are conducted and the spatial resolution of 80 km for the gridded airtemperature reanalysis data 26 and because in permafrost regions, the reanalysis output is more dependent on the model structure and data assimilation methods than in data-rich regions 27 ; local microand secondary climate effects 28 ; and buffering layers at the airground interface 5 that influence the thermal response of permafrost to short-term changes in air temperature.
Previous studies have shown that these surface effects, along with the thermal diffusivity of the underlying materials, act as a buffer that reduces the effect of short-term climate   Fig. 5b), however, do not match the observed strong permafrost warming. This discrepancy is due to large climatic differences between the Antarctic Peninsula and eastern Antarctica 30,31 , the small number of boreholes that fulfill the quality criteria, and the principal climate model bias in Antarctica 32 .
Air temperature trends in the Arctic continuous permafrost zone correspond well with permafrost temperature change rates (Figs. 3a and 4a), suggesting that enhanced warming of permafrost in the High Arctic reflects the polar amplification of recent atmospheric warming 22 . However, in the Arctic discontinuous permafrost zone, air temperatures were statistically unchanged between 2006 and 2014 while permafrost temperatures increased. We found that snow dynamics, the time lag between air and ground temperature, and the latent heat effect serve as concurrent explanations for this phenomenon.
Snow thickness changes. The snow cover reduces the upward transfer of energy from the ground to the air during winter 33,34 . Distinct peaks in the mean snow depth in 2009, 2011 and from 2013 onward (Fig. 4a, b) suggest that the observed continued warming of discontinuous permafrost is facilitated by increasing snow thickness. Compared to the Arctic continuous permafrost zone, the mean snow cover in the discontinuous zone arrived about 1 week later, reached its maximum insulation 1 month earlier, and also disappeared half a month earlier. Compared to 2007-2009 the snow cover in 2014-2016 in the discontinuous zone started to form 13.7 days earlier, reached its maximum insulation effect 37.7 days earlier, and disappeared 9.3 days earlier (Fig. 4f). It was shown previously that a difference of only 10 days caused significant warming in Alaska 35 . Increases of shrub height and density that trap wind drifting snow is likely also a contributing factor 36 . All of these changes provide evidence of increased protection of the ground from low temperatures during winter 37,38 . Snow timing differences within the continuous zone are less distinct but show a generally similar trend (Fig. 4e, f).

Discussion
An important factor that explains the general discrepancy between mean annual temperature changes at Z * in permafrost and the atmosphere is that permafrost progressively with depth "remembers" the surface temperature history of the past several years 25,39 . The temporal dimension of episodes with lower air temperatures between 2009 and 2013 in the Arctic (Fig. 4a, b), and around 2012 in the mountains (Fig. 4c), relative to preceding period of higher air temperatures, however, was not large enough to sustainably impact the general warming trend of permafrost.
We partly attribute the difference in ground temperature change between the continuous permafrost and the discontinuous permafrost zones to the latent heat effect. In this process, the icewater phase change associated with warmer permafrost in the   (Fig. 2a, b) reduces the response of ground temperature to changes in air temperature 4 . Cold permafrost therefore exhibits a greater response to changing air temperature compared to permafrost with a temperature close to 0°C 4,40 .
The warming of permafrost observed since IPY continues the trends documented prior to IPY 41 . Our global analysis suggests that the future increases in air temperature projected under current climate scenarios 42 will result in continued permafrost warming. The duration of our time-series, however, does not yet permit predictive analysis of non-linear climate-permafrost relations as the latent heat effect is stronger near 0°C and surface characteristics are not constant. However, observations of thaw at some of the observation sites demonstrate that the latent heat requirement cannot indefinitely delay permafrost warming down to depths of about 15 m observed in this study (Fig. 6), nor prevent the eventual thawing of permafrost. This could have wide implications in terms of permafrost degradation and release of greenhouse gases from decomposition of organic matter.
The SWIPA 2017 report 41 gave an estimate of 0.5°C warming of permafrost in very cold areas such as the High Arctic since IPY (2007)(2008)(2009)). This is similar to our network observations of strong warming within the Arctic continuous permafrost zone and of continued warming elsewhere. The assessment of permafrost temperature trends presented in this paper can facilitate validation of models to project thawing of permafrost down to the depth of Z * and associated impacts with respect to feedbacks to the climate system.
The current global coverage of permafrost temperature monitoring is not yet ideal, due to the limited sampling in regions such as Siberia, central Canada, Antarctica, and the Himalayan and Andes mountains. Furthermore, even though the data used were quality checked and are as complete as possible, logistical challenges during fieldwork caused gaps in the time series. Better assessments of the evolution of the thermal state of permafrost, including consideration of non-linear system behavior, will benefit from ongoing efforts to enhance the global network spatially and extend the length of the record. Enhancing existing monitoring sites through co-location with meteorological stations could further improve understanding of microclimate and bufferlayer influences, and would also provide the data necessary for a comprehensive assessment of permafrost responses to ongoing climate change.
The newly compiled GTN-P data set has facilitated assessment of trends in permafrost temperatures and can also contribute to improved representation of permafrost dynamics in climate models and the reduction of uncertainty in the prediction of future conditions.

Methods
Field observations of permafrost temperatures. Boreholes were established and temperatures were recorded during annually repeated fieldwork campaigns in polar and high-elevation areas. Temperature was measured either by lowering a calibrated thermistor into a borehole, or recorded using permanently installed multisensor cables 43 . Measurements were recorded either manually with a portable temperature system or by automated continuous data logging. At some borehole sites, permafrost thawed at the measurement depth during the period of observations. The criterion to include non-permafrost sites in the global change calculation was that ground temperatures near the depth of the Z * were below 0°C until the end of the IPY reference period in 2009. They are then transferred to a global data set after a 1-year embargo to allow authors to publish their local findings first. Within the GTN-P Data Management System the data presented were harmonized, quality checked and filtered to generate a standardized global permafrost borehole data set. Data standardization was performed during data entry into the database following international geospatial metadata standards ISO 19115/2 and TC/221. The data management system is based on an object-oriented data model, accessible online at http://gtnpdatabase. org. The GTN-P mean annual ground temperature T compilation is accessible online at https://doi.org/10.1594/PANGAEA.884711.
A total of 154 boreholes with 1264 T values were used in this study. Data analyses of decadal permafrost temperature change were based on 123 boreholes and 1033 T values calculated from > 10 5 sensor observations.
Calculating permafrost temperature change. We used the R environment 44 to calculate the mean permafrost temperature change for every borehole from qualityfiltered T data. The same measurement depth was used each year for a borehole. The depth was chosen to be the nearest available sensor to the depth of Z * , the depth at which seasonal changes in temperature are ≤0.1°C (Fig. 7). The nearest depth to Z * was detected by either an algorithm calculating the difference between annual maximum (summer) and minimum (winter) temperature in the original data starting from the shallowest depth downwards and using cubic spline interpolation between thermistors and a threshold set to sensor accuracy, or by visual inspection of annual maximum and minimum temperature measurements plotted versus depth (Fig. 7). Because the depth of Z * varies over time as temperature changes, we used an average estimated for the observation period. The data revealed that 19.5% of measurements were from above Z * . 59.8% of measurements represented Z * and 20.7% were from below Z * . Measurements from boreholes that had no reliable indication of Z * had a mean depth of 17.1 m, which is well below the average of all indicated Z * values (mean 14.1 m, median 12 m). Thus, the data distribution represents an approximation to Z * which minimizes the potential bias caused by seasonal fluctuations.
We created a data set that reflects long-term climate change and avoids large temperature fluctuations caused by seasonal phenomena, e.g., in Antarctica, by excluding data from shallow boreholes that did not reach Z * . Because Z * could not be determined in all boreholes the minimum depth was set to 10 m. However, five boreholes with depths between 6.7 m and 10 m were included (GTN-P ID's 16 : 137, 860, 861, 877, and 1192), because their depths were equal to Z * , and seasonal fluctuations were less than the instrument precision and accuracy. Boreholes that fulfilled the quality criteria but were not included in this analysis due to depth constraints, represented 22.6% of the original data set. 8.6% were excluded from the Arctic continuous data set; 23.4% from the Arctic discontinuous data set; 30.0% from the mountain data set; and 57.1% from the Antarctic data set. Statistically indifferent temperature trends of the remaining shallow (≤12 m) and deeper (>12 m, max. 40 m) boreholes in the utilized data set confirm that the observed depths near Z * (Fig. 6b) provide a representative sample tracking climate variability coherently.
We applied different methods to extract information on permafrost temperature changes in single years, in single boreholes and for decadal changes in the permafrost regions, described as follows: The last term on the right-hand side of Eq. (1) serves as our mean value for the reference period. We compare this reference period to the latest available mean value period and calculateΔ T b to rank total temperature differences among boreholes.
Equations (1) and (2) require data to be available in each of the observation years.
To calculate the rate of temperature change per decade we follow a third approach using the primary mean annual ground temperature data set T b for all available years in i and perform linear regression, according to the following attribution of our data in the regression equation: where T reg b is the regression estimate of T b , a b is the vertical intercept (the starting temperature in a borehole), c b is the slope of the regression line, and x is the range of years involved.
The requirement to perform linear regression on b was that i included at least one value y in the IPY period (2007, 2008, or 2009), one value in the modern reference period (2015 or 2016) and a minimum of five values in total. We calculated the rate of temperature change in each borehole as the slope of the linear regression c b using the linear model function (lm) in the R environment. To generate decadal change values, we extrapolated 37.7% of the borehole data in the Arctic continuous zone, 47.3% in the Arctic discontinuous zone, 29.3% in the mountain zone and 100% in Antarctica for 1-3 years.
The consistency of temperature time series in boreholes depends on sustained data collection at remote sites. At some boreholes, instrumentation was destroyed, damaged or malfunctioned leading to interruptions in data collection 45 . To avoid broken data runs affecting the annual means, measurements at frequencies greater than monthly (e.g. daily or hourly), were aggregated to monthly means before calculating annual means. Mean annual values were based on at least monthly primary data. Data points based on fewer than one measurement every month were allowed only if the sensor depth was equal to or below the depth of zero annual amplitude. Annual means were calculated from original measurements as calendaryear means in the GTN-P Database. Meteorological years in permafrost areas depend on the onset and termination of the freezing and thaw periods, and in previous studies varied spatially. We therefore indicated the starting month of the period in the data set. Mean values contain only the available valid T data in each year, and thus the number of borehole temperatures included in change-rate calculations varies between years.
To evaluate temperature changes in the Arctic continuous and discontinuous permafrost zones, in the mountain permafrost and in permafrost in Antarctica, we applied a spatial de-clustering prior to calculating mean values of temperature changes from the boreholes. The spatial de-clustering reduces the bias in the calculation of means caused by an inhomogeneous (clustered) spatial distribution of the boreholes. We grouped the boreholes into ten world zones (Fig. 8) and defined the areas underlain by permafrost by correlating the boreholes with the International Permafrost Association (IPA) permafrost zones 46 . Arctic continuous permafrost represents the mean of four different zones: Arctic continuous permafrost West (2.41 × 10 6 km 2 ), Arctic continuous permafrost West islands (1.57 × 10 6 km 2 ), Arctic continuous permafrost Europe (0.22 × 10 6 km 2 ), and Arctic continuous permafrost East (Asia) (6.62 × 10 6 km 2 ). Arctic discontinuous permafrost is averaged over three zones: Arctic discontinuous permafrost West (3.91 × 10 6 km 2 ), Arctic discontinuous permafrost East (Asia) (3.86 × 10 6 km 2 ), and Arctic discontinuous permafrost Europe (0.28 × 10 6 km 2 ). Mountain permafrost is averaged over two zones: Chinese mountains (2.07 × 10 6 km 2 ), and Other mountains (2.33 × 10 6 km 2 ) including the Alps and other sites with high elevations >1000 m a.s.l. such as in Scandinavia and the North American Cordillera. Antarctica is treated as one zone (0.05 × 10 6 km 2 6,47 ). For comparing temperature trends between North American and north Asian permafrost we define two separate data sets by excluding southern, European, and central Asian boreholes. Within the zones, clusters of boreholes close together were grouped when the sum of longitude and latitude differences were <0.1 decimal degree and the T values of adjacent boreholes were averaged before calculating the mean temperature change.
To estimate the mean annual temperature change in each zone we applied areaweighted arithmetic averaging of T values in boreholes. To preserve the signal of local outlier trends showing atypical temperature change directions and magnitudes (e.g., in parts of Antarctica and in Québec, Canada), we did not use medians. To suppress near-surface and geothermal changes indices of boreholes were distributed as three possible integers to multiply the sites before averaging, according to the following criteria: (i) T is available in each year of the reference periods indicated in Eqs. (1) and (2), and (ii) T depth is equal to the depth of Z * and >10 m (few exceptions were made according to the depth of Z * as described above). Calculating air temperature change. The set of air temperature data monitored at borehole sites is incomplete. To develop data comparable to the permafrost temperature data, we calculated mean annual air temperatures (T) from ERA Interim 2 m air temperature data set with 80 km spatial resolution. We derived the reanalysis time series for each borehole from linear interpolation of the four nearest grid points surrounding the borehole coordinates. Mean annual values were calculated from December until November. Given, that the propagation of atmospheric temperature change downward to the depth of Z * takes up to several years 25,37 , depending on the local thermal diffusivity 24 , we extended the time series shown in Fig. 4 backwards to 2000 and used the standard reference period 1981-2010 to estimate anomalies.
We define a set j = {1981,...,2016} to identify the years being considered. We use the coordinates of boreholes b defined in Eq. (1) and calculate the annual difference for specific years y 2 j inT as Based on the average propagation of surface temperature towards Z * of 4 years 25 we calculated 4-year end-point running means to compare air temperature with permafrost temperature changes. To calculate the rate of temperature change over a decade, we apply linear regression onT y;b for all y 2 j using the linear model function in the R environment and the slope of the linear regression in an annual array between 2004 and 2016 and multiplied the annual change rates by 10. Data analyses of air temperature change were based on 137 borehole sites and 4932T values.
Calculating snow thickness change. We calculated the mean annual snow thickness (Ŝ) for the Arctic continuous and the discontinuous permafrost zone from the Canadian Meteorological Centre (CMC) daily snow depth analysis data with 24 km spatial resolution 48 . We derived the reanalysis time series for each borehole from linear interpolation of the four nearest grid points surrounding the borehole coordinates. Mean values were calculated from December until February for each year in the data set. To identify winters we use subsequent years, e.g. in the time series we assign the 1999-2000 winter to 2000. Given that 1999 is the earliest available year in the data set we define a set k = {1999,...,2016} to identify the winter years, where y 2 k. We use the coordinates of boreholes b defined in Eq. (1) and calculate the annual difference inŜ as ΔŜ y;b ¼Ŝ y;b À 1 12 The onset snow has an impact on the ground thermal regime. To assess the onset of snow cover, we assemble a set of snow depths dates at daily resolution between 1 September and 30 April in a set of days l = {1,2,3,...,242} for every year in k. In leap years l = {1,2,3,...,243}. To calculate the onset date of snow SO we use the first day d SO k;b reaching 6 cm in l for which the following 5 days, adding up to a synoptic time scale of 6 days, retain a daily snow cover of at least 6 cm 49 .
The insulation maximum of snow SIM is reached when the snow cover accumulated to a thickness between 40 and 50 cm 33,37 . Accordingly, we set SIM based on the first day d SIM k;b in l reaching 50 cm, or, if it is not reached, take the day representing the maximum snow cover in l (below 50 cm).
To assess the end of snow cover SE, we assemble the snow depth dates at daily resolution between 1 September and 30 August in a set m = {1,2,3,...,365} for every year in k. In leap years m = {1,2,3,...,366}. To calculate SE we use the first day d SE k;b in m reaching down to less than 1 cm after a decreasing gradient of at least 8 cm over 6 days, or, if this gradient is not reached, the first day of at least 6 subsequent snow free (<1 cm) days.
Measurement accuracy. The reported measurement accuracy of our temperature observations, including manual and automated logging systems, varied from ±0.01 to ±0.25°C with a mean of ±0.08°C. Previous tests have shown the comparability of different measurement techniques to have an overall accuracy of ±0.1°C 3 . Thermistors are the most commonly used sensors for borehole measurements. Their accuracy depends on (1) the materials and process used to construct the thermistor, (2) the circuitry used to measure the thermistor resistance, (3) the calibration and equation used to convert measured resistance to temperature, and (4) the aging and resulting drift of the sensor over time. Thermistors are typically calibrated to correct for variations due to (1) and (2). About 20% of the boreholes are visited once per year and measured at or below Z * using single thermistors and a data logger. In this case the system is routinely validated in an ice-bath allowing correction for any calibration drift. The accuracy of an ice-bath is~± 0.01°C 50 . Using the offset determined during this validation to correct the data greatly increases the measurement accuracy near 0°C, an important reference point for permafrost. The remaining systems are permanently installed and typically ice-bath calibrated at 0°C before deployment. The calibration drift is difficult to quantify as thermistor chains are not frequently removed for re-calibration or validation. In many cases removal of thermistor chains becomes impossible some time after deployment, e.g. because of borehole shearing.
The drift rate among bead thermistors from different manufacturers was <0.01°C per year during a 2 year experiment at 0, 30, and 60°C 51 . The calibration drift of glass bead thermistors was found to be 0.01 mK per year 52 , at an ambient temperature of 20°C. A single drifting thermistor in a chain is detectable through its anomalous temporal trend. Such data were excluded from our data set. The absolute accuracy of borehole temperature measurements, in terms of their representativeness of the temperature distribution in undisturbed soil, also depends on the depth accuracy of the sensors' positions in the borehole. This study is concerned with temperatures at Z * , where temperature gradients are typically small (<0.1°C m −1 ). Consequently, mm-level positioning accuracy does not significantly impact measurement accuracy. Finally, as this study is concerned with annual averages, adequate chronometry is ensured.
The above discussion of accuracy relates to the absolute temperature values measured, but the detection of temperature change is more accurate because errors in calibration offset have no impact, sensor nonlinearities are generally small and not of concern. We therefore consider <0.1°C a conservative average estimate of the accuracy of temperature change on an individual sensor basis.
Confidence intervals and statistical significance. Permafrost and air temperature departure from 2008 until 2016 (Δ T i;b and ΔT y;b ) and the regression from 2007 until 2016 of each borehole were used to calculate the 95% confidence intervals within each world zone using a Student t-test in the R environment (52% p < 0.05, 48% p > 0.05, mean |t| = 3.4  were calculated from de-clustered and indexed boreholes. Mean confidence intervals for composite permafrost zones (global, Arctic continuous, Arctic discontinuous, mountain, Asian and American) were area-weighted. Antarctica consists of one zone and thus area-weighting is not applicable. Given a nonnormal, unimodal, and only slightly skewed distribution of data in similarly shaped subsets (regions) gained by eq. 3 (Figs. 5, 6), we performed a Wilcoxon Signed-Rank test and a Kruskal-Wallis test to assess the significance of the difference to zero and the differences between medians, respectively. To consider false positives, we performed a False Discovery Rate adjustment of the p-values, resulting in 43.3% p < 0.05, 56.6% p > 0.05, median 0.08 in a data matrix of 9 years (eq. 1) versus 10 permafrost world zones indicated in Fig. 8. Boxplots represent 25-75% quartiles and whiskers are 1.5 interquartile ranges from the median.

Data availability
The GTN-P global mean annual ground temperature data for permafrost near the depth of zero annual amplitude