Icelandic permafrost dynamics since the Last Glacial Maximum e model results and geomorphological implications

throughout the late Pleistocene. This in ﬂ uenced geomorphological processes and landform generation: the early collapse of the marine-based ice sheet together with the aggradation of permafrost in these zones initiated the formation of abundant and now relict rock glaciers across coastal margins. Permafrost degraded rapidly after the Younger Dryas, with a marked impact on slope stability. Permafrost that formed during the Little Ice Age is again thawing rapidly, and an escalation in slope failure and mass-movement might be currently underway. Our study demonstrates that large regions of Iceland have been underlain by permafrost for millennia, facilitating landform development and in ﬂ uencing the stability of steep slopes. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Periglacial processes related to the dynamics of permafrost exert a strong influence on mass wasting (Hales and Roering, 2007;Matsuoka and Murton, 2008) and slope stability (Gruber and Haeberli, 2007;Krautblatter et al., 2013).They are integral for understanding landscape development over recent geological time scales (Andersen et al., 2015;Berthling and Etzelmuller, 2011;Egholm et al., 2015;Hales and Roering, 2009).Characteristic landforms in permafrost of cold mountains relate to long-term cumulative slope deformation through cohesive flow as a consequence of subsurface ice formation and preservation.Such landforms include rock glaciers derived from slope material accumulations (talus, rock and snow avalanches e "talus-derived rock glaciers") or from perennially frozen ice-cored moraines deposited in slopes, where the debris cover mantles the stagnant ice and prevents it from melting ("moraine-derived rock glaciers", (Anderson et al., 2018b;Berthling, 2011;Jones et al., 2018;Knight et al., 2019;Lilleøren et al., 2013).
In Iceland, permafrost is widespread at present in highlands and mountains over c. 800 m a.s.l. and sporadic in palsa mires in the central highlands (Czekirda et al., 2019;Etzelmüller et al., 2007;Saemundsson et al., 2012), covering around 8% of the total land area (Etzelmüller et al., 2007).Because of the regional climate conditions and high geothermal heat fluxes in Iceland, permafrost is generally thin with temperatures close to the melting point, and is thus sensitive to climate variations (Farbrot et al., 2007b).Ongoing global warming has been attributed as the trigger for several recent landslides from permafrost areas (Morino et al., 2019;Saemundsson et al., 2018).
Cold subglacial conditions, where the ice-bed interface is frozen and underlain by permafrost, can form in marginal areas where glacier cover is thin.The geomorphological imprint of cold or polythermal glacier margins can be different from temperate margins by a dominance of debris accumulation at the glacier front, preservation of buried glacier ice, and the formation of ice-cored moraines and composite ridges (Etzelmüller and Hagen, 2005;Waller et al., 2012).After permafrost thaw and ground ice melt, the landscape becomes dominated by hummocky terrain, often sharply terminated by a collapsed former ice-cored moraine ridge (Midgley et al., 2018;Sollid and Sørbel, 1988).
Iceland's sensitive location in relation to climatic and oceanographic transitions zones in the North Atlantic means that its glacial and periglacial environment, and subsequent imprint on the landscape, varied considerably throughout the Late Pleistocene and Holocene.Paleo-permafrost is a well-known phenomenon, originally described for unglaciated parts of central Europe during Pleistocene glaciations (Vandenberghe and Pissart, 1993), and known from early-deglaciated areas in Fennoscandia, such as in north and north-western Norway (Sollid and Sørbel, 1992;Svensson, 1988).In these areas, local cirque glaciation (Dahl and Nesje, 1992;Mangerud et al., 1979;Sollid and Sørbel, 1979), relict rock glaciers (Lilleøren and Etzelmüller, 2011) and ice-wedges on marine terraces (Svensson, 1988) bear witness to a former cold and mostly unglaciated environment.Finally, hummocky moraine areas, interpreted as down-wasted former ice-cored moraines, indicate ice margins terminating into a permafrost environment (Berthling et al., 2013;Hambrey et al., 1997;Sollid and Sørbel, 1988).Following the Pleistocene, the Holocene thermal maximum (HTM) and the Little Ice Age (LIA) are crucial periods for climate reconstruction.Multiple proxy data exist, showing a higher HTM temperature than displayed by the Greenland ice cores, with several independent studies indicating 2e3 C higher mean annual air temperatures (MAAT) in Iceland compared to the reference period 1961e90 (Caseldine et al., 2006;Flowers et al., 2008;Geirsd ottir et al., 2009).After 5 ka BP, several cool periods occurred, culminating in the LIA (Flowers et al., 2008;Geirsd ottir et al., 2013;St€ otter et al., 1999), with an estimated MAAT lowering of c. 0.8 C compared to recent decades (Anderson et al., 2018a).
To evaluate the dynamics of permafrost in Iceland since the onset of the last deglaciation, we apply the forcing and output of a 3D, coupled climate/ice-sheet model (Hubbard, 2006;Patton et al., 2017b) to a transient permafrost model (CryoGRID 2, Czekirda et al., 2019;Westermann et al., 2017;Westermann et al., 2013) from the Last Glacial Maximum (LGM, c. 20 ka BP) to present.The permafrost model was forced by either modelled subglacial temperatures under conditions of ice coverage, or surface air temperatures under ice-free conditions.The combined modelling approach provides insights into the age of permafrost in Iceland, identifies areas with widespread paleo-permafrost, and enables quantitative assessment of the persistence of permafrost and rates of degradation in the different areas.We discuss our results within the context of glacial and periglacial processes, landforms and landscape development.

Setting
Iceland is located between 63 N and 66 N at the northern part of the North Atlantic mid-ocean ridge.It is the emerged part of the Iceland Basalt Plateau (Thordarson and Larsen, 2007), and the onshore rocks are mainly Neogene and Quaternary basalts (Harðarson et al., 2008).Currently, volcanic activity takes place exclusively in the 15-50 km-wide neo-volcanic zones (Thordarson and H€ oskuldsson, 2008) (Fig. 1).
The Pliocene-Pleistocene volcanic succession of Iceland reveals evidence of numerous major glaciations and that the Pleistocene ice sheets covered most of the island and surrounding shelves (Eiriksson, 2008;Eiriksson, 1994).At the LGM (~20 cal ka BP), the Icelandic Ice Sheet (IIS) extended beyond the present coastline as far as the shelf break in all sectors (Fig. 1A, Andrews et al., 2000;Hubbard, 2006;Norðdahl et al., 2008;Patton et al., 2017b;Spagnolo and Clark, 2009).The widespread extension of the ice sheet into a marine setting, coupled with the highly elastic nature of the lithosphere in the region, meant the IIS was particularly susceptible to eustatic sea level rise during the Bølling-Allerød interstadial, forcing a rapid collapse from the shelf (Norðdahl and Ing olfsson, 2015).Its southerly latitude meant the IIS was also sensitive to atmospheric warming, reflected by the rapid expansion of the zone of surface melting across the low-aspect ice-sheet surface during this interstadial (Patton et al., 2017b).By the end of the Bølling-Allerød interstadial the ice sheet was reduced to a quarter of its LGM extent (Ing olfsson et al., 2010).Glacial advances occurred during the Younger Dryas stadial, reaching the innermost parts of the fjords, and probably included renewed cirque glacier activity in the mountains of Tr€ ollaskagi, while lowland areas were still ice-free (Ing olfsson et al., 1997(Ing olfsson et al., , 2010)).By 10.3 cal ka BP, the ice sheet was rapidly retreating across the highlands, and sub-aerially formed lava flows demonstrate ice-free conditions in the highlands by c. 8.6 cal ka BP (Geirsd ottir et al., 2009;Hjartarson and Ing olfsson, 1988;Kaldal and Víkingsson, 1990).During this period (HTM, c. 8-6 cal ka BP) MAATs 2e3 C were higher than the reference period 1961-90 and most remaining ice caps disappeared (Flowers et al., 2008).Glaciers re-advanced during the Neoglaciation (<6 ka BP) (Guðmundsson, 1997;St€ otter et al., 1999), reaching a maximum Holocene extent during the LIA (Geirsd ottir et al., 2009(Geirsd ottir et al., , 2013)).
The Icelandic climate today is maritime, with precipitation rates of more than 7000 mm a À1 in the south-east, decreasing towards north-west (Fig. 1D).The interior of Iceland is relatively dry, with precipitation rates <1000 mm a À1 .Surface MAATs between 1961 and 1990 are strongly dictated by elevation and latitude, with the southern and south-eastern coast being warmest (>8 C), and the high mountain areas and glacier surfaces being coldest (<-8 C, Fig. 1C).The climate has varied significantly since start of instrumental observations in the late nineteenth century.A warm period during the 1930s was followed by intermittently cooler temperatures until the early 1990s (Hanna et al., 2004) and since then MAATs have increased (Czekirda et al., 2019).
Today, the lower mountain permafrost limits in Iceland vary with climate, with lower limits in the maritime south-east at close to 1000 m a.s.l., decreasing towards 700 m a.s.l. in the more continental north (Czekirda et al., 2019;Etzelmüller et al., 2007).Sporadic permafrost is abundant in palsa mires and small peat plateaus in the central area of Iceland (Friedman et al., 1971;Priesnitz andSchunke, 1978, 1983;Saemundsson et al., 2012).Numerous intact rock glaciers were mapped by Guðmundsson (2000) and have later been discussed in terms of permafrost dynamics of the area (Farbrot et al., 2007a;Kellerer-Pirklbauer et al., 2008;Tanarro et al., 2019;Wangensteen et al., 2006).Currently, four boreholes monitor ground temperatures, of which three record permafrost conditions (Farbrot et al., 2007b).

Methods
The output of a 3D, time-integrated ice sheet model (Hubbard, 2006;Patton et al., 2017b) is used to force a transient permafrost model (CryoGRID 2, Czekirda et al., 2019;Westermann et al., 2013), from the LGM (c.20 ka BP) until today.In the following sections, the principles of the two models and their combination are described.

Principles
The 3D, time-integrated ice sheet model is based on a first-order approximation of the Stokes-equations adopted from Blatter (1995), Hubbard (2000Hubbard ( , 1997)), Marshall et al. (2005) and Pollard and Deconto (2007).It has previously been applied to Iceland (Hubbard, 2006;Patton et al., 2017b), the British Isles (Hubbard et al., 2009), Patagonia (Hubbard et al., 2005) and the Eurasian ice-sheet complex (Patton et al., 2016, Patton et al., 2017a) to investigate the build-up, extent and deglaciation of the paleo-ice sheets that occupied these regions.The approach to solving the membrane stress and associated strain fields equate to the L1L2 classification of higher-order models defined by Hindmarsh ( 2004) that includes longitudinal deviatoric stresses that act to stabilize ice flow over steep relief with high rates of basal lubrication and motion.The model is thermo-mechanically coupled (e.g.Hubbard et al., 2004), and temperature-dependent internal flow (ice deformation) is complimented by basal motion calculated using Weertman's (1964) sliding relation when subglacial temperatures attain pressure melting point.

Forcing
Surface mass balance is determined by a positive degree-day scheme, applied according to Laumann and Reeh (1993), and derives total melt from integrated monthly positive temperatures.Palaeo-climate forcing is implemented from the 50-year interval NGRIP d 18 O record (Andersen et al., 2004), scaled between a maximum prescribed temperature depression and present-day conditions.The inclusion of an additional bulk offset within the temperature scaling calculation controls the degree of fluctuation within the forcing record.Precipitation distributed evenly throughout the year, and accumulates as snow when the surface temperature falls below a threshold of 1 C.Both temperature and precipitation adjust to the evolving ice-sheet surface (corrected for isostatic loading) through applied lapse rates, derived from multiple-regression analyses of meteorological observations over the reference period (1961e1990 provided by the Icelandic Meteorological Office) (Bj€ ornsson et al., 2007;Crochet et al., 2007).To account for reduced precipitation during glacial conditions, a bulk precipitation reduction of 40e50% is imposed during the LGM and Late Glacial, with a further 35% reduction gradient across eastern sectors.These bulk perturbations are returned to the "present-day" type distribution after the YD/Holocene transition.Mass wastage at tidewater margins is calculated according to the calving-front geometry with an empirically derived water-depth and ice thickness algorithm (Brown et al., 1982), which though simplistic, does indirectly capture the effects of both fracture and submarine melting.The pattern of geothermal heat flux used as a basal thermal boundary condition was determined from temperature measurements taken at the base of shallow (30e60 m) boreholes (Flovenz and Saemundsson, 1993) that are interpolated across the model domain using a standard kriging gridding algorithm with a linear variogram model containing no anisotropic weighting.

Principles
The mathematical basis of the model is a one-dimensional transient heat conduction equation that additionally accounts for the latent heat effects due to ice-water phase changes:  effective volumetric heat capacity, thermal conductivity, depth, temperature and time, respectively.Equation ( 1) can be applied for spatially distributed permafrost modelling by calculating it for each grid cell, i.e. there is no interaction between the neighboring cells and the lateral heat flux is thus not considered in the model.Thermal properties of the subsurface layers are estimated as functions of the thermal properties of the individual soil constituents, such as water, ice, organic matter, minerogenic matter, air, and their respective volumetric fractions.The liquid water and ice contents are obtained in CryoGrid 2 from the functions by Dall'Amico et al. ( 2011), which link liquid water content, temperature and soil saturation degree.Variations in soil water content depend exclusively on freezing and thawing processes, i.e. water or water vapor movement in the soil, along with additional external water inputs (meltwater, rain) and their impact on the ground temperature are neglected in the model.

Sub-surface and model domain
Thermal conductivity of the rock matrix or mineral fraction was set to 2 W m À1 K À1 (Flovenz and Saemundsson, 1993).For each subsurface layer volumetric contents of water (q w ), mineral component (q m ), organic matter (q o ) and air (q a ) were defined as means for the whole modelling period.Subsurface stratigraphy was assigned according to the soil map of Iceland (Arnalds, 2008(Arnalds, , 2015)), with more details given in Czekirda et al. (2019).To estimate the overall depth of unconsolidated sediments in each soil class, we used zonal statistics over the global depth to bedrock data set (Shangguan et al., 2017), and based on the results we assumed a uniform regolith depth of 10 m to be an adequate approximation for all the soil classes.In general, the porosity of Icelandic basalt decreases horizontally with the distance from the spreading zone and vertically with depth (Flovenz and Saemundsson, 1993).We assumed a constant bedrock porosity of 8%, which represents average porosity of basaltic lavas in Iceland (Stef ansson et al., 1997).The model depth for each grid cell was 1 km.The discretisation with depth was logarithmic, starting with 5 cm intervals between 0 and 3 m depth, increasing to 1 m to 30 m depth, 10 m to100 m depth and 100 m to the bottom of the domain.The grid spatial resolution was 2 km.

Forcing and runs
The permafrost model is forced either by air temperature, when the surface is not covered by a glacier, or subglacial temperatures provided by the glacier model (Fig. 2a).The air temperature data set is the same used for the glacier forcing, on a spatial resolution of 2 km.We did account for a surface offset due to snow cover in these runs as we do not consider seasonality, and assume that the uncertainty of the forcing is larger than the possible bias due to lacking snow cover.The model results must therefore be considered a first-order-approximation of the true thermal regime.We produced several runs for the period 15 ka BP to present, where between 11 ka and 5 ka, the HTM was increased by 0.5 C, 1 C and 2 C, respectively (Fig. 2b).After 5 ka BP, the temperature forcing was the same for all runs.The model runs with 50-year increments during the whole modelling period, and seasonal temperature variations are not considered in these runs.

Lower boundary conditions
At the lower model boundary at 1000 m depth, we used the published heat flow data by Hjartarson (2015), assuming constant heat flux in time.The heat flux is generally high in Iceland, between 40 mW m À2 to over 300 mW m À2 (Fig. 2c), with the largest values occurring in, or near, the active volcanic zones.

Initialisation and thresholds
Given the long run times, we initialize the model domain with a homogenous ground temperature of 0 C; the first 2000 years of the LGM is regarded as a spin-up period.With high geothermal heat flux imposed across the domain (Hjartarson, 2015), the resulting permafrost sequences are not expected to become excessively thick.For this study, we calculated permafrost thickness over different time steps, along with permafrost duration, defined as the total time a model cell experienced permafrost during its period of sub-aerial exposure.A conservative threshold for permafrost of < -0.5 C at 10 m depth was chosen.

Landform inventory
In order to compare modelling results with landform distribution, landforms developed on slopes, and relevant for the topic were mapped.For simplicity, we distinguished between rock glaciers, intact ice-cored moraines and large slope debris bodies, commonly interpreted as postglacial landslides (J onsson, 1976;Whalley et al., 1983).We used Google Earth and orthophotoservices in ArcGIS and from the Loftmyndir Ltd. (www.map.is) to identify and classify the features, with the help of published maps (Mercier et al., 2017) and existing inventories (Guðmundsson, 2000;Lilleøren et al., 2013).The classification was subjective, and the quality of the imagery varies over Iceland.Thus the inventory is only used qualitatively in this study.

Modelled time series and run validation
Paleo-permafrost distribution in Iceland is directly associated to the extent of the Icelandic Ice Sheet (IIS) at a given time and climate forcing.During glaciation, subglacial temperatures attained a maximum equivalent to the pressure melting point (PMP), determined by the thickness of the ice column above.As temperatures decrease below this PMP a transition to cold-based ice conditions occurs whereby the glacial ice becomes frozen to the substrate.This latter scenario is especially valid for high-mountain cells and marginal areas, where glacier thickness is lower and cold ice can be easily advected towards the bed.At these sites, we model a shallow permafrost body under the glacier (Fig. 3a), normally less than 10 m thick.
The glacier extent decreased from total coverage of the land areas (with the probable exception of some coastal nunataks) during the LGM, to almost ice-free conditions during the HTM.Permafrost coverage varied between 0 km 2 and c. 20,000 km 2 (more than 20% of the land area) according to our model (Fig. 3b).Climate and glacier coverage determine permafrost aggradation and degradation, which may have occurred several times during and after glacier advances.During the late Pleistocene, glacier reduction triggered permafrost aggradation, while after the HTM, glaciers and permafrost increased in area in concert (Fig. 3b).Hence, the dominant factor restricting permafrost development switches from glacial extent and subglacial thermal regime before the Holocene to climate forcing during the Holocene.
The penultimate time slice of the model run (100 yr before the reference period 1961e1990) was compared against results from Czekirda et al. (2019) which used the CryoGrid2 model with a climate forcing independent of the current study (Fig. 3c).The mean difference between the two runs was þ0.7 C, with a standard deviation of ±0.8 C. The spatial distribution of ground temperatures show satisfactory agreement, considering the 2 km resolution and the Greenland ice-core paleo-isotopic proxy data used for forcing.

Paleo-permafrost during the Late Pleistocene
During deglaciation, land areas became sub-aerially exposed.If above sea level, permafrost would aggrade as the MAAT was still much colder than present or during the LIA (for Iceland, c. 8e10 C cooler).The outermost land areas in the Northwest (Westfjords) and the Northeast and Eastern part of Iceland were most affected by this permafrost aggradation (Fig. 4a).During the warm Bølling-Allerød interstadial (c.15.4e13.0cal ka BP), the ice sheet decreased rapidly, and early coastal permafrost degraded close to the coast again (Fig. 4b).Maximum permafrost distribution probably occurred during the onset of the Younger Dryas (YD), following the rapid deglaciation during the Bølling-Allerød interstadial, and ongoing trends of atmospheric cooling towards the YD (Fig. 4c).This period probably experienced the fastest permafrost aggradation (Fig. 3c), with permanently frozen ground expanding especially in the ice-free parts of the mountainous coastal areas around Iceland (Fig. 4c).During this period, permafrost aggradation exceeded 80 m in the coldest areas, which is, however, relatively shallow because of the influence of high ground heat flow (Fig. 4d).Estimations of the duration of a cell being frozen indicate that many areas, especially along the northern, western and eastern coast, accumulated periods of permafrost conditions of more than 4 ka between the LGM and the onset of the Holocene (Fig. 4e).

Paleo-permafrost during the Holocene
When applying an HTM scenario which is >2 C warmer than the last reference period 1961e1990, modelled permafrost disappears almost completely during this period (Fig. 5a and b), while later expansion to a maximum extent occurred during the LIA (Fig. 5c).After the HTM, permafrost probably aggregated at a steady pace in concert with the Neoglacial cooling towards the LIA (Fig. 3a).The modelled distribution of permafrost follows the pattern known from present conditions, with the main permafrost areas in the mountains in central and northern Iceland.In these areas, permafrost has probably been present for more than 4 ka (Fig. 5d).The modelling suggests that the highest areas of the Tr€ ollaskagi peninsula could host relict permafrost of Pleistocene age, which is similar to the Norwegian mountains (Lilleøren et al., 2012).

Uncertainties
The 2-km resolution of the finite-difference grid directly influences the representation of topography, particularly in areas of high relief.It is therefore likely that areas of cool and possibly thick permafrost are underestimated in mountain areas.The insulating effect of snow cover is an important factor for permafrost temperatures and stability in winter, giving rise to considerable surface offsets between air and ground surface temperature.The Greenlandic paleo-climatic forcing also has large uncertainties and is limited by the assumed link to the NGRIP ice-core isotopic record and its temporal resolution.The 1 m depth temperature for the last time step compares well with a model experiment forced by gridded climate data (Czekirda et al., 2019) with similar values, so we are confident that our results demonstrate a realistic pattern within the framework of the forcing data quality.Furthermore, we have not modelled any seasonality.This would normally cool the ground further, depending on water/ice content and substrate characteristics due to differences between frozen and thawed thermal conductivity of the ground substrates and water/ice (thermal offset) (Williams and Smith, 1989).A seasonal forcing also considerably delays permafrost thaw due to the thermal heat capacity of ice (Lilleøren et al., 2012).Further, vegetation significantly alters ground temperatures.Today, vegetation is very low or absent in Iceland's high mountains, and we assume this was also the case during the Pleistocene.Finally, areas such as block fields and coarse slope deposits, including rock glaciers, have a cooling effect on ground temperatures because of air advection (e. g.Ballantyne, 2018).This effect is important for local landform development, but neglected here as such areas are usually small and insufficiently captured in the ground resolution used.
A major control on the pattern and timing of permafrost development is the simulated deglaciation of the modelled IIS.A paucity of geomorphological, chronological and climate proxy data constraining ice evolution after the YD and through the HTM increases uncertainties regarding the accuracy of the modelled ice extents and volumes.Although there is a generally good agreement between modelled and empirically inferred ice extents during the Younger Dryas re-advance, there are regional disparities (Patton et al., 2017b).These include ice-free conditions over the Vestfirðir and Snaefellsnes peninsulas and a restricted ice-extent across NE Iceland, highlighting the climatically sensitive nature of these peripheral ice masses (Harning et al., 2016) in combination with the relatively simplistic climate parameterisation and positive degreeday scheme of the ice sheet model that determines surface mass balance.The final disintegration of the IIS during the Early Holocene proceeded rapidly, potentially to a state where ice caps were smaller than the present day coverage during the HTM or even disappeared (Anderson et al., 2018a;Geirsd ottir et al., 2009).Given the sparsity of geophysical constraints during this time, we have presented a suite of post-Younger Dryas deglaciation/climate scenarios to capture the potential sensitivity of the permafrost distribution to glacial evolution during this time.
The distinct geothermal boundary condition of the Icelandic domain is a primary control determining the basal environment of the ice sheet, effective viscosity and drainage properties (Patton et al., 2017b).Under elevated geothermal fluxes, the thermal regime of the IIS tends towards a state of widespread basal melting, promoting stable and continuous fast-flow drainage, increased frictional heating, and a relatively shallow ice-surface gradient.While perturbations of the geothermal flux have a minor influence on the ice sheet's architecture, it does affect the persistency of coldbased conditions and thus permafrost aggradation at highelevation sites where ice tends to be thinnest (Patton et al., 2017b).
Despite these data-model uncertainties, we regard the predicted spatial pattern and the interaction between glacier retreat and permafrost aggradation and degradation a useful first-order approximation for discussing geomorphological implications of the permafrost history.In the following section, we discuss some major patterns identified through this analysis, which might aid the understanding of periglacial landform development in Iceland.

Modelled permafrost age
Besides using modelling tools, estimating the age of permafrost is only possible by dating certain landforms clearly associated to permafrost, such as ice wedges (Opel et al., 2018;Vasil'chuk and Vasil'chuk, 1997), via paleo-ecological studies (Kjellman et al., 2018) or direct dating of ice samples by using organic material in the ice or carbon in ice bubbles (Ødegård et al., 2017).Ages of permafrost landforms are available from the Arctic permafrost area of Alaska (Hamilton et al., 1988), Canada (Murton, 2005) and Siberia (Vasil'chuk and Vasil'chuk, 1997) while in mountain areas, quantitative permafrost age estimations based on measured information are rare (Haeberli et al., 2003;Krainer et al., 2015).Most studies are related to palsas and peat plateaus, where certain plant assemblages indicate the onset of permafrost in mires (Kjellman et al., 2018;Kuhry, 2008;Oksanen, 2006;Sannel and Kuhry, 2008).These sporadic permafrost sites are mostly dated to the LIA in Scandinavia; however earlier Holocene periods of permafrost development have also been described (Kjellman et al., 2018).Concerning mountain landforms, rock glaciers have been dated, either directly (Haeberli et al., 1999;Ivy-Ochs et al., 2007), indirectly via surface movement analyses (Wangensteen et al., 2006), or relatively by dating tools such as lichenometry (Andr e, 1994) and Schmidt hammer analysis (Kellerer-Pirklbauer et al., 2008;Winkler and Lambiel, 2018).Many of these studies indicate rock glacier activity in the Neoglacial, though in many cases, formation seems to have started during or following Late-Pleistocene deglaciation.(Krainer et al., 2015;Moran et al., 2016).
Our study indicates that between the onset of the last deglaciation and the Pleistocene/Holocene transition, early ice-free areas above the marine limit, especially along the western, northern and eastern mountainous coast, must have experienced subaerial permafrost aggradation over several thousands of years (Fig. 4d), even in the lowlands.This near-coast permafrost degraded after the YD in response to atmospheric temperature rise.These Pleistocene permafrost conditions in mountain areas are well-known from the literature, e.g., from central Europe and areas south of the large ice sheets in North America (Dyke, 2004) and Eurasia (Hughes et al., 2016;Patton et al., 2017a;Stroeven et al., 2016).In northern-most Norway, ice-wedges (Svensson, 1988) and hummocky moraines interpreted as remnants of ice-cored moraines (Sollid and Sørbel, 1988) are found distal to the YD main moraine, indicating cold conditions and permafrost aggradation during deglaciation.Comparisons of recent glacial marginal landform assemblages with YD end moraines in Scotland indicate permafrost conditions outside the ice sheets (Hambrey et al., 1997).Relict rock glaciers in early deglaciated parts of the north Norwegian coast indicate permafrost conditions there, too (Lilleøren and Etzelmüller, 2011).These environments are comparable to Iceland, where we would expect to find remnants of Pleistocene permafrost landforms and glacierpermafrost interaction.
During the Holocene, the Pleistocene permafrost degraded in the outermost coastal areas, while at the same time it developed in the central and eastern Icelandic highlands, following more or less the recent distribution (Fig. 5).Similar to the remnant ice caps, it is likely that very little permafrost survived the HTM (Fig. 5b), while the largest areal distribution thereafter occurred during the LIA (Fig. 5c).With this pattern in mind, we delineated three major probability classes of permafrost in Iceland following the classification proposed by Lilleøren et al. (2012) (Fig. 5d and e): (1) Permafrost which never thawed during the Holocene, and thus could be very old; (2) Permafrost formed after the HTM during the Neoglaciation, and (3) Permafrost formed during the LIA.Permafrost has more recently warmed in all classes, while the youngest LIA permafrost is probably degrading today.
Reliable ages of permafrost ice are sparse to non-existent in most mountain areas.Ice at the bottom of a cold ice-patch has been dated to 7.6 cal ka BP in Jotunheimen in Norway (Ødegård et al., 2017), opening the possibility for permafrost survival through the HTM in comparable settings.However, to our knowledge, permafrost ice in the ground has never been absolutely dated in Scandinavia.The latter class (3) also covers areas with modern palsa mires, which have been attributed to LIA ages (Friedman et al., 1971), and are probably degrading today (Saemundsson et al., 2012;Thorhallsdottir, 1996).This degrading permafrost potentially triggered the frequent landslides observed during recent years (Saemundsson et al., 2018).
It is noteworthy that we model different dominating factors for permafrost aggradation and degradation before and after the onset of the Holocene.During the Late Pleistocene, the dominant factor restricting permafrost development is the glacier extent and basal thermal regime, switching to solely climate forcing during the Holocene (Fig. 3b).This is probably a pattern valid also for other areas influenced by glaciations.

Rock glaciers and rock slide deposits
Rock glaciers, palsas and peat plateaus indicate permafrost in mountainous areas.Palsas and peat plateaus were recognised early-on as permafrost indicators in Iceland (Friedman et al., 1971;Priesnitz and Schunke, 1978;Thorhallsdottir, 1994).However, rock glaciers have been subject to discussions and confusion (Berthling, 2011), both as a landform and as a permafrost indicator, especially in Iceland (Whalley et al., 1995a;Whalley and Martin, 1994).Today, palsas are abundant in the Central Highlands, and give rise to the permafrost classification in the International Permafrost Association (IPA) map (Brown et al., 1995).In this study, the palsa distribution is not further assessed in relation to our results, and for further discussion in the permafrost modelling context, we refer to Czekirda et al. (2019).
Iceland has a very high number of rock glaciers, especially in the mountain areas of the Tr€ ollaskagi and Flateyjarskagi peninsulas (Guðmundsson, 2000;Lilleøren et al., 2013;Tanarro et al., 2019) (Fig. 6).The high abundance of these landforms in Iceland is related to the availability of loose weathering material, and a topography characterized by steep, glacially eroded slopes (Lilleøren et al., 2013).There is consensus today that active rock glaciers contain ice of varying origin, and that these surface structures are expressions of cohesive flow in perennially frozen ice-rich materials.They can therefore be used as indicators of local to regional permafrost occurrence (Barsch, 1996;Berthling, 2011;Haeberli et al., 2006).Recently, studies have focused on rock glaciers as aquifers and their  2012).The results may indicate survival of permafrost in some areas during the HTM (1), permafrost developed during the Neo-glaciation after the HTM (2) and during the LIA (3).Red dots mark recently observed and described landslides in permafrost slopes (Czekirda et al., 2019;Morino et al., 2019;Saemundsson et al., 2018).(F): Forcing temperatures (Run279 and 282) and position of time plots from (A)e(C).All glacier outlines represent glacier coverage today.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)relevance to mountain hydrology (Jones et al., 2018).Relict rock glaciers are indicative of past permafrost conditions, and have been used for paleo-permafrost mapping in mountains through, e. g. statistical prediction methods (Brenning et al., 2007;Brenning and Trombotto, 2006;Marcer et al., 2019;Sattler et al., 2016).
perspective on the distribution of rock glaciers across Iceland.

Active rock glaciers
Active rock glaciers are considered a product of millennia-long creep processes, with their age determined by their size and flow velocity.For relict rock glaciers, a Late Glacial age is often inferred (Barsch, 1996;Frauenfelder and K€ a€ ab, 2000).In addition, rock glaciers often consist of what can be interpreted as several generations of landforms on top of each other, as described, e.g., from H olar in Hjaltadalur in western Iceland (Farbrot et al., 2007a;Kellerer-Pirklbauer et al., 2008) (Fig. 6a, image (2)).The active rock glacier distribution is in agreement with the area modelled as Holocene permafrost (Figs.6a and 7e), both in this and in previous studies (Czekirda et al., 2019;Etzelmüller et al., 2007).The calculated age of permafrost in this study agrees with other studies, indicating an age of around 3e5 cal ka BP (Figs. 5d and 7e), coinciding with the onset of the Neoglaciation in Iceland.

Relict rock glaciers
Relict rock glaciers at lower elevations are older than the currently active ones, and formed during a long period without extensive glacier coverage within a permafrost environment.Therefore, relict rock glaciers have also been used to constrain ice sheets and deglaciation patterns (Bakke et al., 2005).In our inventory, we find relict rock glaciers situated close to sea level, especially at the northernmost coastline and partly in the Westfjords (Fig. 6a).They often terminate in flat areas close to sea level, and this kind of rock glacier-strandflat system is also frequently observed on Svalbard (Berthling et al., 1998(Berthling et al., , 2000;;Sollid and Sørbel, 1992).Near-shore rock glaciers in Iceland are probably older than the HTM.Our interpretation is that there has been a time window of ice-free conditions long and cold enough to facilitate rock glacier development during the deglaciation of the last IIS.Our study clearly demonstrates that permafrost conditions lasted long enough in coastal areas in northern, western and eastern Iceland for such profound landforms to develop between the onset of deglaciation and the HTM (Fig. 4e).Rapid deglaciation of the continental shelf in response to rising eustatic sea level also allowed periglacial conditions to dominate coastal mountain slopes for more than 5 ka (Ing olfsson et al., 2010).Deglaciation was interrupted by the Younger Dryas re-advance, which is poorly constrained in Iceland (Norðdahl et al., 2008), and probably did not affect the Tr€ ollaskagi area beyond local glacier advances in the innermost fjords (Ing olfsson et al., 1997;Van Vliet-Lano€ e et al., 2007).These time constraints should provide enough time to produce rock glaciers (Barsch, 1996).

Slope debris bodies and rock slides
On the steep slopes of Tr€ ollask agi, in the Borgarfj€ orður area in eastern Iceland, and in the Westfjords, large debris bodies cover the slopes and sometimes reach the valley floors (Fig. 6b), with an abundance of these landforms increasing towards the coast.These landforms are all mapped as landslides by J onsson (1976) and later by Whalley (1983) and Mercier et al. (2017), triggered as a result of glacial debutressing and possible seismic activity related to isostatic rebound (Cossart et al., 2014).The Icelandic crust is highly sensitive to glacial loading and unloading, and small and gradual changes in ice mass are compensated for by rapid isostatic movements (Biessy et al., 2008;Ing olfsson and Norðdahl, 2001;Ing olfsson et al., 1995).Recent dating of some of these landforms based on tephrochronology indicates a minimum age of between 12 and 7e8 cal ka BP (Decaulne et al., 2016;Mercier et al., 2017).A possible role of permafrost aggradation and degradation in the development of these spectacular landforms has not been considered in published scientific literature.
However, another hypothesis defines these landforms as rock glaciers, and thus periglacial features which developed over a long time period (Guðmundsson, 2000).Many of these landforms have clear lobe forms, show creep patterns with furrows and ridges at the surface, and show movement even today (J onsson and Agústsson, 2007;Wangensteen et al., 2006).Comparing the position of these debris bodies with our results indicates that (1) many of these landforms lie in areas that deglaciated early and with permafrost beneath the surface (Fig. 6), and (2) the age of these landforms seem to fall into a period where the slopes were warming from a permafrost to a non-permafrost situation after the YD towards the HTM (Fig. 7f).
We suggest that the mechanism of isostatic rebound combined with seismic activity could make large unconsolidated masses available for secondary creep processes in Iceland, and provide one possible explanation for the many large slope deposits close to sea level such as rock glaciers and other creeping slope debris bodies.Based on our observations, we hypothesize that as the deglaciation of mountain areas close to the west, north and east coasts of Iceland rapidly proceeded during the Bølling-Allerød interstadial, there was enough time to aggregate ground ice in slope debris bodies by subsurface freezing processes under the influence of snow meltwater or the burial of avalanche snow (Humlum et al., 2007;Knight et al., 2019), facilitating creep also after deposition of debris.

Holocene permafrost, recent slope dynamics and hazards
Following the HTM, the Holocene is dominated by a cooling trend that culminated in the LIA.Since then, rapid atmospheric warming has occurred, especially during the last two decades (AMAP, 2017;Corell, 2005).Under present trajectories, it is apparent that the youngest permafrost formed during the LIA will continue to degrade.The rate of permafrost degradation is highly dependent on ground substrate and water/ice content (including unfrozen water content).However, high geothermal heat flux may decrease this difference (Nicolsky and Romanovsky, 2018).Thus, in Iceland and other permafrost areas with relatively high geothermal heat fluxes such as the Kamchatka peninsula in eastern Russia (Abramov et al., 2008) or the peaks of Hawaii (Schorghofer et al., 2017), permafrost is thin, permafrost temperatures are close to the melting point, and degradation occurs rapidly in response to climate perturbations (Etzelmüller et al., 2007;Farbrot et al., 2007b).During the last years, several landslides from frozen mountain slopes have been observed in Iceland (Saemundsson et al., 2018), mobilising large bodies of frozen sediment (Fig. 5e).Subsequent thawing of this debris formed a hummocky terrain, called molards, which recently has been discussed as an indicator of permafrost degradation (Morino et al., 2019).We speculate that the present thawing of LIA permafrost is one of the reasons for these events.Numerical modelling at these sites indicates a relation between recent permafrost degradation and weakening of the slope due to permafrost warming (Czekirda et al., 2019).As these landslides have happened in unconsolidated sediments, the temperature increase and subsequent weakening of frozen rock strength (Davis et al., 2001) in the ground probably exacerbated conditions for failure.

Conclusions
From this study, we draw the following conclusions: Permafrost was widespread in the early deglaciated areas of western, northern and eastern Iceland after the LGM.Up to 20% of the Land area of Iceland was underlain by permafrost during the Late Pleistocene, with major occurrences (C): Modelled ground temperatures at 10 m depth and ice coverage at 14.5 cal ka BP in the Tr€ ollaskagi area, (D) modelled 14.2 cal ka BP situation.The outermost landforms are probably older than the innermost, and also here is an initialisation of landform generation proposed to the same time period as in the Westfjords.(E): Permafrost duration since the onset of the Holocene (Fig. 6), related to active permafrost landforms.Most of the landforms cluster in the area of permafrost probably started to form after the HTM, and having an age of <5 ka.(F): Time-depth temperature plot for sites with dated landslides (red arrows, image 8 and 9 in Fig. 6, (Mercier et al., 2017).Upper: site image 8 (Fig. 6) at top of mountain wall, Middle: as upper, but on tow of deposit.Lower: site image 9 (Fig. 6), toe of deposit.The boxes indicate likely ranges of deposition (Mercier et al., 2017).The arrows indicate the modelled end of glacier coverage of the particular cell.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) in the coastal mountain areas of western, northern and eastern Iceland.Due to its unique maritime and geothermal setting, during the last glacial-interglacial transition, Iceland's permafrost was relatively shallow and highly sensitive to climatic perturbations.Before the Younger Dryas (YD), permafrost aggradation happened during phases of deglaciation, after which the YD glacier extent and permafrost varied in concert.This pattern demonstrates that the dominant environmental control restricting permafrost development switches from glacier extent and subglacial thermal regime to sub-aerial climate forcing.Relict rock glaciers are abundant in the coastal areas of Iceland.Their formation is related to long-term debris accumulation associated with the early collapse of the marine-based ice sheet that dominated these areas, and subsequent rapid aggradation of permafrost.Model results demonstrate that large areas were underlain by permafrost over multiple millennia, facilitating long-term landform development and influencing the stability of steep slopes.During the HTM, permafrost and glaciers were mostly absent, except in the highest areas of Tr€ ollaskagi.Permafrost depths attained their Holocene maximum during the LIA and have since been degrading, contributing to slope instabilities and the recent increase in slope failures observed across Iceland.

Fig. 1 .
Fig. 1. (A) Overview map of Iceland identifying the place names used in this study.The black circles are existing boreholes equipped with ground temperature loggers since 2004 (Farbrot et al., 2007b).An empirical reconstruction of the LGM ice margin is shown (after Norðdahl and Ing olfsson, 2015).(B) Map of the volcanic zones of Iceland with indication of fissure swarms and central volcanoes.Modified after Guðmundsd ottir et al. (2018).(C) Mean annual air temperatures (MAAT) for the normal period 1961e1990 (modified after the Icelandic Meteorological Office).(D) Mean annual precipitation (MAP) for the normal period 1961e1990 (modified after the Icelandic Meteorological Office).

Fig. 2 .
Fig. 2. Boundary conditions for permafrost modelling: (A) Upper boundary condition -Surface-forcing temperatures, either sub-aerial if not glacier-covered or subglacial if glaciercovered.Time step ¼ 14.3 cal ka BP.The blue dotted line indicates the modelled glacier boundary.(B) Upper boundary condition e forcing paleo time series transformed and adapted for Iceland on the 50-year interval NGRIP d 18 O record (Andersen et al., 2004).The forcing SAT (y-axis) denotes the linear deviation from the 1961e90 normal period.The runs 279e282 represent different HTM temperature forcings.Blue shaded areas indicate cool periods (LGM, YD, LIA), red areas a indicate warmer period (Bø-Al ¼ Bølling-Allerød, HTM).(C) Lower boundary conditione the ground thermal heat flux (re-drawn based on Hjartarson (2015).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 .
Fig. 3. (A) Example of a 1D plot showing forcing temperature (Run279, upper panel) and resulting ground temperatures with depth (lower panel).The modelled cell is located at Gagnheiði in eastern Iceland (c.800 m a.s.l., position indicated in 3c) in eastern Iceland, where we measure ground temperatures in a borehole (Farbrot et al., 2007b).The land cover is indicated with coloured bars in the upper panel.(B): Modelled glacier and permafrost coverage of Iceland between 15 ka BP and today.The red shaded areas indicate the warm periods of Bølling-Allerød and the HTM, respectively.Note that permafrost aggradation in the Pleistocene is governed by glacier retreat; while after the HTM permafrost aggradation is governed by atmospheric cooling.(C): Upper map: Modelled ground temperatures in 2 m depth after Run279, last time step (À50 yr before 1990), Lower map: Modelled ground temperature based on a transient run from Czekirda et al. (2019), average for 1970e1979.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4 .
Fig. 4. Late Pleistocene glacier and permafrost conditions in Iceland.(A)e(C): Glacier outline and ground temperatures (GT) at 10 m depth for 16 cal ka BP, 14 cal ka BP and 11.7 cal ka BP. (D): Modelled permafrost thickness (PT) at the end of the YD (11.7 ka).(E): Calculated permafrost duration (PD) summarised over the period 20-10 cal ka BP.The glacier outline is modelled for 10 ka BP.We identified cells with T < -0.5 C at 10 m depth, and sum up the time steps with this condition over the modelled period.(F): Forcing temperatures (Run279) and position of time plots from (A)e(C).B-A ¼ Bølling-Allerød interstadial, H ¼ Holocene.

Fig. 5 .
Fig. 5. Holocene permafrost conditions in Iceland.(A), (B): Ground temperatures (GT) at 10 m depth for 7.5 cal ka BP (HTM) for Run279 and Run282, respectively.(C): Ground temperatures at 10 m depth for the LIA (0.25 ka BP) for Run279.(D): Calculated permafrost duration (PD) summarised over the period 10 ka BP to today for the most extreme run Run282.(E): Based on the Run282, three classes of Holocene permafrost distribution are identified, following the classification by Lilleøren et al. (2012).The results may indicate survival of permafrost in some areas during the HTM (1), permafrost developed during the Neo-glaciation after the HTM (2) and during the LIA (3).Red dots mark recently observed and described landslides in permafrost slopes(Czekirda et al., 2019;Morino et al., 2019;Saemundsson et al., 2018).(F): Forcing temperatures (Run279 and 282) and position of time plots from (A)e(C).All glacier outlines represent glacier coverage today.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .
Fig. 7. Ground temperatures (GT), permafrost duration and (PD) permafrost-related landforms in the Westfjords and Trõllaskagi.(A) Modelled ground temperatures at 10 m depth and ice coverage at 14.5 cal ka BP in the Westfjords, in relation to possible relict permafrost landforms (circles, legend see Fig. 6) and (B): modelled permafrost duration of the outermost areas of the Westfjords between LGM and HTM.Most of these landforms are situated outside of the 14,5 ka ice cap and may have initiated during early Bølling-Allerød.(C):Modelled ground temperatures at 10 m depth and ice coverage at 14.5 cal ka BP in the Tr€ ollaskagi area, (D) modelled 14.2 cal ka BP situation.The outermost landforms are probably older than the innermost, and also here is an initialisation of landform generation proposed to the same time period as in the Westfjords.(E): Permafrost duration since the onset of the Holocene (Fig.6), related to active permafrost landforms.Most of the landforms cluster in the area of permafrost probably started to form after the HTM, and having an age of <5 ka.(F): Time-depth temperature plot for sites with dated landslides (red arrows, image 8 and 9 in Fig.6,(Mercier et al., 2017).Upper: site image 8 (Fig.6) at top of mountain wall, Middle: as upper, but on tow of deposit.Lower: site image 9 (Fig.6), toe of deposit.The boxes indicate likely ranges of deposition(Mercier et al., 2017).The arrows indicate the modelled end of glacier coverage of the particular cell.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)