# Modelling the thermal dynamics of permafrost talus slopes: insights from a recently destabilised site (Gislá landslide, October 6 th 2020, Iceland)

Meven Philippe, Florence Magnin, Jean-Yves Josnin, Costanza Morino, Nicolas Monzie, Skafti Brynjólfsson

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Meven Philippe, Florence Magnin, Jean-Yves Josnin, Costanza Morino, Nicolas Monzie, et al. Modelling the thermal dynamics of permafrost talus slopes: insights from a recently destabilised site (Gislá landslide, October 6 th 2020, Iceland). Géomorphologie: relief, processus, environnement, 2025, 30 (3), 10.4000/134ad. hal-05393906

HAL Id: hal-05393906

https://hal.science/hal-05393906v1

Submitted on 2 Dec 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Modelling the thermal dynamics of permafrost talus slopes: insights from a recently destabilised site (Gislá landslide, October 6th 2020, Iceland)

Meven Philippe<fnref n="1" />, Florence Magnin<fnref n="1" />, Jean-Yves Josnin<fnref n="1" />, Costanza Morino<fnref n="2" />, Nicolas Monzie<fnref n="1" />, Skafti Brynjólfsson<fnref n="3" />

## Abstract

The degradation of permafrost enhances the occurrence of landslide in mountain environments. Therefore, predicting permafrost boundaries is crucial to mitigate hazards to populations and infrastructures. In the present study we investigate the Eyjafirði talus slope (Tröllaskagi peninsula, Iceland), from which a landslide occurred in 2020. Although this talus slope is located outside of the permafrost climatic boundaries, geomorphological evidence points towards an involvement of degradation of azonal permafrost in the generation of the landslide. However, the thermal dynamics of talus slopes is currently poorly understood, as air convection (the ‘chimney effect’) can play a role in the persistence of permafrost at the base of talus slopes.

We use the commercial software FEFLOW to run two-dimensional (2D) physical-based simulations of heat transfer within a cross-section of the Eyjafirði talus slope. We create forcing temperature time series by statistical reconstruction (back to November 1881) and modelling (back to – 20,000 years), based on temperature measurements acquired in 2021-2022. We explore the sensitivity of our model to i. the initial porosity/ice content of the talus slope (0.3, 0.5 and 0.8), and ii. the thermal conductivity (TC) of the rock phase (0.75, 1.1 and 1.75 W.m⁻<fnref n="1" />.K⁻<fnref n="1" />).

Our temperature measurements show that a chimney effect indeed occurs at the Eyjafirði talus slope. Although our modelling approach does not simulate air convection itself, permafrost indeed persists at the base of the talus slope in all model scenarios. Increasing the initial porosity/ice content and decreasing the TC of the rock phase generates a larger zone with sub-zero temperatures at the base of the talus slope, which leads to a larger zone where ground ice persists with an overall higher ice content. Our approach is unconventional as we initially know that ground ice was present in the Eyjafirði talus slope at the time of the landslide; it allows us to attest that the permafrost dynamics in the talus slope is best represented by our most ice-conservative scenario – i.e. with a TC = 0.75 W.m⁻<fnref n="1" />.K⁻<fnref n="1" />. Such a low TC is atypical based on the literature; we suggest that this low value could have partly compensated for the cooling effect of air convection (i.e. chimney effect).

<footnote n="1">¹EDYTEM, Université Savoie Mont-Blanc, CNRS UMR 5204, Le Bourget-du-Lac, France.</footnote>

<footnote n="2">²Università degli studi di Padova, Italy.</footnote>

<footnote n="3">³Icelandic Institute of Natural History, Akureyri, Iceland.</footnote>
# 1. Introduction

It has been known for decades that degradation of mountain permafrost can lead to ground destabilisation, and ultimately generate landslides (e.g. Allen et al., 2009; Cathala et al., 2024; Geertsema et al., 2006; Dramis et al., 1995; Fischer et al., 2006, 2011; Gruber and Haeberli, 2007; Karjalainen et al., 2019; Kellerer-Pirklbauer et al., 2012; Lyle et al., 2004; Magnin et al., 2023; Patton et al., 2019; Penna et al., 2023; Ravanel et al., 2017; Sæmundsson et al., 2018). Hence, in a context of global climate change, landslides represent an increasing hazard for populations and infrastructures in mountain environments (e.g. Borgatti and Soldati, 2013; Haeberli et al., 2010, 2017; Harris, 2005; Harris and Vonder Muhll, 2001; Hjort et al., 2022; Lyle et al., 2004; Morino et al., 2019; Petley, 2012; Smith, 1990; Smith and Riseborough, 1996). It is thus crucial to predict the boundaries of permafrost, in order to be able to locate areas that present a landslide susceptibility.

Predicting permafrost boundaries is allowed by numerical modelling of ground temperatures. Notably, statistical permafrost models are calibrated with topoclimatic data, and hence allow to model permafrost distribution at the regional or even global scale (e.g. Boeckli et al., 2012a, 2012b; Czekirda et al., 2019; Deluigi et al., 2017; Gisnås et al., 2013, 2014, 2016; Gruber, 2012; Gruber and Hoelzle, 2001; Guglielmin et al., 2003; Magnin et al., 2015, 2019; Nelson, 1986; Obu et al., 2019; Wang et al., 2019). It should be noted that previous studies developed methods to downscale climatic data, which allows to perform statistical modelling studies at a more local scale (e.g. Fiddes and Gruber, 2012, 2014; Fiddes et al., 2022; Filhol et al., 2023; Levavasseur et al., 2011; Stendel et al., 2007; Yin et al., 2021; Zhang et al., 2020). However, statistical permafrost models are unable to locate azonal permafrost (i.e. permafrost located beyond its climatic boundaries), as its presence depends on local physical processes (e.g. air circulation) that are not taken into account in statistical models.

Physic-based permafrost models differ from statistical models in that they compute physical processes (e.g. heat transfer, hydrology, phase changes of water, air convection) in order to evaluate the presence of permafrost in a given context (e.g. Cathala et al., 2024; Dagenais et al., 2020; Jafarov et al., 2012; Magnin and Josnin, 2021; Magnin et al., 2017; Staub et al., 2015; Stocker-Mittaz et al., 2002; Westermann et al., 2013, 2016; Wicky and Hauck, 2017; Wicky et al., 2024; Zhang et al., 2013). Therefore, they are theoretically able to model the presence of azonal permafrost. Due to the numerical challenge that represents the modelling of such complex physical processes, physic-based models are often restricted to the local scale and simple case studies.

In this scope, the Permolards project studies molards in order to better understanding the link between permafrost degradation and landslide generation. Molards are cones of loose debris found in landslide deposits in periglacial areas, which were recently identified as remnants of ice-cemented blocks of sediment that were transported within a landslide and then degraded in place (formally described in Morino et al., 2019; see also Brideau et al., 2009; Lyle et al., 2004; Milana, 2016). Molards are particularly useful in permafrost landslide studies, as they represent unequivocal evidence of the implication of permafrost degradation in the generation of a given landslide (Morino et al., 2019). It was notably observed that unpredicted molard-bearing landslides in Iceland originated from talus slopes (Figure 1). It suggests that permafrost degradation is involved in the occurrence of such landslides (Morino et al., 2019, 2021), and hence that azonal permafrost could be present within the concerned talus slopes.

However, the dynamics of permafrost within talus slopes is currently poorly understood. Indeed, a talus slope is a complex system where a peculiar air circulation can take place, called the ‘chimney effect’. Such air circulation was suggested to generate, at the base of the talus, a zone of maximum cooling in winter and of minimum warming in summer (Figure 1; e.g. Delaloye and Lambiel, 2005; Germain and Milot, 2024; Lambiel and Pieracci, 2008; Millar et al., 2014; Morard et al., 2010; Phillips et al., 2009; Popescu et al., 2017; Wicky and Hauck, 2017; Wicky et al., 2024). Hence, the chimney effect could allow ground ice to develop and persist beyond the predicted boundaries of permafrost. Talus slopes are thus especially good examples of areas that present a landslide risk, and that large-scale/coarse-resolution permafrost model cannot address. Some studies investigated the present state of permafrost within talus slopes, based on ground temperature measurements, geophysical measurements, and fieldwork observations (e.g. Bommer et al., 2012; Delaloye and Lambiel, 2005; Germain and Milot,
2024; Jahn, 1984; Kenner et al., 2017; Kneisel et al., 2015; Lambiel and Pieracci, 2008; Millar et al., 2014; Morard et al., 2010; Otto and Sass, 2006; Phillips et al., 2009; Popescu et al., 2017; Scapozza et al., 2011; Staub et al., 2015). However, only a handful studies attempted to model numerically the thermal dynamics of talus slopes, at depth and on longer time scales (Myhra et al., 2019; Tanaka et al., 2006; Wicky and Hauck, 2017; Wicky et al., 2024; see Lan et al., 2024 for a review).

<figure><img src="image_3.png" /><figcaption>Figure 1: schematic diagram of the cross-section of a talus slope, and the internal air circulation caused by the chimney effect a. in winter ($T_{air} < T_{talus}$), and b. in summer ($T_{air} > T_{talus}$). Due to convective air circulation, in winter relatively hot air flows from the talus upward which sucks in relatively cold air at the bottom of the talus. In summer, relatively cold air within the talus flows downward. The zone of maximum overcooling (winter) and minimum warming (summer) is located approximately with the white dashed box. On both panels the arrow indicates the direction of air circulation; its colour indicates the relative temperature gradient within the talus.</figcaption></figure>

Therefore, in the present study, we present part of the Permolards project that focuses on the numerical modelling of the thermal dynamics of permafrost. More specifically, we investigate the Eyjafirði talus slope, the source area of the unpredicted Eyjafirði landslide (06/10/2020, Tröllaskagi peninsula, central northern Iceland), which threatened infrastructures and people. Notably, the presence of molards within the deposits of the Eyjafirði landslide indicates the presence of ground ice within the Eyjafirði talus slope at the date of the failure. We model heat transfer and phase changes of water using the finite element method (Courant, 1943) in two-dimensional (2D) cross-sectional models of the Eyjafirði talus slope, for the past ~ 20,000 years, with thermal constraints obtained from field temperature measurements.

We seek to answer the following questions:

1. Can we characterise the occurrence of a chimney effect at the Eyjafirði talus slope?

2. Can we successfully maintain permafrost within the Eyjafirði talus slope until the date of the failure?

3. What was the state of permafrost within the Eyjafirði talus slope at the date of the landslide?

4. What thermal properties of the Eyjafirði talus slope are required to maintain ice?

## 2. Study site

The Eyjafirði landslide (Figure 2a) occurred on October 6th 2020, on the flank of the Hleiðargarðsfjall mountain (65°24'23.60"N, 18°16'5.89"W; Tröllaskagi peninsula, central northern Iceland; Figure 3).
The detachment occurred in the Eyjafirði talus slope (~ 900 m. a.s.l.), within the Tertiary Basalt bedrock formation (Jóhannesson, 2014). The landslide travelled ~ 1,800 m horizontally, stopping less than 100 m from a farm - luckily causing no casualty and nearly no damage to infrastructures. Interestingly, news reports point out that no clear trigger factor can be identified, whether being heavy rainfall or seismic activity (see Frímann, 2020 and Petley, 2020 for news reports in English language).

<figure><img src="image_4.png" /><figcaption>Figure 2: the Eyjafirði landslide; a. photo of the Eyjafirði landslide, with trees for scale (taken on the day of the landslide, October 6th 2020; courtesy of J. P. Jónsson, who obtained it from the local police agents); b. molard example (approximately 2 m high) within the deposits of the Eyjafirði landslide (taken on July 2022).</figcaption></figure>
<figure><img src="image_5.png" /></figure>

Figure 3: location of the Eyjafirði landslide (white dot; 65°24′24.90″N, 18°16′16.05″O), that occurred on October 6th 2020 in the Tröllaskagi peninsula (central northern Iceland). The sub-panel shows the extension of the main panel (white rectangle) in Iceland. Background: Google Earth - Landsat/Copernicus.

The Eyjafirði talus slope is located outside of the predicted boundaries of permafrost (Czekirda et al., 2019; Etzelmüller et al., 2007, 2020; Lilleøren et al., 2013). However, molards located in the landslide deposits indicate that some of the mobilised material was ice-cemented (Figure 2b), suggesting that permafrost degradation was possibly involved in generating the Eyjafirði landslide (Morino et al., 2019). Additionally, we observed in the talus relict creep beads, indicative of the past presence of excess ice within the talus.

# 3. Data & Methods

In this section, we detail i. the acquisition and processing of data (i.e. 2D-geometry of the study site and temperature time series) that serve as inputs for the numerical models, and ii. our modelling approach.

## 3.1. Data

In order to perform our model runs, we have to i. characterise the 2D cross-sectional geometry of the Eyjafirði talus slope, and ii. reconstruct surface temperature time series to be used as forcing data.

### 3.1.1. Site geometry
Our numerical models require a 2D cross-sectional geometry to be functional; hence, we need i. a topographic profile of the talus slope prior to the landslide, and ii. an estimation of the bedrock topography below the talus.

To characterise the topography of the Eyjafirði talus slope, we acquired drone photos to generate a digital elevation model (DEM) and an orthophoto of the study site, with the Structure from Motion photogrammetry method (e.g. Fonstad et al., 2013; Westoby et al., 2012). The final DEM has a resolution of 32 cm/pixel. As the landslide was triggered in 2020, the DEM represents the post-landslide morphology of the Eyjafirði talus slope. However, we seek to model the thermal dynamics of the talus slope prior to the landslide. Therefore, we performed a reconstruction of the pre-landslide topography with the semi-automatic method described in Guimpier et al. (2022), using a GIS software (QGIS version 3.28.15). The reconstruction method summarises as follows (Figure 4):

1. Generation of contour lines at an interval of 2 metres;

2. Delineation of the talus slope, based on contour lines and orthophoto;

3. Cut-out of the contour lines based on the talus slope delineation (Figure 4a);

4. Generation of points every 2 metres along the remaining contour lines (Figure 4b);

5. Extraction of elevation values within each point;

6. Interpolation of the pre-landslide DEM from the points containing elevation data with the Triangulated Irregular Network method.

We then extract topographic data of the pre-landslide DEM along a selected profile parallel to the slope, located approximately in the middle of the talus slope, and crossing the detachment scarp of the landslide (Figure 4c, 4d).

To estimate the bedrock topography below the talus, we extract topographic data from the post-landslide DEM along the selected profile. At the bottom of the talus, we extrapolate upslope the average slope of the bedrock exposed by the landslide; at the top of the talus, we extrapolate downslope the average slope of the bedrock located right above the talus (Figure 4d). In both cases we consider approximately 15m (horizontally) of bedrock to calculate the average slopes. We consider that the bedrock topography below the talus follows the extrapolated profiles, and includes a slope break where the two extrapolations intersect (Figure 4d).
<figure><img src="image_6.png" /><figcaption>Figure 4: main steps of the pre-landslide topographic reconstruction; a. contour lines cut out according to the delineation of the talus slope; b. points generated along the contour lines; c. topographic profile (grey line) selected for defining the surface geometry of the talus slope. Points A and B located on the profile locate the profile shown in panel d; d. reconstructed pre-landslide topography (black line), measured post-landslide topography (grey line), and estimated bedrock topography (black dashed line) obtained by extrapolation of the bedrock topography around the talus.</figcaption></figure>

The extent of this profile is located on the full profile in panel c. In all panels, the white area represents the extent of the talus slope, the white crosses locate the thermal sensors (Section 3.1.2.), and north is up. White arrows indicate overall slope direction. Background (panels a., b., c.): orthophoto.
3.1.2. Temperature time series

In order to create temperature time series to be used as boundary conditions in our models, we i. measured temperature for ~ one year at the Eyjafirði talus slope; ii. compute correlations between our measured data and air temperature datasets that go further back to year 1881; and iii. compute modelled temperature time series back to – 20,000 years based on Etzelmüller et al. (2020).

We installed 18 thermal sensors in the Eyjafirði talus slope and in the surrounding rockwall, that measured temperature with an hourly time step between August 1st 2021 and July 9th 2022. We inserted rockwall sensors into 10 mm depth drilled holes in the bedrock, and we placed talus sensors at various depths (between the surface and into mini-boreholes up to 80 cm deep). Two sensors were placed at a given sensor location, to prevent data loss in case one sensor malfunctions. In some locations, we installed chains of sensors between the surface and the close subsurface. Therefore, the 18 sensors are spread over 6 sensor locations (Figure 5).

In the present study, we need to compute thermal gradients along the rockwall and talus (Section 3.2.3.); therefore, we use 4 of our 18 sensors only – two in the rockwall and two in the talus. We use in the rockwall one sensor above the talus and one below the talus (the ‘RW sensors’, named respectively ‘RW_top’ and ‘RWBot’; Geoprecision® M-Log5W-ROCK, accuracy ± 0.1°C, resolution 0.01°C, Figure 5). We use in the talus one sensor at its top and one at its bottom (the ‘talus sensors’, named respectively ‘Talus_top’ and ‘TalusBot’; model: Mouser® ibutton DS1925L-F5#, accuracy ± 0.5°C, resolution 0.065°C). Talus sensors were taped in order to increase their protection against potential damages and to waterproof them, with few insulating effects. Detailed settings for each sensor are summarised in Table 1.

<figure><img src="image_7.png" /><figcaption>Figure 5: location of the four used thermal sensors (RW_top, RWBot, Talus_top, TalusBot) on the Eyjafirði talus slope (white area). White dots indicate the location of other sensors. The white arrow indicates the direction of the slope. North is up. Background: orthophoto.</figcaption></figure>
<table><caption>Table 1: elevation, depth and geological setting of the emplacements of the four thermal sensors installed in the Eyjafirði talus slope and surrounding bedrock.</caption><thead><tr><th></th><th>Talus_top</th><th>Talus Bot</th><th>RW_top</th><th>RW bot</th></tr></thead><tbody><tr><td>elevation</td><td>887 m</td><td>791 m</td><td>895 m</td><td>821 m</td></tr><tr><td>depth</td><td>70 cm</td><td>45 cm</td><td>surface</td><td>surface</td></tr><tr><td>setting</td><td>coarse debris</td><td>coarse debris</td><td>vertical rockwall</td><td>vertical rockwall</td></tr></tbody></table>

From the measured data, we compute daily temperature averages for all four sensors. We also obtain air temperature time series from nearby meteorological stations from Veðurstofa Íslands, the Icelandic meteorological office. We then compute correlations (i.e. linear regression and R-squared value) between daily averages of both RW sensors, and daily averages from meteorological stations on the measurement period (August 1st 2021 to July 9th 2022). We select the two stations that show the best correlations (i.e. maximum R-squared values), here the automatic stations of Öxnadalsheiði (measurement period: July 29th 1995 to present) and Akureyri (January 1st 1949 to present). We then reconstruct daily temperatures of RW sensors back to January 1st 1949 (earliest available air temperature data), based on the equations of linear regression computed, for each sensor, between the measured data and the meteorological datasets. The Öxnadalsheiði station is preferred as it shows the best correlation with measured data; the Akureyri station is used when Öxnadalsheiði data are missing (e.g. prior to July 29th 1995). Statistics are summarised in Table 2.

We then compute monthly averages of RW sensors for the period January 1st 1949 to present. We obtain from Veðurstofa Íslands monthly air temperature time series from the Akureyri manned station (measurement period: November 1881 to December 2015). As above, we compute correlations (i.e. linear regression and R-squared values) between monthly average temperatures of both RW sensors, and monthly average temperatures from the Akureyri manned station, on the period January 1949 to December 2015 (in both cases, R-squared > 0.99). We then reconstruct monthly temperatures for RW sensors between November 1881 and December 1948, based on the equations of linear regressions.

<table><caption>Table 2: summary of statistics calculated for temperature reconstruction back to January 1st 1949, for all four sensors. For RW sensors, the n-factor is calculated with the Öxnadalsheiði station. For talus sensors, R-squared values are calculated based on the summer period; by comparison, we also provide the R-squared values for the winter period.</caption><thead><tr><th></th><th colspan="2">RW_top</th><th colspan="2">RWBot</th><th>Talus_top</th><th>Talus bot</th></tr></thead><tbody><tr><td>correlation with</td><td>Öxna.</td><td>Aku.</td><td>Öxna.</td><td>Aku.</td><td>RW_top</td><td>RW_top</td></tr><tr><td>R-squared</td><td>0.85</td><td>0.81</td><td>0.88</td><td>0.84</td><td>0.87</td><td>0.82</td></tr><tr><td>R-squared (w)</td><td>/</td><td>/</td><td>/</td><td>/</td><td>0.18</td><td>0.57</td></tr><tr><td>FDD</td><td>-764.332</td><td></td><td>-710.897</td><td></td><td>-75.049</td><td>-306.951</td></tr><tr><td>n-factor</td><td>0.98</td><td></td><td>0.91</td><td></td><td>0.10</td><td>0.40</td></tr></tbody></table>

Contrary to RW sensors, which are located on sub-vertical outcrops, talus sensors are affected by the winter snow cover. Therefore, air and ground surface temperature are decoupled during the snow cover period (Mote, 2008; Namias, 1985; Zhang, 2005). From data indicating the presence of a snow cover at the Akureyri automatic station (qualitatively between 0 = no snow and 4 = continuous snow cover), we calculate the average start and end date of the continuous snow cover based on the period 1998 to 2022. With this method, we approximate the 'summer' as a period with no snow cover, between May 7th and October 13th; and we approximate the 'winter' as a period with a continuous snow cover, between October 13th and May 7th (standard deviation: ± 17 days).
In the absence of local measurements of air temperature, we choose to reconstruct the talus temperature time series based on the RW reconstructed time series. We thus compute daily averages for both talus sensors, and compute correlations between their summer data and summer data from both RW sensors. RW_top showing the best correlations, we use its equations of linear regression to reconstruct temperature of talus sensors, daily (to January 1st 1949) and monthly (to November 1881). For winter periods, as talus sensors are affected by the snow cover, we multiply the results of the temperature reconstruction by their respective *n-factor* (Klene et al., 2001a, 2001b).

The *n-factor* is a quantitative representation of the insulating effect of the snow cover on a given sensor. It is calculated based on the freezing degree days (FDD), which are the sum of negative daily temperatures for a given sensor over the measurement period. In practise, the *n-factor* calculates as a ratio [FDD<sub>talus sensor</sub> / FDD<sub>RW sensor</sub>]. In our case, talus sensors are affected by the snow cover; hence, its insulating effect maintains temperatures around 0 °C, which reduces the number of days with negative temperatures, and ultimately reducing the FDD of talus sensors. On the contrary, RW sensors are not affected by the snow cover and serve as proxies for air temperature; hence, their FDD are relatively high, and equal to that of the local air temperature.

Therefore, the *n-factor* takes values between 0 (no negative daily temperatures measured at a talus sensor, i.e. perfect insulation from the snow cover) and 1 (equal negative daily temperatures measured at the talus and RW sensors, i.e. no insulating effect from the snow cover). In our case, for a given talus sensor, a value close to 0 indicates that temperatures reconstructed from RW_top are strongly attenuated, and thus closely maintain around 0 °C. On the contrary, a value close to 1 indicates that reconstructed temperatures are virtually not attenuated, and thus are nearly equal to temperatures of RW_top.

For all four sensors, we calculate the average temperature for the reference period 1961-1990. We then apply to each sensor the temperature linear deviation modelled in Etzelmüller et al. (2020), with the intermediate scenario of a Holocene Thermal Maximum (HTM) at + 0.5 °C. This linear deviation model was calculated specifically for Iceland (Etzelmüller et al., 2020), and is based on oxygen-isotope data from an ice core acquired in northern Greenland (North Greenland Ice Core Project members, 2004). The linear deviation model allows us to reconstruct, for each sensor, temperatures with a 10-years timestep back to the end of the Last Glacial Maximum (LGM) at -20,000 years, until -100 years before present.

No temperature data are available for the plateau upslope of the Eyjafirði talus slope. Nevertheless, the plateau is included in the selected topographic profile. Therefore, we create an artificial temperature time series for the plateau, assuming that the thermal regime of this part of the site doesn't influence the part of interest (talus). In summer the plateau time series takes the temperature values of the sensor with the highest elevation (i.e. RW_top). To account for the snow cover, in winter the plateau time series takes the temperature values of RW_top multiplied by the lowest *n-factor* we calculated in the talus (i.e. ~ 0.10, for Talus_top).

Ultimately, we create five time series covering the period from -20,000 years to December 13th 2022: RW_top, RWBot, Talus_top, TalusBot, and Plateau time series. A summary of the temporal resolution and sources of data for the temperature time series is provided in Figure 6. As the computation time step of the numerical simulations is automatic (Section 3.2.4.), FEFLOW interpolates linearly in-between temperature values of the time series when the data time step is larger than the computation time step.
<figure><img src="image_8.png" /><figcaption>Figure 6: schematic diagram of the resolution and sources of temperature data along the time series.</figcaption></figure>

As a comparison, we model temperatures at the TalusBot sensor over the year of measurement with the Cryogrid community model (Westermann et al., 2023; see also Aga et al., 2023). Cryogrid is a one-dimensional energy-balance model that uses downscaled topoclimatic data for simulating ground thermal dynamics and water/ice balance (Gisnås et al., 2013; Westermann et al., 2013; Westermann et al., 2016). As input air temperature data, we use ERA5 global reanalysis data (Hersbach et al., 2020) downscaled with the methods described in Filhol et al. (2023).

## 3.2. Modelling

We use the commercial software FEFLOW version 7.2 (for Finite Element subsurface FLOW and transport system; DHI-WASY GmbH; DHI, 2017; Diersch, 2014) to model the thermal dynamics of permafrost within the Eyjafirði talus slope. FEFLOW uses the finite element method (Courant, 1943; see also Argyris et al., 1979; Bathe, 2008; Holland, 1974; Strang, 1973; Williamson, 1980) to model heat transfer, mass transport and groundwater flow in various case scenarios (e.g. 2D or 3D, variably porous medium, saturated or not, transient or steady flow); its plug-in piFreeze (DHI, 2017) accounts for phase changes of water. Details on the mathematical approach of FEFLOW are given in Section 3.2.2. FEFLOW with piFreeze was used for less than a decade in geological studies, which investigated various scenarios: rockwalls in contexts of rock avalanche events (Cathala et al., 2024; Magnin and Josnin, 2021; Magnin et al., 2017), degrading permafrost plateau (Langford et al., 2020), thermal effect of hot buried pipeline in permafrost ground (Nagare et al., 2021), or heat transport within multi-layered compacted clay (Dalla Santa et al., 2022).

More specifically, the present study is a continuation of previous work done by Magnin et al. (2017) and Magnin and Josnin (2021), which used FEFLOW/piFreeze to investigate the relationship between permafrost dynamics and rockfall events in rockwall contexts within the Mont Blanc massif (French Alps). Notably, they found that the results of their numerical models match with borehole temperature and electrical resistivity data (Magnin et al., 2017) – validating the use of FEFLOW/piFreeze to study natural permafrost contexts. We innovatively adapt their methodology to investigate permafrost dynamics within a talus slope, considering conductive heat transfer only (i.e. no water flow, as in Magnin et al., 2017 and in one study case in Magnin and Josnin, 2021).

### 3.2.1. Site geometry and discretisation
We import in FEFLOW the topographic data (surface and underlying bedrock) described in Section 3.1.1. to create a 2D cross-sectional model of the Eyjafirði study site. Additionally, we create a 5 km-deep layer below the Eyjafirði model, following the methodology in Magnin et al. (2017) and Magnin and Josnin (2021). The geometry hence contains three polygons (Figure 7): talus, surface bedrock, and deep bedrock (the 5 km-deep layer). We mesh the geometry with the Triangulated Irregular Network (TIN) method, which discretises all polygons into triangular elements which have a node at each vertex. We use FEFLOW default parameters (notably maximum element angles = 30°), and refine the meshing at polygon edges in order to accommodate the computational complexity of the transition zone between talus and bedrock. The modelling spatial resolution required is, relatively speaking, high for the talus; intermediate for the surface bedrock; low for the deep bedrock. Therefore, we manually increase or decrease the relative meshing density for each polygon in order to match the required differential modelling resolution. Final relative meshing densities and number of elements per polygon are:

* deep bedrock: 1 (292 elements);

* surface bedrock: 25 (487 elements);

* talus: 2000 (222 elements).

Proceeding this way optimises the ratio [modelling resolution/processing time] of our numerical model runs.

<figure><img src="image_9.png" /><figcaption>Figure 7: geometry of the Eyjafirði study site, divided into three polygons (talus, surface bedrock, deep bedrock), each one discretised in triangular elements at various spatial resolutions. Thermal boundary conditions are indicated in italics. Arrows indicate the nodes where the base thermal boundary condition (geothermal heat flux of – 0.1 W.m⁻²) is applied. Crosses indicate the nodes where the sensor time series (RW_top, RWBot, Talus_top, TalusBot) are applied. The solid black line spans the nodes where the Plateau time series is applied. The dashed black lines span the nodes where values of sensor time series are extrapolated or interpolated.</figcaption></figure>

### 3.2.2. Numerical approach
We consider the Eyjafirði study site (bedrock and talus) as constituted entirely of basalt. Only porosity varies between the bedrock and talus. Bedrock and talus are individually considered homogeneous and isotropic. As we work in saturated conditions, the medium is composed of a solid (rock and ice) and a liquid phase (water). As in Magnin et al. (2017), we do not simulate water flow but heat transfer only. Therefore, heat transfer is by conduction only, and the energy conservation equation writes as follow (DHI, 2017; Diersch, 2014):

$$ \frac{\partial[\varphi\rho_L c_L + (1 - \varphi)\rho_S c_S]T}{\partial t} = -\nabla \cdot (\rho_L c_L q T - \Lambda \nabla T) \quad (1) $$

with $\varphi$ the porosity of the medium; $\rho_L c_L$ and $\rho_S c_S$ are the volumetric heat capacities of the liquid and solid phases respectively (J.m⁻³.K⁻¹, with $\rho$ the phase density and $c$ its specific heat capacity); $q$ is the flow velocity (set to 0 m.s⁻¹ in our study case, as we do not consider water flow); $T$ is the temperature (K); $\Lambda$ is the hydrodynamic thermal dispersion tensor (J.m⁻¹.s⁻¹.K⁻¹; which includes the thermal conductivity, W.m⁻¹.K⁻¹).

The piFreeze plugin then accounts for the presence of ice in the solid phase, by computing the ice bulk volumetric fraction as $\varepsilon_a + \varepsilon_w + \varepsilon_r + \varepsilon_i = 1$ (with $\varepsilon$ the bulk volumetric fractions of air, water, rock and ice respectively; as we work in saturated conditions, $\varepsilon_a = 0$ in our study case). Then, the freezing function (DHI, 2017) is computed as:

$$ F = \frac{\varepsilon_w \rho_w}{\varepsilon_w \rho_w + \varepsilon_i \rho_i} \quad (2) $$

with $\varepsilon$ and $\rho$ the volumetric bulk fractions and densities of the corresponding phase respectively (water or ice). This equation links the bulk mass of unfrozen water to the total bulk mass of water (frozen and unfrozen), and decreases with increasing presence of ice. Ice forms gradually within a temperature interval we define, which is $[T_0 - \frac{\Delta T}{2}, T_0 + \frac{\Delta T}{2}]$ with $T_0$ the freezing point of water (set at -1 °C; Magnin and Josnin, 2021; Magnin et al., 2017), and $\Delta T$ the length of the temperature interval (set at 2 °C; Magnin and Josnin, 2021; Magnin et al., 2017).

Finally, values of thermal conductivity and heat capacity of the solid phase in Equation 1 are modified by the ice bulk volumetric fraction $\varepsilon_i$ and the freezing function $F$, respectively as follows:

$$ \lambda_s = \lambda_{s,0} + \frac{\varepsilon_i(\lambda_i - \lambda_s)}{1 - \varepsilon} \quad (3) $$

and

$$ \rho_S c_S = \rho_{S,0} c_{S,0} + \frac{\rho_i c_i - \rho_S c_S}{1 - \varepsilon} - \left[ \frac{\rho_i L_f}{1 - \varepsilon} \cdot \frac{\partial \varepsilon_i}{\partial F} \cdot \frac{\partial F}{\partial T} \right] \quad (4) $$

with $L_f$ the latent heat of ice formation.

### 3.2.3. Thermal parameters, initial conditions and boundary conditions

The thermal parameters of the material the user can set in FEFLOW are: porosity, volumetric heat capacity of solid (i.e. the rock phase), thermal conductivity (TC) of solid, volumetric heat capacity of water, and TC of water. piFreeze adds the possibility to set: water density, ice density, volumetric heat capacity of ice, TC of ice, latent heat, and freezing function (freezing temperature, temperature range of the phase change, and residual liquid fraction). The values of thermal parameters we use are detailed in Table 3, with the associated reference(s) for each parameter value.

Notably, the thermal parameter that plays the most important role in the modelling of thermal dynamics is the TC of the rock phase. We perform model runs with values of 0.75, 1.1, and 1.75 W.m⁻¹.K⁻¹. Values
of ~ 1.75 W.m⁻¹.K⁻¹ are commonly found in the literature for Icelandic basalts (1.6 to 1.9 in Flóvenz and Saemundsson, 1993; 1.5 to 2.0 in Oxburgh and Agrell, 1982; 1.7 in Pálmason and Saemundsson, 1979). A more recent work suggests a value of ~ 1.1 W.m⁻¹.K⁻¹ (Ruether, 2011 - master's thesis), from 810 measurements made on Icelandic basalts. For reasons detailed in Section 5.1.2., we decided to perform model runs with a lower value than those typically found in the literature (i.e. 0.75 W.m⁻¹.K⁻¹). Nevertheless, the lowest TC value Ruether (2011) measured on Icelandic basalts is 0.62 W.m⁻¹.K⁻¹, which gives us credibility in using a TC value of 0.75 W.m⁻¹.K⁻¹.

For the talus, we use porosity values of 0.5 (as in Wicky and Hauck, 2017) and 0.8. We chose to investigate a porosity value as high as 0.8, as fieldwork observations revealed the possibility of excess ice within the talus (see Section 2.)

<table><caption>Table 3: summary of thermal parameters, initial conditions, and boundary conditions used to perform our model runs; associated reference is provided (when applicable). Letters in brackets indicate the phase (s solid, w liquid, i ice) the parameter applies to.</caption><thead><tr><th>parameter</th><th>value(s)</th><th>unit</th><th>reference(s)</th></tr></thead><tbody><tr><td>porosity</td><td>0.5 / 0.8 (talus)<br>0.05 (bedrock)</td><td>/</td><td>Wicky and Hauck (2017)</td></tr><tr><td>vol. heat cap. (s)</td><td>2.5</td><td>MJ.m<sup>-3</sup>.K<sup>-1</sup></td><td>Goranson (1942)<br>Laloui and Rotta Loria (2020)<br>Waples and Waples (2004)</td></tr><tr><td>TC (s)</td><td>0.75 / 1.1 / 1.75</td><td>W.m<sup>-1</sup>.K<sup>-1</sup></td><td>Flóvenz and Saemundsson (1993)<br>Oxburgh and Agrell (1982)<br>Pálmason and Saemundsson (1979)<br>Ruether (2011)</td></tr><tr><td>vol. heat cap. (w)</td><td>4.2</td><td>MJ.m<sup>-3</sup>.K<sup>-1</sup></td><td></td></tr><tr><td>TC (w)</td><td>0.65</td><td>W.m<sup>-1</sup>.K<sup>-1</sup></td><td></td></tr><tr><td>vol. heat cap. (i)</td><td>1.86</td><td>MJ.m<sup>-3</sup>.K<sup>-1</sup></td><td>same as</td></tr><tr><td>TC (i)</td><td>2.2</td><td>W.m<sup>-1</sup>.K<sup>-1</sup></td><td>Magnin et al. (2017) and</td></tr><tr><td>water density</td><td>1000</td><td>kg.m<sup>-3</sup></td><td>Magnin and Josnin (2021)</td></tr><tr><td>ice density</td><td>917.5</td><td>kg.m<sup>-3</sup></td><td></td></tr><tr><td>latent heat</td><td>0.334</td><td>MJ.kg<sup>-1</sup></td><td></td></tr><tr><td colspan="4" style="text-align: center;">Freezing function:</td></tr><tr><td>freezing temperature</td><td>-1</td><td>°C</td><td>same as</td></tr><tr><td>temperature range of<br>phase change</td><td>-2 to 0</td><td>°C</td><td>Magnin et al. (2017) and</td></tr><tr><td>residual liquid fraction</td><td>0.0001</td><td>/</td><td>Magnin and Josnin (2021)</td></tr><tr><td>temperature</td><td>-3 (all nodes)</td><td>°C</td><td>/</td></tr><tr><td>ice bulk volumetric<br>fraction</td><td>equal to porosity<br>(0.5 / 0.8)</td><td>/</td><td>/</td></tr><tr><td>surface temperatures</td><td>temperature time<br>series (Section 3.1.2.)</td><td>°C</td><td>/</td></tr><tr><td>geothermal heat flux</td><td>-0.1</td><td>W.m<sup>-2</sup></td><td>Jóhannesson et al. (2020),<br>from Hjartarson (2015)</td></tr></tbody></table>

All model runs have a temperature initial condition of – 3 °C applied at all nodes – i.e. homogeneous frozen conditions. Setting an initial temperature below the freezing temperature (– 1 °C) ensures that no
ice melting occurs in the talus during the first part of the runs (that serves as an initialisation of the models, Section 3.2.4.), when boundary conditions in surface temperatures are themselves several degrees below zero. Using frozen initial conditions is justified by the fact that our model runs start at the end of the LGM, before the northern hemisphere deglaciation (19,000 to 20,000 years before present; Clark et al., 2009). Depending on the model run, the initial ice bulk volumetric fraction within the talus varies between 0.5 (when porosity = 0.5) and 0.8 (when porosity = 0.8). Details on initial conditions are given in Table 3.

We use two types of temperature boundary conditions (Figure 7). First, we apply at all nodes at the base of the model a geothermal heat flux of –0.1 W.m⁻² (Jóhannesson et al., 2020, from Hjartarson, 2015). Second, we force temperatures at the surface of the models with the temperature time series we reconstructed (Section 3.1.2.; Figure 7). We apply the Plateau time series uniformly at all surface nodes on the plateau. We apply the RW_top time series at the first surface node above the talus; we apply the RWBot time series at the first surface node below the talus. We use the 1D linear interpolation FEFLOW tool to compute interpolated (between RW_top and RWBot nodes) and extrapolated time series (above RW_top and below RWBot) at all surface nodes, based on the values of RW_top and RWBot time series. We apply the Talus_top time series at the first surface node at the top of the talus; we apply the TalusBot time series at the first surface node at the bottom of the talus. We use the 1D linear interpolation to compute interpolated time series at all surface nodes on the talus. Details on boundary conditions are given in Table 3.

In the present study, we seek to investigate only the heat transfer within the Eyjafirði talus slope; hence we are not considering groundwater flow in our model runs. To ensure that meltwater from ice degradation does not flow, we apply at all nodes a hydraulic head of ~ 975 m (i.e. the highest elevation value of the topographic profile) and at all elements a hydraulic conductivity of 0 m.s⁻¹.

### 3.2.4. *Simulation process & scenarios*

We run numerical models on a computer with two processors (Intel® Xeon® CPU E5-2640 v4; 2.4 GHz) of 10 cores each, and 256 Go RAM. Each model runs from – 20,000 years to December 13th 2022, with an initial time step of one day, which adapts automatically to accommodate the simulation complexity that can arise during the model runs (e.g. when temperatures are around the freezing point and phase changes occur). The temporal resolution of temperature time series is shown in Figure 6. Each model run takes approximately 4 to 5 hours to complete. All simulation parameters (e.g. distribution of temperatures or ice bulk volumetric fraction) are automatically recorded at specific time steps we selected prior to all model runs.

The first part of the simulations (from – 20,000 to – 100 years) serves as an initialisation of the model, i.e. generates a plausible temperature zonation before simulating contemporary temperature distributions (further referred as the ‘LGM-initialisation’ method). In Magnin et al. (2017) and Magnin and Josnin (2021), the initialisation is performed individually by running models with the geothermal heat flux as the base boundary condition, and the first time steps of surface time series as the surface boundary conditions, until a thermal equilibrium is reached (further referred as the ‘equilibrium-initialisation’ method). These two methods giving similar results in our study case (Section 4.2), we choose to use the LGM-initialisation method as it reduces computation time.

Apart from a test run (Section 4.2.), we dedicate 6 simulation scenarios to the investigation of the sensitivity of our model to i. the initial porosity/ice content of the talus (Section 4.4.), and ii. the thermal conductivity (TC) of the rock phase (Section 4.5.). We investigate those two parameters as they have the greatest effect on the simulation results. Therefore, we perform:

i. three runs with a fixed TC (1.1 W.m⁻¹.K⁻¹) and an initial porosity/ice content varying between 0.3, 0.5 and 0.8;

ii. three runs with a fixed initial porosity/ice content (0.8) and a TC varying between 0.75, 1.1 and 1.75 W.m⁻¹.K⁻¹.
For each parameter (initial porosity/ice content and TC), we investigate values consistent with the literature, as well as extreme values that show the boundary behaviours of our model.

# 4. Results

## 4.1. Analysis of temperature data

Temperature measurements we acquired over ~ one year show that talus sensors are affected by the snow cover, which maintains temperatures around 0 °C at those sensors for part of the winter due to a zero-curtain effect (Outcalt et al., 1990; see also Hanson and Hoelzle, 2004; Figure 8a). However, in the middle of winter (i.e. between ~ November and April), both talus sensors show temperature oscillations – when they are supposed to be covered by the snow. Those oscillations reflect the temperature oscillations measured by RW_top (which serves as a proxy for air temperature), despite being significantly attenuated. Additionally, temperatures measured by the TalusBot sensor are also lower than the temperatures measured by the Talus_top sensor (Figure 8a).

Temperatures we modelled with Cryogrid adequately match measured temperatures in summer (Figure 8b). However, winter modelled temperatures do not show temperature oscillations that appear on measured data (Figure 8).

<figure><img src="image_10.png" /><figcaption>**Figure 8:** a. daily average temperatures, measured by RW_top (black line), Talus_top (red line) and TalusBot (green line) between August 1st 2021 and July 9th 2022. The approximate summer period for this year is greyed out. Example periods showing the effect of the snow cover and air coupling on temperatures measured by talus sensors are located; b. measured and Cryogrid-modelled temperatures at the TalusBot sensor. Modelled data do not show temperature oscillations due to the chimney effect in winter.</figcaption></figure>

The average temperatures for each sensor and over different periods (1961-1990, 1991-2020, year of measurement) is given in Table 4. We see from that table that i. in the talus, average temperatures were above 0 °C for both reference periods (1961-1990 and 1991-2020); and ii. despite being higher in altitude, temperatures reconstructed for Talus_top (887 m) are higher than temperatures reconstructed for TalusBot (791 m) for all three time periods: +0.73°C (1961-1990), +0.66 °C (1991-2020), +1.04 °C (year of measurement).
<table><caption>Table 4: average temperatures (in °C) over different time periods (1961-1990, 1991-2020, year of measurement) for all four sensors.</caption><thead><tr><th>period</th><th>RW_top</th><th>RWBot</th><th>Talus_top</th><th>TalusBot</th></tr></thead><tbody><tr><td>1961-1990</td><td>-0.30</td><td>0.20</td><td>1.21</td><td>0.48</td></tr><tr><td>1991-2020</td><td>0.88</td><td>1.47</td><td>1.69</td><td>1.03</td></tr><tr><td>08/2021-07/2022</td><td>1.03</td><td>1.64</td><td>2.21</td><td>1.17</td></tr></tbody></table>

## 4.2. Initialisation methods

We perform a test run to compare the two LGM-initialisation and the equilibrium-initialisation methods (Section 3.2.4.). These methods give similar final temperature zonation (Figure 9), justifying the use of the one requiring less computation time (i.e. LGM-initialisation method) for all further models presented in this study.

<figure><img src="image_11.png" /><figcaption>Figure 9: temperature zonation at the date of the Eyjafirði landslide (October 6th 2020), with a. the LGM-initialisation method, and b. the equilibrium-initialisation method, with a TC of the rock phase of 1.1 W.m⁻¹.K⁻¹ and an initial porosity/ice bulk volumetric fraction of 0.8.</figcaption></figure>

## 4.3. Ground temperature evolution over the Holocene

All model runs follow the same temperature dynamics (Figure 10): surface temperatures rise from –20,000 years to the HTM (between –10,000 to –8,000 years), leading to an onset of ground ice degradation within the talus during the HTM around –11,000 years (Figure 10a and b); ground ice entirely melts around –9,500 years. Temperatures then decrease from –8,000 years to contemporary times and stabilise slightly above 0 °C at the surface of the talus, maintaining its inside unfrozen (Figure 10c and d). The water within the talus then freezes again before the start of contemporary time series (i.e. November 1881). From the beginning of contemporary time series surface temperatures steadily increase, causing ground ice within the talus to degrade progressively (1950: Figure 12e and f), hence reaching a minimum at the date of the Eyjafirði landslide (Figure 12g and h).
<figure><img src="image_12.png" /></figure>

Figure 10: evolution of the temperature zonation and the ice bulk volumetric fraction over a model run, with a TC of the rock phase of 1.1 W.m⁻¹.K⁻¹ and an initial porosity/ice bulk volumetric fraction of 0.8. **a.** temperature zonation and **b.** ice bulk volumetric fraction at – 11,000 years; **c.** temperature
zonation and d. ice bulk volumetric fraction at – 8,000 years; e. temperature zonation and f. ice bulk volumetric fraction on 01/01/1950; g. temperature zonation and h. ice bulk volumetric fraction on the date of the Eyjafirði landslide (06/10/2020).

4.4. Sensitivity to porosity/ice content

Increasing the porosity and initial ice bulk volumetric fraction within the talus decreases the final temperatures at the base of the talus, leading to a larger zone where ground ice persists, with an overall higher ice bulk volumetric fraction (Figure 11).

<figure><img src="image_13.png" /></figure>

Figure 11: temperature zonation and ice bulk volumetric fraction at the date of the Eyjafirði landslide (October 6th 2020) with a TC of the rock phase of 1.1 W.m⁻¹.K⁻¹. a. temperature zonation and b. ice bulk volumetric fraction, for an initial porosity/ice bulk volumetric fraction of 0.3; c. temperature zonation and d. ice bulk volumetric fraction, for an initial porosity/ice bulk volumetric fraction of 0.5; e. temperature zonation and f. ice bulk volumetric fraction, for an initial porosity/ice bulk volumetric fraction of 0.8.
## 4.5. Sensitivity to the thermal conductivity

Decreasing the TC of the rock phase also decreases the final temperatures at the base of the talus, leading to a larger zone where ground ice persists (Figure 12). Notably, the core of persisting ground ice seems to shift toward the edge of the talus with decreasing TC of the rock phase.

<figure><img src="image_14.png" /></figure>

Figure 12: temperature zonation and ice bulk volumetric fraction at the date of the Eyjafirði landslide (October 6th 2020) with an initial porosity/ice bulk volumetric fraction of 0.8. a. temperature zonation and b. ice bulk volumetric fraction, for TC of the rock phase of 1.75 W.m⁻¹.K⁻¹; c. temperature zonation and d. ice bulk volumetric fraction, for TC of the rock phase of 1.1 W.m⁻¹.K⁻¹; e. temperature zonation and f. ice bulk volumetric fraction, for TC of the rock phase of 0.75 W.m⁻¹.K⁻¹.

# 5. Discussion
## 5.1. Thermal regime of the Eyjafirði talus slope

### 5.1.1. Persistence of permafrost within the Eyjafirði talus slope

The temperatures we measured at talus sensors over the period August 1st 2021 to July 9th 2022 show oscillations over a continuous period of time during winter (Figure 8), despite the presence of a snow cover. Such oscillations can be explained by the onset of an air circulation inside the talus slope - i.e. the occurrence of a chimney effect (Figure 1; Delaloye and Lambiel, 2005; Germain and Milot, 2024; Lambiel and Pieracci, 2008; Millar et al., 2014; Morard et al., 2010; Phillips et al., 2009; Popescu et al., 2017; Sawada et al., 2003; Wicky and Hauck, 2017; Wicky et al., 2024). In winter, air temperatures decrease below the inner talus temperatures, causing the relatively warm air within the talus to flow upward; by convection, relatively cold air from outside penetrates in the talus by its base, forming a zone of maximum overcooling. In summer, air temperatures increase above the inner talus temperatures, causing the relatively cold air within the talus to flow downward toward the base of the talus, forming a zone of minimum warming. Moreover, temperatures measured by the TalusBot sensor are generally lower than the temperatures measured by the TalusTop sensor, over the measurement period as well as on longer timescales (Table 4; Figure 8). Colder measured temperatures at the base of the talus also support the contemporary occurrence of a chimney effect within the Eyjafirði talus slope.

Our modelling approach is simplified, in the sense that we do not explicitly model the air convection (i.e. the chimney effect) within the talus. Nevertheless, we consider that our simulations model successfully the expected permafrost occurrence in the Eyjafirði talus slope. Indeed, the evolution of temperature zonation shows a zone with sub-zero temperatures that maintains perennially at depth, at the base of the talus slope (Figure 10). The direct consequence of such temperature zonation is that ground ice also maintains at the base of the Eyjafirði talus slope (Figure 10), as observed in previous studies (e.g. Delaloye and Lambiel, 2005; Lambiel and Pieracci, 2008; Morard et al., 2010; Phillips et al., 2009; Popescu et al., 2017), which is consistent with geomorphological observations of molards within the Eyjafirði deposits. Note that permafrost maintains although i. the Eyjafirði talus slope is located outside of the predicted boundaries of permafrost (Czekirda et al., 2019; Etzelmüller et al., 2007, 2020; Lilleøren et al., 2013), and ii. the average air temperature reconstructed at the RW bot sensor (proxy for the air temperature close to the talus slope) is positive for the reference periods 1961-1990 (+0.20 °C) and 1991-2020 (+1.47 °C). We suggest that reproducing the cooling effect of air convection is allowed in our simulations by the combination of our field temperature measurements (i.e. colder temperatures measured by TalusBot compared to TalusTop) and thermal parameters we use in the simulations (i.e. particularly low TC; Section 5.1.2.).

### 5.1.2. Thermal parameters of the Eyjafirði talus slope

The presence of molards in the deposits of the Eyjafirði landslide indicates that permafrost degradation was possibly among the trigger factors of the landslide (Section 2; Morino et al., 2019); hence, ground ice was present in the Eyjafirði talus slope at the time of the landslide. More specifically, Brideau et al. (2009) estimated that observed molard initial blocks are constituted of ~50% of ice. However, at the time of the Eyjafirði landslide, none of our model scenarios, even the most ice-conservative one (initial porosity/ice content = 0.8; TC = $0.75 W.m^{-1}.K^{-1}$; Figure 12e and f), is able to reach this 50% ice content. In addition, ground ice must be present relatively close to the surface of the talus slope to generate molards, as the uppermost parts of talus are primarily mobilised within landslides. However, in our models, ground ice maintains at the base of the talus only. Therefore, we can assume that, in the Eyjafirði talus slope at the time of the landslide, i. the actual temperature zonation shows a larger zone with sub-zero temperatures, and ii. the actual ground ice content is higher at the base of the talus and non-zero in the close subsurface - compared to our most ice-conservative scenario (initial porosity/ice content = 0.8; TC = $0.75 W.m^{-1}.K^{-1}$; Figure 12e and f).

We obtain this ice-conservative scenario with a TC of the rock phase of $0.75 W.m^{-1}.K^{-1}$, which is significantly lower than the TC values for Icelandic basalt usually found in the literature (~1.75 W.m$^{-1}.K^{-1}$ in Flóvenz and Saemundsson, 1993; Oxburgh and Agrell, 1982; Pálmason and Saemundsson,
1979). To our knowledge, only Ruether (2011) reported TC values as low as 0.75 W.m⁻¹.K⁻¹, but in less than 1% of their 810 TC measurements. Note also that, although the thermal effect of the chimney effect can be identified in the temperature records of our talus sensors (Figure 8; Section 5.1.1.), air convection itself is not simulated in our model runs. Wicky and Hauck (2017) and Wicky et al. (2024), who specifically modelled convective heat transfer within talus slopes, observed from sensitivity tests that air circulation between the inner talus and the outside plays a major role in (or even is necessary for) the persistence of permafrost in talus slopes. Air convection should theoretically enhance persistence of ground ice within the Eyjafirði talus slope, compared to our model results. Therefore, we assume that decreasing the rock-phase TC to a value considered as very low with respects to the literature (which seems to be necessary in our model runs for significant amounts of ground ice to persist closer to the surface; Figure 12) might have partly compensated for the effect of air convection, enhancing the persistence of ground ice within the talus.

## 5.2. Limitations of our model runs

We acknowledge several limitations in our modelling approach. First, the contemporary temperature time series (Section 3.1.2.) we use as boundary conditions in our models are not measured, but reconstructed from the air temperature measured by two meteorological stations (automatic stations in Öxnadalsheiði and Akureyri), both being located at ~ 30 km from the Eyjafirði study site. Hence, neither station represents perfectly the air temperature conditions in Eyjafirði, which is reflected in the imperfection of the correlations between temperatures measured at our sensors and at meteorological stations (R-squared values between 0.81 and 0.88; Table 2). Moreover, the Holocene temperature time series (i.e. from – 20,000 to present times) are inherently inaccurate: we modelled them based on a simple temperature shift, and at a 10-years time step only (Etzelmüller et al., 2020). Overall, our reconstructed time series represent inaccurately the actual surface temperature dynamics at the Eyjafirði talus slope. However, forcing data for such long-term numerical modelling studies necessarily have to be reconstructed or modelled, and thus are inherently subject to such inaccuracies.

Second, we do not take into account the dynamics of snow cover in our reconstruction method for temperature time series, as we consider the snow cover as a binary state only (present/not present). Such a method does not account for variations in snow depth, which lead to varying insulation effects depending on the snow depth (e.g. Ge and Gong, 2010; Luetschg et al., 2008; Park et al., 2015; Rödder and Kneisel, 2012; Zhang, 2005). The timing of snow cover in our temperature time series is also an approximation, as we calculate its average based on the period 1998-2022 and consider the resulting winter period as constant over the contemporary part of the temperature time series. This timing may not be representative of the whole time series. Nevertheless, different timings of snow cover imply different durations of its insulating effect, and can lead to ground cooling or warming depending on complex factors (e.g. Goodrich, 1982; Ling and Zhang, 2003; Stieglitz et al., 2003; Zhang, 2005). Additionally, the *n*-factors that are applied to every winter period all over the contemporary part of the time series are calculated based on the year of measurement only. However, we have no mean to know whether the calculated *n*-factors are representative of the average conditions at the Eyjafirði talus slope over the whole reconstruction period; reconstructed temperature data may then be under- or over- attenuated during winter periods. Note that the effect of snow cover insulation on surface temperatures are more adequately taken into account in energy-balance models (e.g. Cryogrid). However, we saw that such models cannot reproduce the temperature oscillations due to the onset of the chimney effect in winter (Figure 8b).

Third, we use a simple semi-automatic method to reconstruct the pre-landslide surface topography of the Eyjafirði talus slope (Section 3.1.1.; Guimpier et al., 2022). It was proved that this method gives results that are less accurate than other reconstruction methods (Guimpier et al., 2022), but also that it is supposed to be adapted to landslide with a detachment area that is relatively smooth. It should be the case with the Eyjafirði study site, as talus slopes are debris that lie at their angle of repose. For the underlying bedrock topography, we also use a simplistic estimation based on geometrical extrapolations. These topographic reconstructions may not represent perfectly the reality of the Eyjafirði talus slope. However, the methods we use have the advantage of being systematic, hence being easily reproducible
from one site to another, and reducing to a minimum biases arising from methods requiring more user intervention (e.g. the two other pre-landslide reconstruction methods shown in Guimpier et al., 2022).

Fourth, our numerical models do not represent the complexity of a natural talus slope. We consider a material that is homogeneous on the whole talus (e.g. in terms of porosity and initial ice bulk volumetric fraction), when natural talus slopes are constituted of debris of varying grain size, heterogeneously distributed within the talus (e.g. Curry and Black, 2003; Hinchliffe et al., 1998; Sanders, 2010). We also run our models in saturated conditions, and only heat transfer by conduction is considered (i.e. no water flow). However, the present study is innovative with respects to previous modelling studies (Myhra et al., 2019; Tanaka et al., 2006; Wicky and Hauck, 2017; Wicky et al., 2024) as i. our contemporary temperature time series are based on field measurements and have a relatively high temporal resolution (1 month to 1 day; Figure 6); and ii. the geometry of the talus is constrained by topographic observations and measurements. Therefore, considering a simplified geological setting allows us to evaluate the sensitivity of the model to other parameters (e.g. thermal conductivity; Section 5.1.2.).

Additionally, the fact that our model runs do not simulate air convection (i.e. the chimney effect) also represents a limitation: temperatures are certainly lower at the base of the talus compared to our model results, enhancing the persistence of permafrost (Section 5.1.2.).

# 6. Conclusions & Outlook

In the present study, we investigated through numerical modelling the permafrost dynamics at the Eyjafirði talus slope (Tröllaskagi peninsula, Iceland). This talus slope was the source area of a landslide which involved the degradation of azonal permafrost. We sought to understand the conditions under which permafrost could persist within the Eyjafirði talus slope. Our study leads to the following conclusions:

1. One-year near-surface temperature measurements allowed us to characterise the occurrence of a chimney effect in the Eyjafirði talus slope, and to gain relevant data to build forcing time series to model permafrost dynamics within this talus;

2. With a simplistic approach that does not explicitly model air convection, we successfully modelled the expected permafrost occurrence within the Eyjafirði talus slope. Notably, although the average air temperature over the reference period 1991-2020 is + 1.47 °C, our models account for the contemporary persistence of permafrost at the base of the talus;

3. Field observations of molards in the landslide deposits gave additional constraints on the ground ice content within the Eyjafirði talus slope at the time of the landslide. Therefore, we adopted an unconventional modelling approach where we sought for the combination of parameters that allow ground ice to persist;

4. The most ice-conservative model scenario (i.e. lowest rock phase TC) is required to only approach this constraint on ground ice content. However, all scenarios allow ground ice to persist within the Eyjafirði talus slope.

Our study also leads to several key point information for future work. First, thanks to our unconventional approach, we showed in this study that the thermal properties of the Eyjafirði talus slope (more specifically the thermal conductivity) that allow the persistence of permafrost differ significantly from the values that can typically be found in the literature. We highlight for future modelling studies that adopting such unconventional approaches, based on ground-truth evidence, can help in defining additional modelling scenarios with significantly different results. Overall, such approaches highlight the occurrence of physical phenomena other than those modelled: in our case, reducing the thermal conductivity may have partly compensated for the effect of air convection. Hence, future studies could investigate the sensibility of this assumption and compare the results of models with air convection vs. models without air convection and with reduced thermal conductivity. The results of such a study could make easier future numerical modelling studies of talus slopes, by providing a quantitative way to compensate for the effect of air convection.

Additionally, we provided all the necessary information to easily reproduce our workflow (reconstruction of pre-landslide topography, reconstruction of temperature time series, numerical
modelling). Therefore, similar studies could emerge on different study sites. Notably, Iceland has other talus-sourced landslides that bear evidence for the implication of permafrost degradation, and the Permolards project also aims to inventory molards around the globe – which will provide additional study site in periglacial regions worldwide.

Overall, the greatest challenge in realising more similar studies is logistics for temperature data acquisition. A way to overcome this challenge lies in energy-balance models (e.g. CryoGRID), which can produce ground temperature time series, downscaled from climatic data in order to match the scale of a site-specific study. However, we showed in the present study that energy-balance models are not suitable to study talus slopes specifically – and more generally sites where air convection influences the near-surface temperatures. We suggest that an interesting outlook would be to perform a comparative study between model runs forced with time series based on measured data and model runs forced with time series modelled by energy-balance models. Such a study would question the need for forcing field data in numerical modelling studies. Comparing the results of energy-balance models that include freeze/thaw processes (like CryoGRID) to the results of 2D heat transfer models (such as ours) for a given study site could shed light on the pros and cons of those two methods. Using ground-truth evidence for the presence of permafrost could highlight physical processes that are challenging to model for energy-balance models, and potentially ways to compensate it – similarly to our approach with the Eyjafirði talus slope.

## Acknowledgement

Authors acknowledge funding from the french Agence Nationale de la Recherche in the framework of the projects ANR-19-CE01-0010 PERMOLARDS and ANR-19-CE01-0018 WISPER. Authors are grateful towards Emmanuel Malet who helped with the retrieving of thermal sensors.

## References

Aga, J., Boike, J., Langer, M., Ingeman-Nielsen, T., Westermann, S., 2023. Simulating ice segregation and thaw consolidation in permafrost environments with the CryoGrid community model. The Cryosphere 17, 4179–4206. https://doi.org/10.5194/tc-17-4179-2023

Allen, S.K., Gruber, S., Owens, I.F., 2009. Exploring steep bedrock permafrost and its relationship with recent slope failures in the Southern Alps of New Zealand: Permafrost and Recent Slope Failures, New Zealand. Permafrost Periglac. Process. 20, 345–356. https://doi.org/10.1002/ppp.658

Argyris, J.H., Balmer, H., Doltsinis, J.St., Dunne, P.C., Haase, M., Kleiber, M., Malejannakis, G.A., Mlejnek, H.-P., Müller, M., Scharpf, D.W., 1979. Finite element method — the natural approach. Computer Methods in Applied Mechanics and Engineering 17–18, 1–106. https://doi.org/10.1016/0045-7825(79)90083-5

Bathe, K., 2008. Finite Element Method, in: Wah, B.W. (Ed.), Wiley Encyclopedia of Computer Science and Engineering. Wiley, pp. 1–12. https://doi.org/10.1002/9780470050118.ecse159

Boeckli, L., Brenning, A., Gruber, S., Noetzli, J., 2012a. A statistical approach to modelling permafrost distribution in the European Alps or similar mountain ranges. The Cryosphere 6, 125–140. https://doi.org/10.5194/tc-6-125-2012

Boeckli, L., Brenning, A., Gruber, S., Noetzli, J., 2012b. Permafrost distribution in the European Alps: calculation and evaluation of an index map and summary statistics. The Cryosphere 6, 807–820. https://doi.org/10.5194/tc-6-807-2012

Bommer, C., Fitze, P., Schneider, H., 2012. Thaw-Consolidation Effects on the Stability of Alpine Talus Slopes in Permafrost. Permafrost & Periglacial 23, 267–276. https://doi.org/10.1002/ppp.1751

Borgatti, L., Soldati, M., 2013. 7.30 Hillslope Processes and Climate Change, in: Treatise on Geomorphology. Elsevier, pp. 306–319. https://doi.org/10.1016/B978-0-12-374739-6.00180-9

Brideau, M., Stead, D., Lipovsky, P., Jaboyedoff, M., Hopkinson, C., Demuth, M., Barlow, J., Evans, S., Delaney, K., 2009. Preliminary description and slope stability analyses of the 2008 Little Salmon Lake and 2007 Mt. Steele landslides, Yukon. Yukon Exploration and Geology 119–134.
Cathala, M., Bock, J., Magnin, F., Ravanel, L., Ben Asher, M., Astrade, L., Bodin, X., Chambon, G., Deline, P., Faug, T., Genuite, K., Jaillet, S., Josnin, J., Revil, A., Richard, J., 2024. Predisposing, triggering and runout processes at a permafrost-affected rock avalanche site in the French Alps (Étache, June 2020). Earth Surf Processes Landf esp.5881.
https://doi.org/10.1002/esp.5881

Clark, P.U., Dyke, A.S., Shakun, J.D., Carlson, A.E., Clark, J., Wohlfarth, B., Mitrovica, J.X., Hostetler, S.W., McCabe, A.M., 2009. The Last Glacial Maximum. Science 325, 710–714.
https://doi.org/10.1126/science.1172873

Courant, R., 1943. Variational methods for the solution of problems of equilibrium and vibrations.
Bull. Amer. Math. Soc. 49, 1–23.

Curry, A.M., Black, R., 2003. Structure, sedimentology and evolution of rockfall talus, Mynydd Du, south Wales. Proceedings of the Geologists’ Association 114, 49–64.
https://doi.org/10.1016/S0016-7878(03)80027-5

Czekirda, J., Westermann, S., Etzelmüller, B., Jóhannesson, T., 2019. Transient Modelling of Permafrost Distribution in Iceland. Front. Earth Sci. 7, 130.
https://doi.org/10.3389/feart.2019.00130

Dagenais, S., Molson, J., Lemieux, J.-M., Fortier, R., Therrien, R., 2020. Coupled cryo-hydrogeological modelling of permafrost dynamics near Umiujaq (Nunavik, Canada).
Hydrogeol J 28, 887–904.
https://doi.org/10.1007/s10040-020-02111-3

Dalla Santa, G., Cola, S., Galgaro, A., 2022. Assessment of frost and heating penetration in compacted clay layers of landfill top covers in temperate climate. Bull Eng Geol Environ 81, 498.
https://doi.org/10.1007/s10064-022-02999-z

Delaloye, R., Lambiel, C., 2005. Evidence of winter ascending air circulation throughout talus slopes and rock glaciers situated in the lower belt of alpine discontinuous permafrost (Swiss Alps).
Norsk Geografisk Tidsskrift - Norwegian Journal of Geography 59, 194–203.
https://doi.org/10.1080/00291950510020673

Deluigi, N., Lambiel, C., Kanevski, M., 2017. Data-driven mapping of the potential mountain permafrost distribution. Science of The Total Environment 590–591, 370–380.
https://doi.org/10.1016/j.scitotenv.2017.02.041

DHI, 2017. FEFLOW 7.2 Documentation.

Diersch, H.-J.G., 2014. FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media. Springer Berlin Heidelberg, Berlin, Heidelberg.
https://doi.org/10.1007/978-3-642-38739-5

Dramis, F., Govi, M., Guglielmin, M., Mortara, G., 1995. Mountain permafrost and slope instability in the Italian Alps: The Val Pola Landslide. Permafrost and Periglacial Processes 6, 73–81.
https://doi.org/10.1002/ppp.3430060108

Etzelmüller, B., Farbrot, H., Guðmundsson, Á., Humlum, O., Tveito, O.E., Björnsson, H., 2007. The regional distribution of mountain permafrost in Iceland. Permafrost Periglac. Process. 18, 185–199.
https://doi.org/10.1002/ppp.583

Etzelmüller, B., Patton, H., Schomacker, A., Czekirda, J., Girod, L., Hubbard, A., Lilleøren, K.S., Westermann, S., 2020. Icelandic permafrost dynamics since the Last Glacial Maximum – model results and geomorphological implications. Quaternary Science Reviews 233, 106236.
https://doi.org/10.1016/j.quascirev.2020.106236

Fiddes, J., Aalstad, K., Lehning, M., 2022. TopoCLIM: rapid topography-based downscaling of regional climate model output in complex terrain v1.1. Geosci. Model Dev. 15, 1753–1768.
https://doi.org/10.5194/gmd-15-1753-2022

Fiddes, J., Gruber, S., 2014. TopoSCALE v.1.0: downscaling gridded climate data in complex terrain.
Geosci. Model Dev. 7, 387–405.
https://doi.org/10.5194/gmd-7-387-2014

Fiddes, J., Gruber, S., 2012. TopoSUB: a tool for efficient large area numerical modelling in complex topography at sub-grid scales. Geosci. Model Dev. 5, 1245–1257.
https://doi.org/10.5194/gmd-5-1245-2012

Filhol, S., Fiddes, J., Aalstad, K., 2023. TopoPyScale: A Python Package for Hillslope ClimateDownscaling. JOSS 8, 5059.
https://doi.org/10.21105/joss.05059

Fischer, L., Eisenbeiss, H., Kääb, A., Huggel, C., Haeberli, W., 2011. Monitoring topographic changes in a periglacial high-mountain face using high-resolution DTMs, Monte Rosa East Face,
Italian Alps: Topographic Investigation of Steep Periglacial Rock Walls. Permafrost Periglac. Process. 22, 140–152. https://doi.org/10.1002/ppp.717

Fischer, L., Kääb, A., Huggel, C., Noetzli, J., 2006. Geology, glacier retreat and permafrost degradation as controlling factors of slope instabilities in a high-mountain rock wall: the Monte Rosa east face. Nat. Hazards Earth Syst. Sci. 6, 761–772. https://doi.org/10.5194/nhess-6-761-2006

Flóvenz, Ó.G., Saemundsson, K., 1993. Heat flow and geothermal processes in Iceland. Tectonophysics 225, 123–138. https://doi.org/10.1016/0040-1951(93)90253-G

Fonstad, M.A., Dietrich, J.T., Courville, B.C., Jensen, J.L., Carbonneau, P.E., 2013. Topographic structure from motion: a new development in photogrammetric measurement: TOPOGRAPHIC STRUCTURE FROM MOTION. Earth Surf. Process. Landforms 38, 421–430. https://doi.org/10.1002/esp.3366

Frímann, J., 2020. Large landslide far south of Akureyri town. icelandgeology.net.

Ge, Y., Gong, G., 2010. Land surface insulation response to snow depth variability. J. Geophys. Res. 115, 2009JD012798. https://doi.org/10.1029/2009JD012798

Geertsema, M., Clague, J.J., Schwab, J.W., Evans, S.G., 2006. An overview of recent large catastrophic landslides in northern British Columbia, Canada. Engineering Geology 83, 120–143. https://doi.org/10.1016/j.enggeo.2005.06.028

Germain, D., Milot, J., 2024. An overcooled coarse-grained talus slope at low elevation: New insights on air circulation and environmental impacts, Cannon Cliff, New Hampshire, USA. Earth Surf Processes Landf 49, 1705–1720. https://doi.org/10.1002/esp.5792

Gisnås, K., Etzelmüller, B., Farbrot, H., Schuler, T.V., Westermann, S., 2013. CryoGRID 1.0: Permafrost Distribution in Norway estimated by a Spatial Numerical Model. Permafrost & Periglacial 24, 2–19. https://doi.org/10.1002/ppp.1765

Gisnås, K., Westermann, S., Schuler, T.V., Litherland, T., Isaksen, K., Boike, J., Etzelmüller, B., 2014. A statistical approach to represent small-scale variability of permafrost temperatures due to snow cover. The Cryosphere 8, 2063–2074. https://doi.org/10.5194/tc-8-2063-2014

Gisnås, K., Westermann, S., Schuler, T.V., Melvold, K., Etzelmüller, B., 2016. Small-scale variation of snow in a regional permafrost model. The Cryosphere 10, 1201–1215. https://doi.org/10.5194/tc-10-1201-2016

Goodrich, L.E., 1982. The influence of snow cover on the ground thermal regime. Can. Geotech. J. 19, 421–432. https://doi.org/10.1139/t82-047

Goranson, R.W., 1942. Heat Capacity; Heat of Fusion, in: Chairman, F.B., Schairer, J.F., Spicer, H.C. (Eds.), Handbook of Physical Constants. Geological Society of America, pp. 223–242. https://doi.org/10.1130/SPE36-p223

Gruber, S., 2012. Derivation and analysis of a high-resolution estimate of global permafrost zonation. The Cryosphere 6, 221–233. https://doi.org/10.5194/tc-6-221-2012

Gruber, S., Haeberli, W., 2007. Permafrost in steep bedrock slopes and its temperature-related destabilization following climate change. J. Geophys. Res. 112, F02S18. https://doi.org/10.1029/2006JF000547

Gruber, S., Hoelzle, M., 2001. Statistical modelling of mountain permafrost distribution: local calibration and incorporation of remotely sensed data. Permafrost & Periglacial 12, 69–77. https://doi.org/10.1002/ppp.374

Guglielmin, M., Aldighieri, B., Testa, B., 2003. PERMACLIM: a model for the distribution of mountain permafrost, based on climatic observations. Geomorphology 51, 245–257. https://doi.org/10.1016/S0169-555X(02)00221-0

Guimpier, A., Conway, S.J., Pajola, M., Lucchetti, A., Simioni, E., Re, C., Noblet, A., Mangold, N., Thomas, N., Cremonese, G., 2022. Pre-landslide topographic reconstruction in Baetis Chaos, mars using a CaSSIS Digital Elevation Model. Planetary and Space Science 218, 105505. https://doi.org/10.1016/j.pss.2022.105505

Haeberli, W., Noetzli, J., Arenson, L., Delaloye, R., Gärtner-Roer, I., Gruber, S., Isaksen, K., Kneisel, C., Krautblatter, M., Phillips, M., 2010. Mountain permafrost: development and challenges of a young research field. J. Glaciol. 56, 1043–1058. https://doi.org/10.3189/002214311796406121
Haeberli, W., Schaub, Y., Huggel, C., 2017. Increasing risks related to landslides from degrading permafrost into new lakes in de-glaciating mountain ranges. Geomorphology 293, 405–417.
https://doi.org/10.1016/j.geomorph.2016.02.009

Hanson, S., Hoelzle, M., 2004. The thermal regime of the active layer at the Murtèl rock glacier based on data from 2002. Permafrost & Periglacial 15, 273–282. https://doi.org/10.1002/ppp.499

Harris, C., 2005. Climate Change, Mountain Permafrost Degradation and Geotechnical Hazard, in: Huber, U.M., Bugmann, H.K.M., Reasoner, M.A. (Eds.), Global Change and Mountain Regions, Advances in Global Change Research. Springer Netherlands, Dordrecht, pp. 215–224. https://doi.org/10.1007/1-4020-3508-X_22

Harris, C., Vonder Muhll, D., 2001. Permafrost and Climate in Europe. Climate Change, Mountain Permafrost Degradation and Geotechnical Hazard, in: Visconti, G., Beniston, M., Iannorelli, E.D., Barba, D. (Eds.), Global Change and Protected Areas, Advances in Global Change Research. Springer Netherlands, Dordrecht, pp. 71–82. https://doi.org/10.1007/0-306-48051-4_9

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., De Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., Thépaut, J., 2020. The ERA5 global reanalysis. Quart J Royal Meteoro Soc 146, 1999–2049.
https://doi.org/10.1002/qj.3803

Hinchliffe, S., Ballantyne, C.K., Walden, J., 1998. The structure and sedimentology of relict talus, Trotternish, northern Skye, Scotland. Earth Surf. Process. Landforms 23, 545–560.
https://doi.org/10.1002/(SICI)1096-9837(199806)23:6<545::AID-ESP876>3.0.CO;2-E

Hjartarson, Á., 2015. Heat flow in Iceland, in: Proceedings of the World Geothermal Congress. Presented at the Worl Geothermal Congress, Melbourne, Australia.

Hjort, J., Streletskiy, D., Doré, G., Wu, Q., Bjella, K., Luoto, M., 2022. Impacts of permafrost degradation on infrastructure. Nat Rev Earth Environ 3, 24–38.
https://doi.org/10.1038/s43017-021-00247-8

Holland, I., 1974. Fundamentals of the finite element method. Computers & Structures 4, 3–15.
https://doi.org/10.1016/0045-7949(74)90074-1

Jafarov, E.E., Marchenko, S.S., Romanovsky, V.E., 2012. Numerical modeling of permafrost dynamics in Alaska using a high spatial resolution dataset. The Cryosphere 6, 613–624.
https://doi.org/10.5194/tc-6-613-2012

Jahn, A., 1984. Periglacial talus slopes. Geomorphological studies on spitsbergen and in Northern Scandinavia. Polar Geography and Geology 8, 177–193.
https://doi.org/10.1080/10889378409377224

Jóhannesson, H., 2014. Geological map of Iceland at 1: 600 000 (bedrock geology).

Jóhannesson, T., Pálmason, B., Hjartarson, Á., Jarosch, A.H., Magnússon, E., Belart, J.M.C., Gudmundsson, M.T., 2020. Non-surface mass balance of glaciers in Iceland. J. Glaciol. 66, 685–697. https://doi.org/10.1017/jog.2020.37

Karjalainen, O., Aalto, J., Luoto, M., Westermann, S., Romanovsky, V.E., Nelson, F.E., Etzelmüller, B., Hjort, J., 2019. Circumpolar permafrost maps and geohazard indices for near-future infrastructure risk assessments. Sci Data 6, 190037. https://doi.org/10.1038/sdata.2019.37

Kellerer-pirklbauer, A., Lieb, G.K., Avian, M., Carrivick, J., 2012. Climate change and rock fall events in high mountain areas: numerous and extensive rock falls in 2007 at mittlerer burgstall, central austria. Geografiska Annaler: Series A, Physical Geography 94, 59–78.
https://doi.org/10.1111/j.1468-0459.2011.00449.x

Kenner, R., Phillips, M., Hauck, C., Hilbich, C., Mulsow, C., Bühler, Y., Stoffel, A., Buchroithner, M., 2017. New insights on permafrost genesis and conservation in talus slopes based on observations at Flüelapass, Eastern Switzerland. Geomorphology 290, 101–113.
https://doi.org/10.1016/j.geomorph.2017.04.011

Klene, A.E., Nelson, F.E., Shiklomanov, N.I., 2001a. THE N-FACTOR AS A TOOL IN GEOCRYOLOGICAL MAPPING: SEASONAL THAW IN THE KUPARUK RIVER
BASIN, ALASKA. Physical Geography 22, 449–466.
https://doi.org/10.1080/02723646.2001.10642754

Klene, A.E., Nelson, F.E., Shiklomanov, N.I., Hinkel, K.M., 2001b. The N-factor in Natural Landscapes: Variability of Air and Soil-Surface Temperatures, Kuparuk River Basin, Alaska, U.S.A. Arctic, Antarctic, and Alpine Research 33, 140–148.
https://doi.org/10.1080/15230430.2001.12003416

Kneisel, C., Emmert, A., Polich, P., Zollinger, B., Egli, M., 2015. Soil geomorphology and frozen ground conditions at a subalpine talus slope having permafrost in the eastern Swiss Alps. CATENA 133, 107–118. https://doi.org/10.1016/j.catena.2015.05.005

Laloui, L., Rotta Loria, A.F., 2020. Analysis and Design of Energy Geostructures. Elsevier.
https://doi.org/10.1016/C2017-0-04166-6

Lambiel, C., Pieracci, K., 2008. Permafrost distribution in talus slopes located within the alpine periglacial belt, Swiss Alps. Permafrost & Periglacial 19, 293–304.
https://doi.org/10.1002/ppp.624

Lan, S., Cao, B., Hu, Y., Sun, Z., Ma, R., Li, X., 2024. Talus and its cooling effects on the thermal regime of permafrost: A review. Permafrost & Periglacial 35, 60–75.
https://doi.org/10.1002/ppp.2213

Langford, J.E., Schincariol, R.A., Nagare, R.M., Quinton, W.L., Mohammed, A.A., 2020. Transient and Transition Factors in Modeling Permafrost Thaw and Groundwater Flow. Groundwater 58, 258–268. https://doi.org/10.1111/gwat.12903

Levavasseur, G., Vrac, M., Roche, D.M., Paillard, D., Martin, A., Vandenberghe, J., 2011. Present and LGM permafrost from climate simulations: contribution of statistical downscaling. Clim. Past 7, 1225–1246. https://doi.org/10.5194/cp-7-1225-2011

Lilleøren, K.S., Etzelmüller, B., Gärtner-Roer, I., Kääb, A., Westermann, S., Guðmundsson, Á., 2013. The Distribution, Thermal Characteristics and Dynamics of Permafrost in Tröllaskagi, Northern Iceland, as Inferred from the Distribution of Rock Glaciers and Ice-Cored Moraines. Permafrost & Periglacial 24, 322–335. https://doi.org/10.1002/ppp.1792

Ling, F., Zhang, T., 2003. Impact of the timing and duration of seasonal snow cover on the active layer and permafrost in the Alaskan Arctic. Permafrost & Periglacial 14, 141–150.
https://doi.org/10.1002/ppp.445

Luetschg, M., Lehning, M., Haeberli, W., 2008. A sensitivity study of factors influencing warm/thin permafrost in the Swiss Alps. J. Glaciol. 54, 696–704.
https://doi.org/10.3189/002214308786570881

Lyle, R.R., Hutchinson, D.J., Preston, Y., 2004. Landslide processes in discontinuous permafrost, Little Salmon Lake (NTS 105L/1 and 2), south-central Yukon. Yukon Exploration and Geology 193–204.

Magnin, F., Brenning, A., Bodin, X., Deline, P., Ravanel, L., 2015. Modélisation statistique de la distribution du permafrost de paroi : application au massif du Mont Blanc. geomorphologie 21, 145–162. https://doi.org/10.4000/geomorphologie.10965

Magnin, F., Etzelmüller, B., Westermann, S., Isaksen, K., Hilger, P., Hermanns, R.L., 2019. Permafrost distribution in steep rock slopes in Norway: measurements, statistical modelling and implications for geomorphological processes. Earth Surf. Dynam. 7, 1019–1040.
https://doi.org/10.5194/esurf-7-1019-2019

Magnin, F., Josnin, J.-Y., 2021. Water Flows in Rock Wall Permafrost: A Numerical Approach Coupling Hydrological and Thermal Processes. JGR Earth Surface 126, e2021JF006394.
https://doi.org/10.1029/2021JF006394

Magnin, F., Josnin, J.-Y., Ravanel, L., Pergaud, J., Pohl, B., Deline, P., 2017. Modelling rock wall permafrost degradation in the Mont Blanc massif from the LIA to the end of the 21st century. The Cryosphere 11, 1813–1834. https://doi.org/10.5194/tc-11-1813-2017

Magnin, F., Ravanel, L., Ben-Asher, M., Bock, J., Cathala, M., Duvillard, P.-A., Jean, P., Josnin, J.-Y., Kaushik, S., Revil, A., Deline, P., 2023. From Rockfall Observation to Operational Solutions: Nearly 20 years of Cryo-gravitational Hazard Studies in Mont-Blanc Massif. rga.
https://doi.org/10.4000/rga.11703
Milana, J.P., 2016. Molards and Their Relation to Landslides Involving Permafrost Failure: Molards and permafrost degradation. Permafrost and Periglac. Process. 27, 271–284.
https://doi.org/10.1002/ppp.1878

Millar, C.I., D. Westfall, R., Delany, D.L., 2014. Thermal Regimes and Snowpack Relations of Periglacial Talus Slopes, Sierra Nevada, California, U.S.A. Arctic, Antarctic, and Alpine Research 46, 483–504. https://doi.org/10.1657/1938-4246-46.2.483

Morard, S., Delaloye, R., Lambiel, C., 2010. Pluriannual thermal behavior of low elevation cold talus slopes in western Switzerland. Geogr. Helv. 65, 124–134. https://doi.org/10.5194/gh-65-124-2010

Morino, C., Conway, S.J., Balme, M.R., Helgason, J.K., Sæmundsson, Þ., Jordan, C., Hillier, J., Argles, T., 2021. The impact of ground-ice thaw on landslide geomorphology and dynamics: two case studies in northern Iceland. Landslides 18, 2785–2812.
https://doi.org/10.1007/s10346-021-01661-1

Morino, C., Conway, S.J., Sæmundsson, Þ., Helgason, J.K., Hillier, J., Butcher, F.E.G., Balme, M.R., Jordan, C., Argles, T., 2019. Molards as an indicator of permafrost degradation and landslide processes. Earth and Planetary Science Letters 516, 136–147.
https://doi.org/10.1016/j.epsl.2019.03.040

Mote, T.L., 2008. On the Role of Snow Cover in Depressing Air Temperature. Journal of Applied Meteorology and Climatology 47, 2008–2022. https://doi.org/10.1175/2007JAMC1823.1

Myhra, K.S., Westermann, S., Etzelmüller, B., 2019. Modeling Conductive Heat Flow Between Steep Rock Walls and Talus Slopes – Thermal Processes and Geomorphological Implications. Front. Earth Sci. 7, 192. https://doi.org/10.3389/feart.2019.00192

Nagare, R.M., Mohammed, A.A., Park, Y.-J., Schincariol, R.A., 2021. Modeling shallow ground temperatures around hot buried pipelines in cold regions. Cold Regions Science and Technology 187, 103295. https://doi.org/10.1016/j.coldregions.2021.103295

Namias, J., 1985. Some Empirical Evidence for the Influence of Snow Cover on Temperature and Precipitation. Mon. Wea. Rev. 113, 1542–1553. https://doi.org/10.1175/1520-0493(1985)113<1542:SEEFTI>2.0.CO;2

Nelson, F.E., 1986. Permafrost Distribution in Central Canada: Applications of a Climate-Based Predictive Model. Annals of the Association of American Geographers 76, 550–569.
https://doi.org/10.1111/j.1467-8306.1986.tb00136.x

North Greenland Ice Core Project members, 2004. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature 431, 147–151.
https://doi.org/10.1038/nature02805

Obu, J., Westermann, S., Bartsch, A., Berdnikov, N., Christiansen, H.H., Dashtseren, A., Delaloye, R., Elberling, B., Etzelmüller, B., Kholodov, A., Khomutov, A., Kääb, A., Leibman, M.O., Lewkowicz, A.G., Panda, S.K., Romanovsky, V., Way, R.G., Westergaard-Nielsen, A., Wu, T., Yamkhin, J., Zou, D., 2019. Northern Hemisphere permafrost map based on TTOP modelling for 2000–2016 at 1 km2 scale. Earth-Science Reviews 193, 299–316.
https://doi.org/10.1016/j.earscirev.2019.04.023

Otto, J.C., Sass, O., 2006. Comparing geophysical methods for talus slope investigations in the Turtmann valley (Swiss Alps). Geomorphology 76, 257–272.
https://doi.org/10.1016/j.geomorph.2005.11.008

Outcalt, S.I., Nelson, F.E., Hinkel, K.M., 1990. The zero-curtain effect: Heat and mass transfer across an isothermal region in freezing soil. Water Resour. Res. 26, 1509–1516.
https://doi.org/10.1029/WR026i007p01509

Oxburgh, E.R., Agrell, S.O., 1982. Thermal conductivity and temperature structure of the Reydarfjordur borehole. J. Geophys. Res. 87, 6423–6428.
https://doi.org/10.1029/JB087iB08p06423

Pálmason, G., Sæmundsson, K., 1979. Summary of Conductive Heat Flow in Iceland, in: Čermák, V., Rybach, L. (Eds.), Terrestrial Heat Flow in Europe. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 218–220. https://doi.org/10.1007/978-3-642-95357-6_22

Park, H., Fedorov, A.N., Zheleznyak, M.N., Konstantinov, P.Y., Walsh, J.E., 2015. Effect of snow cover on pan-Arctic permafrost thermal regimes. Clim Dyn 44, 2873–2895.
https://doi.org/10.1007/s00382-014-2356-5
Patton, A.I., Rathburn, S.L., Capps, D.M., 2019. Landslide response to climate change in permafrost regions. *Geomorphology* 340, 116–128. https://doi.org/10.1016/j geomorph.2019.04.029

Penna, I.M., Magnin, F., Nicolet, P., Etzelmüller, B., Hermanns, R.L., Böhme, M., Kristensen, L., Nöel, F., Bredal, M., Dehls, J.F., 2023. Permafrost controls the displacement rates of large unstable rock-slopes in subarctic environments. *Global and Planetary Change* 220, 104017. https://doi.org/10.1016/j.gloplacha.2022.104017

Petley, D., 2020. Eyjafordur (Eyjafirði): a large landslide in Iceland yesterday. AGU Blogosphere.

Petley, D., 2012. Global patterns of loss of life from landslides. *Geology* 40, 927–930. https://doi.org/10.1130/G33217.1

Phillips, M., Mutter, E.Z., Kern-Luetschg, M., Lehning, M., 2009. Rapid degradation of ground ice in a ventilated talus slope: Flüela Pass, Swiss Alps. *Permafrost & Periglacial* 20, 1–14. https://doi.org/10.1002/ppp.638

Popescu, R., Vespremeanu-Stroe, A., Onaca, A., Vasile, M., Cruceru, N., Pop, O., 2017. Low-altitude permafrost research in an overcooled talus slope–rock glacier system in the Romanian Carpathians (Detunata Goală, Apuseni Mountains). *Geomorphology* 295, 840–854. https://doi.org/10.1016/j geomorph.2017.07.029

Ravanel, L., Magnin, F., Deline, P., 2017. Impacts of the 2003 and 2015 summer heatwaves on permafrost-affected rock-walls in the Mont Blanc massif. *Science of The Total Environment* 609, 132–143. https://doi.org/10.1016/j.scitotenv.2017.07.055

Rödder, T., Kneisel, C., 2012. Influence of snow cover and grain size on the ground thermal regime in the discontinuous permafrost zone, Swiss Alps. *Geomorphology* 175–176, 176–189. https://doi.org/10.1016/j geomorph.2012.07.008

Ruether, J., 2011. The validation of the LMC device: analysis of Icelandic basaltic rocks.

Sæmundsson, Þ., Morino, C., Helgason, J.K., Conway, S.J., Pétursson, H.G., 2018. The triggering factors of the Móafellshyrna debris slide in northern Iceland: Intense precipitation, earthquake activity and thawing of mountain permafrost. *Science of The Total Environment* 621, 1163–1175. https://doi.org/10.1016/j.scitotenv.2017.10.111

Sanders, D., 2010. Sedimentary facies and progradational style of a Pleistocene talus-slope succession, Northern Calcareous Alps, Austria. *Sedimentary Geology* 228, 271–283. https://doi.org/10.1016/j.sedgeo.2010.05.002

Sawada, Y., Ishikawa, M., Ono, Y., 2003. Thermal regime of sporadic permafrost in a block slope on Mt. Nishi-Nupukaushinupuri, Hokkaido Island, Northern Japan. *Geomorphology* 52, 121–130. https://doi.org/10.1016/S0169-555X(02)00252-0

Scapozza, C., Lambiel, C., Baron, L., Marescot, L., Reynard, E., 2011. Internal structure and permafrost distribution in two alpine periglacial talus slopes, Valais, Swiss Alps. *Geomorphology* 132, 208–221. https://doi.org/10.1016/j geomorph.2011.05.010

Smith, M., 1990. Potential Responses of Permafrost to Climatic Change. J. Cold Reg. Eng. 4, 29–37. https://doi.org/10.1061/(ASCE)0887-381X(1990)4:1(29)

Smith, M.W., Riseborough, D.W., 1996. Permafrost monitoring and detection of climate change. Permafrost Periglac. Process. 7, 301–309. https://doi.org/10.1002/(SICI)1099-1530(199610)7:4<301::AID-PPP231>3.0.CO;2-R

Staub, B., Marmy, A., Hauck, C., Hilbich, C., Delaloye, R., 2015. Ground temperature variations in a talus slope influenced by permafrost: a comparison of field observations and model simulations. *Geogr. Helv.* 70, 45–62. https://doi.org/10.5194/gh-70-45-2015

Stendel, M., Romanovsky, V.E., Christensen, J.H., Sazonova, T., 2007. Using dynamical downscaling to close the gap between global change scenarios and local permafrost dynamics. *Global and Planetary Change* 56, 203–214. https://doi.org/10.1016/j.gloplacha.2006.07.014

Stieglitz, M., Déry, S.J., Romanovsky, V.E., Osterkamp, T.E., 2003. The role of snow cover in the warming of arctic permafrost. *Geophysical Research Letters* 30, 2003GL017337. https://doi.org/10.1029/2003GL017337

Stocker-Mittaz, C., Hoelzle, M., Haeberli, W., 2002. Modelling alpine permafrost distribution based on energy-balance data: a first step. *Permafrost & Periglacial* 13, 271–282. https://doi.org/10.1002/ppp.426

Strang, G., 1973. PIECEWISE POLYNOMIALS AND THE FINITE ELEMENT METHOD. Bull. Amer. Math. Soc. 79.
Tanaka, H.L., Nohara, D., Byun, H.-R., 2006. Numerical Simulation of Wind Hole Circulation at Ice Valley in Korea Using a Simple 2D Model. Journal of the Meteorological Society of Japan 84, 1073–1084. https://doi.org/10.2151/jmsj.84.1073

Wang, T., Wu, T., Wang, P., Li, R., Xie, C., Zou, D., 2019. Spatial distribution and changes of permafrost on the Qinghai-Tibet Plateau revealed by statistical models during the period of 1980 to 2010. Science of The Total Environment 650, 661–670. https://doi.org/10.1016/j.scitotenv.2018.08.398

Waples, D.W., Waples, J.S., 2004. A Review and Evaluation of Specific Heat Capacities of Rocks, Minerals, and Subsurface Fluids. Part 1: Minerals and Nonporous Rocks. Natural Resources Research 13, 97–122. https://doi.org/10.1023/B:NARR.0000032647.41046.e7

Westermann, S., Ingeman-Nielsen, T., Scheer, J., Aalstad, K., Aga, J., Chaudhary, N., Etzelmüller, B., Filhol, S., Kääb, A., Renette, C., Schmidt, L.S., Schuler, T.V., Zweigel, R.B., Martin, L., Morard, S., Ben-Asher, M., Angelopoulos, M., Boike, J., Groenke, B., Miesner, F., Nitzbon, J., Overduin, P., Stuenzi, S.M., Langer, M., 2023. The CryoGrid community model (version 1.0) – a multi-physics toolbox for climate-driven simulations in the terrestrial cryosphere. Geosci. Model Dev. 16, 2607–2647. https://doi.org/10.5194/gmd-16-2607-2023

Westermann, S., Langer, M., Boike, J., Heikenfeld, M., Peter, M., Etzelmüller, B., Krinner, G., 2016. Simulating the thermal regime and thaw processes of ice-rich permafrost ground with the land-surface model CryoGrid 3. Geosci. Model Dev. 9, 523–546. https://doi.org/10.5194/gmd-9-523-2016

Westermann, S., Schuler, T.V., Gisnås, K., Etzelmüller, B., 2013. Transient thermal modeling of permafrost conditions in Southern Norway. The Cryosphere 7, 719–739. https://doi.org/10.5194/tc-7-719-2013

Westoby, M.J., Brasington, J., Glasser, N.F., Hambrey, M.J., Reynolds, J.M., 2012. ‘Structure-from-Motion’ photogrammetry: A low-cost, effective tool for geoscience applications. Geomorphology 179, 300–314. https://doi.org/10.1016/j.geomorph.2012.08.021

Wicky, J., Hauck, C., 2017. Numerical modelling of convective heat transport by air flow in permafrost talus slopes. The Cryosphere 11, 1311–1325. https://doi.org/10.5194/tc-11-1311-2017

Wicky, J., Hilbich, C., Delaloye, R., Hauck, C., 2024. Modeling the Link Between Air Convection and the Occurrence of Short-Term Permafrost in a Low-Altitude Cold Talus Slope. Permafrost & Periglacial 35, 202–217. https://doi.org/10.1002/ppp.2224

Williamson, F., 1980. Richard courant and the finite element method: A further look. Historia Mathematica 7, 369–378. https://doi.org/10.1016/0315-0860(80)90001-4

Yin, G.-A., Niu, F.-J., Lin, Z.-J., Luo, J., Liu, M.-H., 2021. Data-driven spatiotemporal projections of shallow permafrost based on CMIP6 across the Qinghai–Tibet Plateau at 1 km2 scale. Advances in Climate Change Research 12, 814–827. https://doi.org/10.1016/j.accre.2021.08.009

Zhang, T., 2005. Influence of the seasonal snow cover on the ground thermal regime: An overview. Reviews of Geophysics 43, 2004RG000157. https://doi.org/10.1029/2004RG000157

Zhang, Y., Qian, B., Hong, G., 2020. A Long-Term, 1-km Resolution Daily Meteorological Dataset for Modeling and Mapping Permafrost in Canada. Atmosphere 11, 1363. https://doi.org/10.3390/atmos11121363

Zhang, Y., Wang, X., Fraser, R., Olthof, I., Chen, W., Mclennan, D., Ponomarenko, S., Wu, W., 2013. Modelling and mapping climate change impacts on permafrost at high spatial resolution for an Arctic region with complex terrain. The Cryosphere 7, 1121–1137. https://doi.org/10.5194/tc-7-1121-2013