A ground temperature map of the North Atlantic permafrost region based on remote sensing and reanalysis data

Permafrost is a key element of the terrestrial cryosphere which makes mapping and monitoring of its state variables an imperative task. We present a modeling scheme based on remotely sensed land surface temperatures and reanalysis products from which mean annual ground temperatures (MAGT) can be derived at a spatial resolution of 1 km at conti-5 nental scales. The approach explicitly accounts for the uncertainty due to unknown input parameters and their spatial variability at subgrid scale by delivering a range of MAGTs for each grid cell. This is achieved by a simple equilibrium model with only few input parameters which for each grid cell allows scanning the range of possible results by running many realizations with different parameters. The approach is applied to the unglacierized 10 land areas in the North Atlantic region, an area of more than 5 million km 2 ranging from the Ural mountains in the East to the Canadian Archipelago in the West. A comparison to in-situ temperature measurements in 143 boreholes suggests a model accuracy better than 2.5 ◦ C, with 139 considered boreholes within this margin. The statistical approach with a large number of realizations facilitates estimating the probability of permafrost occurrence 15 within a grid cell so that each grid cell can be classiﬁed as continuous, discontinuous and sporadic permafrost. At its southern margin in Scandinavia and Russia, the transition zone between permafrost and permafrost-free areas extends over several hundred km width with gradually decreasing permafrost probabilities. The study exempliﬁes the unexploited potential of remotely sensed data sets in permafrost mapping if they are employed in multi-sensor 20 multi-source data fusion approaches.

Permafrost shapes approximately a quarter of the landmass of the Northern Hemisphere (Brown et al., 1997) and is thus one of the largest elements of the terrestrial cryosphere.Permafrost occurs mainly in arctic regions which will be strongly impacted by global warming.The recent and predicted future warming of the global climate may lead to widespread thawing which can trigger climatic feedbacks from local to global scale and severely impact ecosystems, infrastructure and communities in the Arctic.Thawing organic-rich permafrost is projected to spark substantial emissions of the greenhouse gases CO 2 , CH 4 and N 2 O (Walter et al., 2006;Schuur et al., 2008;Elberling et al., 2010Elberling et al., , 2013) ) which may be relevant for future climate projections and hereof derived mitigation strategies (e.g.Schaefer et al., 2011;Schneider von Deimling et al., 2012).These processes are poorly represented in the General Circulation Models (GCMs) used for the Intergovernmental Panel on Climate Change report (IPCC, 2013), but considerable efforts are dedicated to better capturing the effects of permafrost thaw in future global assessments.Due to the large thermal inertia of the frozen ground and associated long time periods required for thawing, future projections critically depend on the knowledge of current thermal ground conditions.However, the more than 15 years old, but still widely used global permafrost map from the International Permafrost Association (Brown et al., 1997) features a very coarse scale and lacks quantitative information on permafrost state variables, in particular of ground temperatures.In order to improve such shortcomings, Gruber (2012) derived a global high-resolution data set of permafrost probabilities based on downscaled air temperatures from reanalysis data.
Permafrost is a thermally defined subsurface phenomenon.While satellite sensors can map surface indicators of permafrost presence, such as landforms or vegetation types, remote sensing technologies can not directly measure its physical state variables, in particular ground temperature (Westermann et al., 2015a).Therefore, monitoring and mapping of the ground thermal state is restricted to either direct point observations or coarse-scale modeling using atmospheric circulation models.However, there exist a variety of satellite data sets which can be employed in numerical permafrost models to compute the thermal state of the ground, in particular Land Surface Temperature (LST) and Snow Water Equivalent (SWE).In a first proof-of-concept study, Langer et al. (2013) demonstrated the potential of such concepts by forcing a transient ground thermal model at a spatial scale of 1 km by time series of MODIS LST and downscaled GlobSnow SWE (Takala et al., 2011).While the agreement with measured in-situ data of ground temperatures and thaw depth is striking, the study also highlights the considerable sensitivity of physically-based approaches to model parameters, which are generally not known for very large spatial domains, e.g. an entire continent.On the other hand, field studies suggest a considerable variability of ground temperatures at spatial scales much smaller than 1 km (e.g.Schmid et al., 2012;Gisnås et al., 2014), which questions the validity of a single model scenario per grid cell.
Multi-scenario-runs with e.g.different snow depths and ground thermal properties (similar to tiling approaches employed in GCMs) would be an elegant way to represent the small-scale variability in a statistical way, but this implies a significant increase in the computational demands if all combinations of model parameters should be explored.In such a context, simple schemes with few parameters, such as the modified Kudryavtsev (Sazonova and Romanovsky, 2003) or the TTOP approach (Smith and Riseborough, 1996) can offer advantages over more sophisticated models.They have been developed to a state of operational maturity and are employed for large-scale mapping of permafrost properties with a wide variety of input data sets (e.g.Luo et al., 2014;Panda et al., 2014).
In this study, we employ a simple semi-empirical equilibrium model (CryoGrid 1, based on the TTOP approach, Gisnås et al., 2013) with a small number of model parameters in a statistical framework which is computationally efficient enough to be applied for large spatial domains.We demonstrate high-resolution statistical mapping of ground temperatures based on a combination of remotely sensed land surface temperatures and reanalysis data.The model results are benchmarked against measurements of ground temperatures in boreholes from the "Thermal State of Permafrost" (TSP) data base (Romanovsky et al., 2010a;International Permafrost Association, 2010).The results suggest that mapping of ground temperatures at continental scales is feasible, a variable highly demanded by the scientific community, as well as policymakers and industries in affected countries.

Study region and in-situ data
The study focuses on the permafrost areas bordering the North Atlantic Ocean which comprise a large gradient of ground temperatures, from permafrost-free areas in Scandinavia to some of the coldest permafrost on Earth in N Greenland and Ellesmere Island, Canada.
In detail, the region ranges from the Ural mountains, Novaja Zemlja and Franz Josef Land in the East over Scandinavia, Svalbard, Iceland and Greenland to the Eastern parts of the Canadian Archipelago.The area contains a wide range of permafrost environments, from high-relief mountain ranges to wetlands and subarctic mires underlain by permafrost, which makes it a well-suited test region for a modeling scheme.In-situ monitoring programs of ground temperatures in boreholes (Christiansen et al., 2010;Romanovsky et al., 2010b;Smith et al., 2010) have created a excellent record to which the model results can be compared.In total, we employ 143 borehole temperatures in N America, the Nordic region and W Russia (see Fig. 1) from the "Thermal State of Permafrost (TSP) Snapshot Borehole Inventory" (International Permafrost Association, 2010).In addition, regional-scale modeling studies have compiled maps of ground temperatures and permafrost distribution at spatial scales of 1 km (e.g.Gisnås et al., 2013;Westermann et al., 2013Westermann et al., , 2015b)).The study is based on data sets from the ten-year period 2003 to 2012, in which most of the validation data sets have been compiled.

MODIS LST
The "Moderate Resolution Imaging Spectroradiometer" (MODIS) operationally delivers remotely sensed land surface temperatures on global scale at a spatial resolution of 1 km.
The MODIS sensor is carried on board of two satellites of NASA's "Earth Observing System", Terra (launched in 2000) and Aqua (launched in 2002).As such, a time series of more than ten years is available for the combined Terra/Aqua data set.For this study, we make use of the level 3 products MOD11A1 (Terra) and MYD11A1 (Aqua) in the version 5.They contain a day and a night value for LST each, so that four values per day are available, thus in principle capturing the diurnal temperature cycle.Cloudy regions are automatically detected and removed by the MODIS cloud mask (Frey et al., 2008).In practice, frequent cloudiness leads to a drastic reduction of the data density with prolonged periods without measurements.To cover the study region, in total 16 MODISL LST tiles with an area of 1200 km×1200 km each were employed (h15v01, h15v02, h15v03, h16v00, h16v01, h16v02, h17v00, h17v01, h17v02, h18v00, h18v01, h18v02, h18v03, h19v01, h19v01, h19v02, h20v02).For each of the tiles, a ten year period from 2003 to 2012 was evaluated, so that more than 116 000 individual files were processed.
While the target accuracy of individual LST measurements is 1 K (Wan et al., 2002), several studies from arctic region suggest a significantly reduced accuracy both for individual measurements and hereof derived long-term LST averages (Westermann et al., 2011(Westermann et al., , 2012;;Østby et al., 2014).This is in particular attributed to prolonged cloudy periods with systematically different average surface temperatures that are not captured by the satellite measurements.Furthermore, detection of clouds is imperfect in particular during polar night conditions (Liu et al., 2004), so that cloud top temperatures are contained in the MODIS LST time series.Valdiation studies have shown that these effects can lead to a systematic coldbias of up to 3 K in seasonal averages (Westermann et al., 2012;Østby et al., 2014).It is therefore questionable if MODSIS LST alone can be used to model ground temperatures.To overcome this problem, we synthesize a composite of MODIS LST and ERA-interim reanalysis data as input for the ground thermal model CryoGrid 1 (see Sect. 2.4).
As exemplified by Østby et al. (2014) for the MODIS LST tile h18v00, the land mask employed in the MODIS LST production algorithm is shifted by about 5 km for the tiles N of 80 • N, so that land surface temperatures are produced for sea areas, while they are missing for some land areas.In the study region, this mainly affects the N parts of Svalbard, Greenland and Ellesmere Island.Reanalysis products provide global datasets of meteorological variables based on atmospheric models in which a range of observations, such as meteorological surface observations, remotely sensed sea surface temperatures, or vertical temperature and humidity profiles from radio soundings, are assimilated.In this study, we make use of air temperatures of the ERA-interim reanalysis (Simmons et al., 2007;Dee et al., 2011), which is available from 1979 until now at a spatial resolution of 0.75 • × 0.75 • .ERA-interim is the state-ofthe-art reanalysis product from the European Centre for Medium-Range Weather Forecast (ECMWF), and it has been successfully employed in previous permafrost modeling studies (Fiddes et al., 2015).While the atmospheric model can provide a gap-free time series with four values per day (at 00:00, 06:00, 12:00, and 18:00 UTC), the coarse spatial resolution requires a down-scaling scheme to correct for the difference between the true altitude of a point and the coarse-scale ERA orography, which in the study region can be as much as 1500 m in mountainous terrain.To obtain a downscaled air temperature for a point in space and time, Fiddes and Gruber (2014) propose a four-dimensional interpolation between different reanalysis grid cells, pressure levels and time steps.Since only long-term freezing and thawing degree days are required as input for CryoGrid 1 (Sect.2.7), a simplified version of this scheme (similar to Gruber, 2012) is employed: 1.An average atmospheric lapse rate is derived for each day of year (DOY) and MODIS grid cell, based on ERA-interim surface fields and the 700 mbar pressure level, which corresponds to an altitude of approximately 3000 m, slightly higher than the highest non-glacierized areas within the study region.For each DOY, all available ERA-interim data from 2003 to 2012 are employed.
2. The downscaled air temperature is computed from the lapse rate and the difference between the ERA orography (interpolated to the center position of a MODIS grid cell) and the "true" altitude.The "true" altitude of the 1 km MODIS LST grid cell is derived by interpolation of the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) at 30 arcsec resolution (Danielson and Gesch, 2011) to the center position of each MODIS 1 km grid cell.This data set is specifically compiled for continental-scale applications and includes the best available global elevation data.This procedure ensures that the spatial and seasonal variability of the lapse rate is accounted for so that e.g.inversions of the air temperature during winter in continental regions can be reproduced.

Data fusion of MODIS LST and ERA-interim
Due to frequent cloudiness, the MODIS LST time series contains a large number of gaps and only clear-sky LST averages can be derived from remote sensing products alone.To eliminate or at least moderate the effect of clouds, the MODIS LST time series is merged with the gap-free time series of downscaled air temperatures from the ERA reanalysis (see Sect. 2.3).For this purpose, the acquisition time of each MODIS day and night LST measurement (which varies by several hours from day to day) is converted to UTC, and the ERA value closest in time to each LST measurement subsequently removed.From the resulting time series, daily averages are computed from which freezing and thawing degrees are derived.Depending on the cloudiness, the fraction of MODIS LST measurements in the composite MODIS/ERA data set is variable (Fig. 1).In areas with frequent cloudiness, such as Iceland, less than a quarter of the data points are derived from MODIS LST, while the fraction of MODIS LST is more than 50 % in other areas.

ERA snowfall
The average yearly snowfall from 2003 to 2012 was determined for each grid cell by interpolation of the snowfall surface field coarse-scale ERA-interim reanalysis data.These ), it captures large-scale precipitation patterns, such as the high snowfall in the S Norwegian mountains close to the Atlantic Ocean.On the other hand, it cannot account for small-scale variations of winter precipitation.However, since only the range of possible n f factors is assigned in a statistical modeling framework with many different model realizations (Sect.2.7), utilizing large-scale precipitation patterns may indeed be adequate.

MODIS land cover
The land cover has strong implications for the small-scale distribution of the snow cover: in areas with low vegetation or bare ground, redistribution of snow due to wind drift can lead to a strong spatial variability of snow depths and ground temperatures (e.g.Gisnås et al., 2014).In contrast, a more uniform snow depth is expected in forested areas.To represent such effects in the modeling (Sect.2.7), the MODIS land cover product MCD12Q1 for the year 2006 at a spatial resolution of 500 m is employed.Due to the poor performance of the landcover classification in some parts of the study region (see Sect. 4.1), the 500 m grid cells are aggregated to the 1000 m MODIS LST grid and the land cover classes merged to two units: high vegetation/forest ("International Geosphere-Biosphere Programme" -IGBP classes 1 to 9, 11) and bare ground/low vegetation (IGBP classes 10, 15, 16).The classes 12 to 14 (cropland/urban) occur only sparsely in permafrost areas and are treated like the high vegetation class.

Statistical modeling with CryoGrid 1
CryoGrid 1 is an equilibrium model for ground temperatures (Gisnås et al., 2013) based on the TTOP approach (e.g.Smith and Riseborough, 1996).It computes mean annual ground temperature (MAGT) based on thawing and freezing degree days (TDD s and FDD s ) of surface temperatures, according to where τ is the period for which TDD and FDD are evaluated, while r k , n f and n t are semiempirical parameters which aim to capture a variety of key processes in a single variable.
The winter n f factor relates the freezing degree days at the surface (here derived from MODIS/ERA) to the freezing degree days at the ground surface as (2) It thus captures the insulating effect of the snow cover which is a result of snow depth, the thermal properties of the snow cover, and the ground thermal regime itself.For bare surfaces, this variable is unity, while n f factors as low as 0.2-0.3 are reported from in-situ measurements (e.g.Gisnås et al., 2013) for high snow cover in Scandinavia.The summer n t factor is defined in a similar way as TDD gs = n t TDD s . (3) If air temperatures are employed as surface forcing TDD s , variable n t factors are considered to account for differences in air and ground surface temperatures (e.g.Gisnås et al., 2013).
In this study, we assume remotely sensed LST in conjunction with ERA-interim to be a satisfactory representation of ground surface temperature and set n t to unity.In case of a dense canopy, where MODIS LST may rather represent top-of-canopy instead of surface temperatures, a different n t may be required, but such conditions rarely occur in permafrost areas in the study region.The active layer damping factor r k gives rise to the thermal offset between average ground surface and ground temperatures (Osterkamp and Romanovsky, 1999).It is related to the average thermal conductivities in the active layer in fully thawed and frozen states (k t , k f ) as r k = k t / k f , so that it is closely related to the water/ice content of the soil.The thermal conductivities of water and ice (k For each of the land cover units "high vegetation/forest" and "bare ground/low vegetation" (Sect.2.6), upper and lower bounds for n f and r k are defined (Table 1), oriented at physical constraints and published values from previous studies.However, the upper and lower bounds of n f and r k are also tuning parameters which can be adjusted to achieve the best possible match between model and observations (Sect.3.1) for each land cover unit.We emphasize that remotely sensed data sets on precipitation or snow depth from which the n f factors could be estimated (Gisnås et al., 2013) are not available at 1 km scale and the interval of possible values for n f was chosen according to large-scale ERA-interim snowfall data (Sect.2.5).For the "bare ground/low vegetation" class, the upper bound for n f is set to 1, since bare-blown spots without insulating snow cover can occur.The lower bound decreases with increasing snow fall, i.e. for areas with high snow fall, the spread of possible snow depths and thus n f -values is larger than for areas with low snowfall.For the mountain areas in Norway, Gisnås et al. (2013) assumed a minimum n f of 0.3, while the minimum lower bound of n f is around 0.1 for Norway in this study.However, Gisnås et al. (2013) employed grid-cell averages of snow depths to calculate n f , while in this study the lower bound of n f aims to represent the areas with maximum snow depths within each grid cell (e.g.snow drifts).The parameterization of n f vs. snow depth applied by Gisnås et al. (2014) to represent such small-scale variability of snow depths yields n f values between 0.1 and 0.2 for snow depths of several meters, which regularly occur in the Norwegian mountain areas, in agreement with the choice in this study.For the "high vegetation/forest" class, both upper and lower bound for the winter n f -factor decrease linearly with snow fall, while the difference between minimum and maximum n f is held constant.The r k values are chosen close to 1 for both classes, which implies that large changes in the thermal conductivities do not occur.Areas with potentially large thermal offsets (i.e.low r k values), such as wetlands or extensive block field areas (Gruber and Haeberli, 2007), are clearly not represented by this choice.However, with the currently available land cover products, it is not possible to reliably detect such areas (see Sect. 4.1.3),so that they cannot be accounted for by the employed scheme.Furthermore, variations of altitude and exposition within grid cells are an additional source of spatial variability of ground temperatures which is not accounted for For each grid cell, CryoGrid 1 is run for 30 values of n f and r k (in equally spaced intervals from the lower to the upper limit), yielding a total 900 different values for the mean annual ground temperature (MAGT).The statistics of these values is a representation of both the small-scale (i.e.< 1 km) variability of model parameters and the model uncertainty due to the generally unknown "true" parameter distribution of a grid cell.As final result, we calculate the average of all 900 realizations of MAGT, as well as the standard deviation, and the maximum and minimum MAGT.

Comparison to in-situ data
We compare the mean of all model realizations to in-situ measurements of ground temperatures from the "IPA-IPY Thermal State of Permafrost (TSP) Snapshot Borehole Inventory" (International Permafrost Association, 2010) documented in more detail for N America (Smith et al., 2010), Russia (Romanovsky et al., 2010b), and the Nordic areas (Christiansen et al., 2010).While this data basis provides only a single value for ground temperatures measured at different depths and different points in time, it is the most extensive compilation of permafrost temperatures available for model validation (see Sect. 4.1.1).For comparison, we generally select the grid cell closest to the borehole location.For boreholes located close to a larger water body, a nearby grid cell located at least 1 km from the shore is employed to avoid contamination of the MODIS LST record with surface temperatures from the water body.For the boreholes north of 80 • N (five boreholes near Alert, Ellesmere Island, Canada), no MODIS LST measurements are available due to the erroneous land mask (Sect.2.2).We therefore use the closest grid cell featuring MODIS LST measurements E of Alert which is located on land.The results of the comparison are displayed in Fig. 2, with 12 boreholes located in N America (Canadian Archipelago), 69 in  (Scandinavia, Greenland, Svalbard, Iceland) and 62 in Russia.For each borehole, the error bars represent one SD of the results of all model realizations for a grid cell which depends on the ranges of n f and r k .For the large majority of the boreholes, the agreement between measured and modeled ground temperatures is better than 2 • C, with 67 (47 %) contained within one, 108 (75 %) within two and 128 (90 %) within three SDs of all model realizations from the mean (see Supplement).
A histogram of deviations between measured and the mean modeled ground temperatures is displayed in Fig. 3.It roughly follows a Gaussian distribution with width of 1.2 • C.
For 60 % of the boreholes, the mean of all model realizations is within 1 • C of the measurement, while the agreement is better than 2 • C for 93 % and better than 2.5 • C for 97.5 % of the boreholes.From the comparison with in-situ measurements, we conclude that the accuracy of the modeled ground temperatures is on the order of 2 to 2.5 • C for the study region.There is no significant regional bias for Russia and for the Nordic areas, while there is a slight cold-bias for the available boreholes in N America.However, the latter cannot be fully secured due to the comparatively weak data basis of only 12 boreholes, all of which are assigned large model standard deviations of 1.9 to 2. between minimum and maximum corresponds well to the range of measured AGST (−4.5 to +0.5 • C).Finally, the spatial variability of ground temperatures has been modeled for the high-Arctic Zackenberg site in NE Greenland, taking into account the spatial variability of the snow cover, ground surface and ground properties (Westermann et al., 2015b).For the period 2002-2012, modeled annual average ground temperatures at 1 m depth range from −9.0 to −3.5 • C, while the model realizations of CryoGrid 1 yield a temperature range from −11.0 to −2.9 • C. While not a validation in a strict sense, the satisfactory to good agreement for these few examples suggests that the ensemble of all model realizations can represent the small-scale spatial variability of ground temperatures, at least for the land cover class "bare ground/low vegetation".Similar field data sets or modeling studies do not exist for the "high vegetation/forest" class, so that is not possible to benchmark the modeled temperature range at this point.

Modeled ground temperatures in the N Atlantic permafrost region
The model approach facilitates large-scale mapping of the ground thermal regime on a continental scale.As shown in the previous section, the comparison to in-situ measurements suggests that the mean of all model realizations reproduces the ground thermal regime within an accuracy of 2 to 2.5 • C. regions of the Pechora Sea E of the Kanin peninsula.In the high elevations of the Ural mountains, subzero temperatures occur as far south as 60 • N. In the lowlands E of the Ural mountains, negative modeled ground temperatures reach significantly farther south compared to the W side.In N Scandinavia and Russia, large areas with modeled ground temperatures just above 0 • C exist, for which it is not possible to establish permafrost presence or absence considering the accuracy of 2 to 2.5 • C accuracy limit (see Sect. 3.3).The modeled permafrost distribution in the N Atlantic region is in good qualitative agreement with the "classic" permafrost map of the IPA based on field evidence (Brown et al., 1997), including prominent features like the asymmetry of the permafrost extent E and W of the Urals.
The standard deviation of all model realizations is displayed in Fig. 5.The highest values of around 3 • C are reached in N and E Greenland, as well as on Baffin Island, areas with a large number of freezing degree days that are classified as "bare ground/low vegetation" and hence are assigned a wide range of n f -values.The large spread therefore stems from the strong variability of winter ground surface temperatures reflected in the term n f F DD in Eq. 1. Permafrost areas further south, e.g. in Scandinavia, feature lower numbers of F DD and hence a smaller standard deviation, generally between 1 and 2 • C. Areas classified as "high vegetation/forest" have standard deviations of 1 • C or less, mainly due to the smaller range of winter n f -factors, which reflect the more homogeneous snow cover in such areas.
At 1 km spatial resolution, it is meaningful to "zoom" into subregions for a detailed assessment of the ground thermal regime.Figure 6  In Iceland, subzero ground temperatures are restricted to the interior, ranging from the central Sprengisandur between Hofsjökull and Vatnajökull north to the mountain ranges of the Tröllaskagi peninsula, where active rock glaciers suggest a periglacial environment (Lilleøren et al., 2013).The modeled permafrost extent is in good agreement with previous regional estimates based on mean annual air temperatures (Etzelmüller et al., 2007;Farbrot et al., 2007).
On Svalbard, modeled ground temperatures in the ice-free areas range from −2 to −3 • C at the W coast to −5 to −6 • C in the ice-free parts of Nordaustlandet.In Scandinavia, the southernmost areas with subzero modeled ground temperatures are located at high elevations in the S Norwegian mountains.Within the accuracy, the modeled ground temperatures agree well with the spatially distributed modeling of Westermann et al. (2013), showing rather warm permafrost with ground temperatures generally warmer than −3 • C. In general, more areas are mapped with subzero temperatures compared to the previous study (especially in the Hardangervidda area), but the differences in the absolute modeled temperatures are mostly less than 1.5 • C (see also Sect.3.3), and patchy permafrost is known to exist in these regions (Etzelmüller et al., 2003).In N Scandinavia, permafrost is modeled in higher elevations of the Swedish and Norwegian mountains.In the latter, the pattern is again in agreement with the regional assessments of Farbrot et al. (2013) and Gisnås et al. (2013) who employed CryoGrid 1 with gridded air temperature and snow depth products.
In the rolling plains of the Finnmarksvidda, the model shows temperatures just above 0 • C, while sporadic permafrost is found in high-lying fell areas and palsa mires (Farbrot et al., 2013).In this area, the model approach clearly fails to reproduce the permafrost patterns.Most likely, a main reason is the employed MODIS landcover product, in which almost the entire area is uniformly classified as "open shrubland" (and thus as "high vegetation/forest" in our approach), although it is characterized by a complex pattern of fell areas, tree-less mires and areas covered by mountain birch forest.If the model is applied to the area with the parameter set for the "bare ground/low vegetation class" (Ta- ), subzero ground temperatures are modeled for a large part of the area.Therefore, the model approach would indeed yield a pattern of permafrost and permafrost-free areas if a landcover product capable of distinguishing mountain birch forest from bare fell and tundra areas is employed.The same issue applies to modeled ground temperatures on the Kola peninsula in Russia, where ground temperatures around 0 • C are modeled for the N and E parts.

Permafrost probability maps from statistical modeling
With an estimated accuracy of 2 to 2.5 • C, it is impossible to delineate exact boundaries for the permafrost extent of the study region based on the mean of all model realizations.Around the thaw threshold, the statistical model approach features model realizations below and above 0 • C, depending on the input parameters n f and r k .However, permafrost and permafrost-free areas coexist in this transition zone also in nature.This long-known fact (e.g.Baranov, 1959) gave rise to the classification system "continuous" (> 90 % of an area underlain by permafrost), "discontinuous" (50-90 %), and sporadic (< 50 %) permafrost, as defined by the IPA long before the advent of sophisticated modeling techniques (Brown et al., 1997).
In contrast to model approaches with only a single model realization for each grid cell, the ensemble of realizations of CryoGrid 1 facilitates deriving a similar permafrost zonation from the modeling.Figure 8 displays the fraction of model realizations with subzero ground temperatures in the same classes as the IPA permafrost map for Scandinavia and W Russia.The potential permafrost extent is larger than in Figs. 4 and 7, with the most significant differences in lowland areas in Scandinavia and Russia.The interior of N Scandinavia (for which the mean of all model realizations is above 0 • C, Sect.3.2) is classified as sporadic permafrost in this approach, which is in good agreement with observed field evidence for this area.In Russia, sporadic to discontinuous permafrost is modeled from the Kola and Kanin peninsulas to the mouth of the Pechora River, east of which the permafrost is classified as continuous.Also in this probability map, the individual permafrost zones extend significantly further south on the east side of the Ural Mountains than on the west side.This asymmetric permafrost distribution around the Ural mountains is a prominent feature in the IPA permafrost map (Brown et al., 1997), which the presented satellite-based modeling approach is capable of reproducing.In contrast, the global map of Gruber (2012), which also delivers comparable probabilities of permafrost occurrence, shows a more or less symmetric permafrost extent around the Urals.Further pronounced differences occur on the eastern tip of the Kola peninsula where the map by Gruber (2012) shows no permafrost, other than the IPA permafrost map and the satellite-based modeling approach.Compared to the IPA map, the transition zone from permafrost-free areas to the continuous permafrost zone in NW Russia is broader in the satellite-based map, which is an indication that some model scenarios (i.e.combinations of model parameters) do not occur in nature.A similar broadening of the transition zones is visible in the probability map of Gruber (2012).

The equilibrium model
The TTOP-approach is one of the older modeling schemes conceived to obtain ground temperatures and thus permafrost occurrence from near-surface meteorological variables (e.g.Smith andRiseborough, 1996, 2002).It delivers a "mean annual ground temperature" at the depth corresponding to the top of the permanently frozen ground under the strong assumption that the ground thermal regime is in thermal equilibrium with the applied surface forcing data (Appendix A).Since this assumption is violated to a certain degree for realworld conditions, a comparison to in-situ measurements, as conducted in Sect.3.1, is not a straight-forward task.For a typical geothermal gradient of 0.025 • Cm −1 , measured M AGT s would be less than 1 • C warmer than modeled M AGT s for borehole depths of up to 40 m.
On the other hand, ground temperatures have been warming by a similar magnitude even at depths of 10 to 20 m in the model period (Romanovsky et al., 2010a), so that recorded borehole temperatures may represent colder climate conditions compared to the model period.With the equilibrium approach, it is thus principally impossible to exactly reproduce ground temperatures that can be measured in the field.However, we note that the combined uncertainty due to measurement depth and transient effects is not necessarily systematic, and presumably well below the estimated accuracy of the modeling scheme of 2 K for most boreholes.

Land surface temperature
In this study, land surface temperatures are derived from satellite measurements during clear-sky conditions, while air temperatures from reanalysis products are employed for cloud-covered periods when satellite measurements are unavailable.This approach combines the strengths of both products -the capacity of satellite sensors to provide actual measurements of the footprint area at a high spatial resolution, and the dense, gap-free time series of the reanalysis product.On the other hand, it strongly reduce the uncertainty associated with both input data sets when used separately: even when downscaled with altitude, reanalysis products only provide a large-scale temperature field which does not account for the heterogeneity of the surface temperature caused by spatially variable surface conditions.Employing only MODIS LST and thus clear-sky LST to calculate long-term averages can result in a serious bias since cloudy periods with a significantly different surface temperature regime are not contained (e.g Westermann et al., 2012;Østby et al., 2014).
However, these validation studies also concluded that a significant number of highly erroneous measurements due to undetected clouds is contained in the MODIS LST time series.These measurements are not removed in the synthesized MODIS/ERA time series, which could in principle introduce a bias in the freezing and thawing degree days.Since including ERA reanalysis data strongly increases the total number of values in the time series, these biased values will attain a smaller weight in the FDD s and TDD s calculation compared to employing MODIS LST only, so that their influence is moderated.Nevertheless, it should be investigated if highly erroneous MODIS LST can be detected by a consistency check In this study, the ERA-interim reanalysis is employed to complement satellite-based LST measurements.Since a large number of operational meteorological observations are assimilated in the reanalysis (Dee et al., 2011), independent validation data sets are rare especially in Arctic regions.Over the central Arctic ocean during summer, Jakobson et al. (2012) found near-surface temperatures to be warm-biased by up to 2 • C, but the performance improved at higher altitudes.For Ireland, the performance for near-surface temperatures was excellent with a root mean square error of generally less than 0.5 • C, but a slight warm bias during the winter months (Mooney et al., 2011).While these examples give a hint on the expected accuracy, the performance should be investigated further, especially in regions with pronounced topography.The coarse-scale temperature field of the ERA-interim reanalysis is downscaled using high-resolution topography and a seasonal lapse-rate calculated from the surface and the 700 mbar levels.The similar downscaling scheme presented by Fiddes and Gruber (2014), which is based on a 4-D interpolation between all available pressure levels and time steps, did not provide satisfactory results for grid cells with elevations below the ERA orography.This case is common e.g. at the coasts of Svalbard and Greenland, where the elevation difference can be more than 1000 m.In these cases, the interpolation relies on the "near-surface lapse rate", i.e. the modeled temperature difference between the surface and the first pressure level, which is subsequently extrapolated to altitudes far below the surface level.In particular if an ERA grid cell is snow-covered, this near-surface lapse-rate can feature strong temperature inversions, so that extrapolation to lower altitudes leads to a strong cold-bias of the downscaled temperature.To moderate the effect of near-surface stratifications, we employ an "average lapse rate" for the entire ERA air column between the surface and the 700 mbar pressure level located at approximately the elevation of the highest non-glacierized areas in the study region.We emphasize that the procedure cannot account for many regional and local climate and weather conditions, such as temperature inversions in valley systems.Downscaling of reanalysis data using e.g. the Weather Research and Forecasting (WRF) Model, though computationally expensive, may be a way to overcome such difficulties (Aas et al., 2015).

Land cover
For each 1 km grid cell, the range of possible values for n f and r k is assigned based on the MODIS landcover product and large-scale snowfall data sets.Hereby, the main goal is to distinguish areas with bare ground or low vegetation, where strong redistribution of snow due to wind drift can occur, from areas with trees and high vegetation where the snow cover is more uniform.At least in some areas, the performance of the MODIS landcover product in distinguishing between the two classes seems to be poor: in the rolling plains of Finnmarksvidda, N Norway, mountain birch forest alternates with bare fell areas, while the entire area is classified as "open shrubland" in the MODIS landcover and thus as "high vegetation/forest" in the modeling.In the studied region, this problem is restricted to N Scandinavia and Russia, where permafrost occurs close to and within forested areas.In regions with more extensive forest cover in permafrost areas, this issue should be investigated in more detail.The new global land cover classification delivered by the "Climate Change Initiative" of the "European Space Agency" (ESA CCI, Hollmann et al., 2013) may constitute an improved product for permafrost modeling, but it's performance in distinguishing areas with low trees and shrubs from bare tundra should be investigated in more detail.However, with this product, it may be feasible to distinguish additional classes which differ in their parameter ranges employed in the modeling (Table 1), at least in regional applications.Ideally, a pan-arctic land cover classification should be compiled which is specifically tailored for permafrost applications.Such a product should not only deliver information on the vegetation and surface state, but also on subsurface properties, which could e.g.be employed to better estimate the r k factors in the presented model scheme.The presented ground temperature map is compiled without employing fine-scale data sets on precipitation or snow depth.Instead, only plausible ranges of values for n f are assigned, which account for both total winter precipitation (i.e.spatially averaged snow depth) and wind redistribution of snow.In a recent study for Svalbard and Norway, Gisnås et al. (2014) demonstrated that annual average ground surface temperatures within areas of approximately one modeling grid cell (1 km) vary by up to 5 • C in areas with pronounced wind drift, while a borehole can only provide one point value sampled from the distribution of ground temperatures.The modeling approach explicitly accounts for this spatial variability by computing ground temperatures for a whole range of n f value.We find that the agreement with borehole temperatures, despite using large-scale data sets of winter precipitation, is around ±2 • C and thus similar to natural spatial variability of ground temperatures.

Deterministic vs. statistical modeling of ground temperatures
In the version employed for this study (Eq.1), the TTOP approach accounts for the insulating effect of the snow cover, as well as the thermal offset between mean annual ground surface and ground temperatures.The complex physical processes which give rise to these effects, are lumped into one variable each (n f and r k ) which must be adjusted to achieve adequate results.While both n f and r k are in principle empirical constants, they are subject to physical constraints which confines them to a certain range.Other than deterministic modeling schemes which seek to employ a single set of model parameters for each grid cell, we use the cumulative information content of many realizations, i.e. combinations of different parameters selected from the physically possible range (Sect.4.1.1).Due to its simplicity and the small number of input parameters, the TTOP-approach is computationally efficient enough to compute a large number of realizations.The scheme can be expected to deliver adequate results if the uncertainty due the simplified model (Appendix A) is smaller than the uncertainty due to the unknown model parameters.We note that the model approach can not separate the spatial variability of the input parameters within a model grid cell from the uncertainty due to generally unknown values of these parameters.Ideally, it delivers a range of possible MAGT's in which the true range of MAGTs (as documented by measurements, e.g.Gisnås et al., 2014) is contained.While the agreement between modeled and measured ranges is excellent for the few in-situ data sets on small-scale spatial variability of ground temperatures (Sect.3.1), more work and in particular improved validation data sets are needed for a robust assessment of modeled and true ranges of MAGTs.
More sophisticated modeling schemes (e.g.Jafarov et al., 2012;Zhang et al., 2012;Fiddes et al., 2015;Westermann et al., 2013) feature significantly more model parameters which are required to achieve a realistic description of physical processes.While these parameters can generally be estimated to a certain level of confidence in regional-scale applications (e.g.Fiddes et al., 2015), their variation and spatial distribution is hard to access on continental to global scale.It is therefore questionable whether the performance of such computationally expensive schemes is better than the simple CryoGrid 1 model.Furthermore, the significant number of input parameters and the computational costs make sensitivity studies with many realizations unpractical, at least for many grid cells.In the CryoGrid 1 based scheme, however, a sensitivity is performed for each grid cell so that confidence limits for the results can be given.On the point scale, Langer et al. (2013) performed a sensitivity analysis for a transient permafrost model.They found modeled ground temperatures to vary by up to 4 • C if the input parameters related to snow and soil properties were allowed to vary within physically plausible limits.This suggests that on a continental scale the performance of more sophisticated schemes could be on the order of the much simpler TTOP approach.

Towards a global high-resolution ground temperature map
The "classic" permafrost map compiled by the IPA more than 15 years ago (Brown et al., 1997) is based on available field evidence for permafrost occurrence and properties.Since then, new global data sets have become available which so far have only rarely been used for permafrost mapping on large spatial scales.An exception is the global approach by Gruber (2012) who derived the probability of permafrost occurrence at 1 km scale from In this study, we have demonstrated the continental-scale application of a statistical model approach with global satellite data sets and a reanalysis product as input data.In addition to probabilities of permafrost occurrence (as in Gruber, 2012), the approach facilitates mapping of ground temperatures and thus quantifying the thermal state of the permafrost.
The model was successfully applied to an area of more than 5 million km 2 with strong spatial differences in ground temperatures and a wide variety of permafrost landscapes which suggests that the approach is scalable to include the entire contiguous permafrost extent on the Northern Hemisphere, an area of approximately 23 million km 2 (Brown et al., 1997).
Hereby, it may become necessary to refine the parameterizations for the ranges of n f and r k (Table 1), possibly by introducing functional dependences on other environmental variables.A main uncertainty is the performance in forested regions which could not be sufficiently investigated in the N Atlantic study region.At least in case of dense forests or strong offnadir angles remotely sensed LST may represent top-of-canopy temperatures instead of skin temperatures at the ground or snow surface, as required by CryoGrid 1.Thus, it should be investigated whether this effect leads to a reduced accuracy or even a systematic bias of modeled ground temperatures.
In addition to a global mapping, the presented method is also feasible for regional-scale applications, especially if improved data sets on e.g.surface cover, air temperature or snow depth are available.An example is the ground temperature map of Norway compiled by Gisnås et al. (2013) through application of CryoGrid 1 with gridded air temperature data sets, which could be improved both by integrating remotely sensed LST and by accounting for the subgrid variability of the snow cover in a statistical framework.However, we emphasize the semi-empirical nature of the model approach so that it should only be applied in areas with sufficient in-situ data sets for validation.Furthermore, in mountain areas with strong topographic variations, satellite-based LST measurements at 1 km scale cannot sufficiently capture variations of altitude and exposition, so that the modeling scheme at best can be expected to deliver a first-order approximation of the permafrost distribution.In this study, we present a modeling approach for ground temperatures using remote sensing data and reanalysis products as input data that could be operationally applied on continental to global scale at a scale of 1 km.Based on the equilibrium permafrost model Cryo-Grid 1, remotely sensed land surface temperatures and land cover are employed in conjunction with air temperatures and snowfall from the ERA-interim reanalysis.The coarse spatial resolution of reanalysis data is compensated by a statistical modeling framework which scans a range of possible model parameters and thus facilitates estimating the uncertainty.The approach is applied to the permafrost areas around the North Atlantic, and area with a wide variety of permafrost environments and strong differences in ground temperatures.
-The mean of all model realizations is within 2.5 • C of in-situ measurements for 139 of 143 boreholes contained in the "Thermal State of Permafrost" network.The differences between modeled and measured ground temperatures can be approximated by a Gaussian distribution with width 1.2 • C centered around zero.
-For tundra areas, the modeled spread of ground temperatures within the ensemble of model realizations is roughly in agreement with measurements of the spatial variability of the ground thermal regime caused by variable snow depths due to wind drift of snow.
-The ensemble approach with many model realizations facilitates assigning fractions of permafrost or permafrost probabilities for each grid cell.It is therefore possible to revamp the "classic" permafrost zonation system with sporadic, discontinuous and continuous permafrost at a scale of 1 km.
-Within the accuracy established by comparison to in-situ measurements, the modeled permafrost extent agrees well with field observations and previous studies in all investigated regions.The study demonstrates the potential of remotely sensed data sets for permafrost modeling if they are employed in multi-source data fusion approaches.It constitutes an important step towards satellite-based mapping of the ground thermal regime at unprecedented spatial resolutions, but uncertainties related to satellite-based land cover maps and associated ranges of model parameters must be investigated and resolved prior to global application.Such novel ground temperature products bear significant potential for permafrost science, e.g. by improving the validation of climate model results.Furthermore, they can be of high value for societies and industries in permafrost-affected countries, e.g. for the planning of infrastructure.

Appendix A: Alternative derivation of the TTOP equation
A derivation of the TTOP equation (Eq. 1) is found in Romanovsky and Osterkamp (1995).Here, we present an alternative derivation which allows insight in the simplifications that are applied in the TTOP model.
We consider the permafrost case, M AGT ≤ 0 • C.During the summer, the active layer thaws and features a thermal conductivity k t .At the bottom of the active layer (i.e. at the top of the permafrost) at depth d, we assume a constant temperature of 0 • C. Furthermore, d and k t are assumed constant during the entire thaw season, which is a gross simplification of the natural processes.At each time t, the heat flux in the permafrost is then computed as a linear interpolation between the ground surface temperature T gs (t) as During winter, the active layer is frozen with thermal conductivity k f and a temperature of Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 1 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2 Modeling tools and satellite data sets Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2.3 Downscaled near-surface air temperatures from ERA-interim reanalysis

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | snowfall data are employed to determine the range of the winter n f factors for the statistical modeling (in conjunction with a remotely sensed land cover data, see below).While there are well documented biases in the precipitation fields of the reanalysis (e.g.Schmidli et al., Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2006 1) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | confine the range of physically possible values of r k between one for dry soil or rock and approximately 0.2 for pure water.The statistical modeling is based on the assumption that FDD s and TDD s can be exactly determined from the MODIS/ERA-composite, while the two remaining variables in Eq. (1), n f and r k are unknown and only confined by physically plausible limits (note that n t = 1 is Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | assumed).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | in the presented model scheme.In mountain areas with strong topographic variations, the spatial variability of ground temperatures within a grid cell is most likely underestimated.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the Nordic area 7 • C (see Supplement) which indicate a potentially large spatial variability of ground temperatures.While the borehole data of the TSP network are well suited to validate the large-scale performance of the modeling, only very few quantitative studies on the spatial variability of ground temperatures within areas of a model grid cell exist to which the ensemble of all model realizations could be compared.Recently,Gisnås et al. (2014) presented a one-year data set of ground surface temperature measurements based on 40 to 100 temperature sensors distributed in three areas of approximately 0.5 km 2 , located in Svalbard and S Norway.Since the thermal offset between average ground surface and ground temperatures is considered small for the sites, we compare the measured distribution of average ground surface temperatures (AGST) to the ensemble of modeled ground temperatures.For Finse (60 • 34 ′ N, 7 • 32 ′ E), 30 % of the temperature sensors displayed AGSTs below 0 • C, compared to 51 % of the model relization for the MODIS grid cell comprising the study site (for ground temperatures).The measured range of AGST was −2 to +2.5 • C, while minimum Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and maximum modeled ground temperatures were −2.2 and +2.1 • C. For the Juvvasshøe site (61 • 41 ′ N, 8 • 23 ′ E), 76 % of the measured AGSTs were below 0 • C, compared to 77 % of the model realizations.The modeled temperature range (−3.8 to 1.2 • C) was slightly larger than the measured range of AGST (−2 to +1.5 • C).For Ny-Ålesund (78 • 55 ′ N, 11 • 50 ′ E), the modeled temperatures are slightly cold-biased, with measured borehole temperatures of −2.3 • C (International Permafrost Association, 2010) compared to a mean of all model realizations of −3.8 • C. Nevertheless, the modeled temperature range of 4.9 • C Figure 4 displays the resulting ground temperature map for the permafrost regions bordering the N Atlantic Ocean.The modeled ground temperatures span a wide range, from −15 • C in N Greenland and Ellesmere Island to more than +5 • C in S Scandinavia.Ground temperatures below 0 • C are modeled in the Canadian Archipelago, Greenland, Iceland, Svalbard, the Scandinavian mountains and the coastal Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | displays a ground temperature map for Greenland and Iceland, with a significantly finer resolution compared to the 25 km scale assessment ofDaanen et al. (2011).In the ice-free parts of N Greenland, modeled ground temperatures range from −10 and −15 • C, while they gradually increase southwards in the ice-free regions of NE Greenland to around −3 • C in the Scoresbysound region.Further South, only few ice-free land areas exist on the E coast, but subzero ground temperatures are modeled in the Tasilaq region at 65.5 • N. On the W coast, the permafrost extends from the outer coast line to the ice sheet in the extensive ice-free land regions S of Disko Bay, with modeled ground temperatures warmer than −5 • C. Starting just N of Nuuk, the coastal Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |areas become permafrost-free, but subzero ground temperatures are still modeled at higher elevations in S Greenland.

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | ble 1 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | with downscaled ERA air temperatures, e.g. by setting a threshold value for the difference between the two temperatures.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 4.1.4Precipitation and snow depth Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | downscaled air temperatures (obtained from atmospheric model output) using empirical probability functions.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5 Conclusions

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

TFigure 1 .Figure 3 .
Figure 1.Fraction of MODIS LST measurements in the combined MODIS/ERA surface temperature product.Borehole sites from International Permafrost Association (2010) used for validation of the modeled ground temperatures, note that a dot can represent several boreholes.Further boreholes outside the displayed regions are employed.