Rock Slope Temperature Evolution and Micrometer‐Scale Deformation at a Retreating Glacier Margin

In deglaciating environments, rock mass weakening and potential formation of rock slope instabilities is driven by long‐term and seasonal changes in thermal‐ and hydraulic‐ boundary conditions, combined with unloading due to ice melting. However, in‐situ observations are rare. In this study, we present new monitoring data from three highly instrumented boreholes, and numerical simulations to investigate rock slope temperature evolution and micrometer‐scale deformation during deglaciation. Our results show that the subsurface temperatures are adjusting to a new, warmer surface temperature following ice retreat. Heat conduction is identified as the dominant heat transfer process at sites with intact rock. Observed non‐conductive processes are related to groundwater exchange with cold subglacial water, snowmelt infiltration, or creek water infiltration. Our strain data shows that annual surface temperature cycles cause thermoelastic deformation that dominate the strain signals in the shallow thermally active layer at our stable rock slope locations. At deeper sensors, reversible strain signals correlating with pore pressure fluctuations dominate. Irreversible deformation, which we relate with progressive rock mass damage, occurs as short‐term (hours to weeks) strain events and as slower, continuous strain trends. The majority of the short‐term irreversible strain events coincides with precipitation events or pore pressure changes. Longer‐term trends in the strain time series and a minority of short‐term strain events cannot directly be related to any of the investigated drivers. We propose that the observed increased damage accumulation close to the glacier margin can significantly contribute to the long‐term formation of paraglacial rock slope instabilities during multiple glacial cycles.

boundary conditions (McColl & Draebing, 2019). The present work is focused on field investigations of the transient thermal regime induced by glacier retreat in adjacent valley slopes, and the related subsurface deformations.
As glaciers retreat, the stable thermal boundary conditions below glaciers are perturbed, which results in a "thermal shock" and long term cyclic thermal loading, that can induce long term rock mass damage (Grämiger et al., 2018). Until now, investigations on this topic were mainly based on numerical models and relied on assumptions unconstrained by in-situ data (e.g., Grämiger et al., 2020;Riva et al., 2017). The data presented herein is key to the understanding of heat transfer mechanisms in the subsurface, the timescales over which valley flanks adjust to the new thermal boundary conditions following deglaciation, as well as the magnitudes of thermomechanical stresses, which can lead to long term rock mass damage. Investigating these phenomena requires high resolution measurements of subsurface temperatures and deformation, which have never been recorded in alpine regions that feature temperate valley glaciers. Additionally, thermomechanical deformation signals can be overprinted for example, by hydromechanical deformation (cf., Grämiger et al., 2020), and disentangling the two requires an understanding of slope hydrology, as well as measurements of pore pressure and precipitation events. Detailed investigations of hydraulic processes and hydromechanical deformations in our test site are presented in a companion publication.
When considering the thermal regime of a rock slope, the glacier-rock slope interaction is controlled by both thermal and hydraulic conditions at the ice-rock interface. Bedrock in contact with temperate ice stays at the pressure melting point of the ice, which is around 0°C (Wegmann et al., 1998). When ice retreats, the previously covered bedrock surface is newly exposed to ambient temperature variations and solar radiation. These warmer temperatures then propagate to depth through a combination of conduction and advection.
A conduction dominated subsurface temperature regime, which is expected in intact crystalline bedrock, is controlled by: (a) the surface temperature, which is subject to variability at daily to millennial timescales, (b) an outward flow of heat from the earth's interior, which can be considered constant at the relevant timescales, (c) the thermal diffusivity of the medium, (d) latent heat effects from melting ice or freezing water in pores/cracks, and (e) the topography (Guéguen & Palciauskas, 1994;Pollack et al., 1998). In the case of homogenous subsurface thermal properties and flat topography, a steady-state temperature profile shows a constant rate of temperature increase with depth (dependent on the thermal conductivity and the heat flow from below). Surface temperature fluctuations then propagate as thermal waves at a pace governed by the thermal diffusivity, with amplitudes diminishing exponentially with depth (Pollack & Huang, 2000).
In reality, rock masses contain discontinuities, where fracture characteristics such as density, aperture, persistence, and connectivity control non-conductive heat transfer processes (e.g., advection of air or water), and heterogeneities (e.g., in the lateral thermal conductivity) potentially disturb the idealized heat conduction dominated subsurface temperature regime (Ge, 1998;Moore et al., 2011;Rybach & Pfister, 1994;Wegmann & Gudmundsson, 1999;Welch & Allen, 2014). Of particular importance to the present work are advective heat transfer processes driven by seasonal changes in rock slope hydrology in a glaciated valley, which influence the thermal regime. Transient pore pressures in such rock slopes are controlled by direct recharge to the fractured bedrock slope from snow melt in late spring and rainstorms in summer and fall, as well as by the englacial water pressure variations. Temperate glaciers carry large amounts of liquid water, and a thin water film with pressures dependent on the glacial hydrological system is assumed to be present at greater depth at the glacier-rock interface (Lappegard et al., 2006). Different types of temperature-coupled processes such as frost weathering, freeze-thaw cycles, and permafrost degradation contribute to progressive rock mass weakening in alpine regions (Deprez et al., 2020;Draebing & Krautblatter, 2019;Draebing et al., 2017;Hales & Roering, 2007). Additionally, thermomechanical stresses themselves can drive deformation and progressive rock mass damage in rock slopes (Gischig et al., 2011a(Gischig et al., , 2011bMarmoni et al., 2020;Mufundirwa et al., 2011). In the context of post-and late-glacial ice retreat, numerical investigations suggest that thermomechanical effects induced by long-term temperature variations related to glacial cycles and seasonal temperature variations induce significant thermomechanical stresses to the valley flanks which can promote progressive rock mass damage (Baroni et al., 2014;Grämiger et al., 2018). These studies have shown that, during a glacial cycle, thermomechanically driven deformation and rock mass damage can be substantially higher than purely mechanical loading and unloading effects from glacier retreat or advance, and that the degree of fracturing (and the criticality of the fracture system in terms of its Mohr-Coulomb strength) is a an important factor controlling the amount of induced thermomechanical damage (Grämiger et al., 2018).
The aforementioned studies (Gischig et al., 2011a(Gischig et al., , 2011bMarmoni et al., 2020;Mufundirwa et al., 2011) have shown that thermomechanical stress perturbations (or stress perturbation of any origin) can drive reversible (elastic) deformation and irreversible (plastic) deformation related to progressive damage. The deformation characteristics strongly depend on the magnitude of the stress perturbation, as well as the rock mass properties and spatially varying in-situ stress levels that control the criticality of fractures in the system (Gischig et al., 2011a(Gischig et al., , 2011bStock et al., 2012;Styles et al., 2011). A critically stressed fracture has a ratio between shear stress and normal effective stress which approaches the Mohr-Coulomb failure envelope. Stress perturbations that occur at levels exceeding critical stress thresholds can trigger instantaneous plastic (irreversible) deformation for example, in form of fracture slip (shear failure). If a stress perturbation occurs at a level below the critical stress threshold, it either mainly causes elastic (reversible) deformation or if it exceeds a so-called fatigue limit it can lead to slow subcritical fatigue crack growth and hence minor plastic deformation (Attewell & Farmer, 1973;Brain et al., 2014). Processes contributing to subcritical crack growth can be both of cyclic nature (cyclic fatigue) or static (stress corrosion) (Cerfontaine & Collin, 2018;Eppes & Keanini, 2017). This has been observed in laboratory experiments and is described by several authors for example, Attewell and Farmer (1973); Brown and Hudson (1973); Cerfontaine and Collin (2018). Thus, such low stresses can also increase the criticality of a fracture to a level exceeding the fracture toughness, which enables the fracture to fail at an apparently random, low magnitude stress perturbation, complicating correlations with environmental drivers (Eppes & Keanini, 2017).
In this study we use novel borehole monitoring data and numerical simulations to investigate the transient subsurface temperature regime and related micrometer-scale deformation in a rock slope adjacent to a temperate valley glacier undergoing rapid downwasting. We first analyze the borehole temperature and strain data and perform a statistical analysis of the strain data to highlight correlations between potential drivers and deformation events. We then use the borehole temperature data to estimate the depth-dependent apparent thermal diffusivity at our borehole locations, which is used to parameterize continuum 2D finite element models for heat conduction in the subsurface. Comparing these modeling results with temperature monitoring data provides insights into the transient subsurface temperature regime that follows glacier retreat, as well as the dominant mechanisms governing heat transfer at our study site (conduction vs. advection). Finally, we expand our finite element models to compute the expected thermoelastic deformation resulting from annual temperature cycles. These analyses provide crucial information regarding the relative importance of thermomechanical deformation on the total strain recorded at our instrumented rock slope, and hence its relevance for the overall observed deformation and progressive rock mass damage.

Geology and Geomorphology
Our study site is located on the southern valley slope adjacent to the retreating glacier tongue of the Great Aletsch Glacier (Canton Valais, Switzerland) (see Figure 1). The bedrock consists of fractured crystalline rocks (mainly gneisses and granites with thin schist layers) of the Aare-Massif (Steck, 2011). The rock slope is dipping NNW with an average slope angle of about 30° and the morphology is characterized by glacially smoothed bedrock surfaces with ridges and furrows. Furrows often follow weak schist layers or brittle-ductile shear zones normally oriented subparallel to the alpine foliation. This foliation strikes valley subparallel with a subvertical dip angle normally dipping into the slope. The rock mass is intersected by three persistent joint sets with a normal spacing around 1-3 m at the ground surface (for joint set orientations refer to Figure  4 in Hugentobler et al., 2020). 10.1029/2021JF006195 4 of 33 A thin and patchy recent till layer, surrounding several rock outcrops, covers the lower study area. Due to the recent glacial retreat, no noteworthy vegetation is present in the study area. The Moosfluh instability, a well-documented deep-seated gravitational slope deformation (DSGSD) with superimposed shallower rockslides (Glueer et al., 2019(Glueer et al., , 2020, is situated in the southwest of the study area. The glacial history of the Great Aletsch Glacier, the largest ice mass in the European Alps, is summarized in Holzhauser et al. (2005) and glacial extents in the vicinity of the study area have been investigated in detail by Grämiger et al. (2017). Additionally, glacial retreat from the last local maxima in 1860, the Little Ice Age (LIA), is well documented by GLAMOS (1881-2019). Since the LIA, the glacier tongue lost more than 3,000 m of length and over 300 m in ice thickness at the location of our study site. This leads to a mean ice height loss of about 2 m per year since LIA. However, annually acquired, high-resolution digital elevation models of the Great Aletsch Glacier area from the years 2012 until 2019 (source: Swiss Federal Office of Topography) show that the current melting rates are higher, with about 10 m ice height loss per year.  Hugentobler et al., 2020). The former glacial extents of years 1927-2014 were digitized from historical aerial photos by Glueer et al. (2020) and the Little Ice Age (LIA) and Egesen maximum extents are well constraint by lateral moraines in the map area. We extended the geological map from Steck (2011) in the area that had been glacier covered during their mapping campaign. The numbered solid black lines are the extents of the 2D profiles used for the numerical study.

Subsurface Monitoring System and Previous Work
In summer 2017, a subsurface monitoring system has been installed to investigate and quantify the effects of thermo-hydro-mechanical rock slope processes that occur along with downwasting and retreat of glacial ice (Hugentobler et al., 2020). The monitoring system is installed in three vertical, 44-50 m deep boreholes specifically drilled for research purposes and with a diameter of 10 cm (B2, B4, and B6, see Figure 1). The study site is considered as representative for many recently deglaciated rock slopes in inner Alpine valleys with its moderately to steep dipping slope, a crystalline lithology intersected by slope subparallel joints and alpine tectonic joint sets and faults. Similar geological preconditioning factors can be found in very many Alpine valleys, and presumably also in many other young mountain belts, such as the Himalaya and Southern Alps. The three borehole locations were selected to have varying lateral distance to the glacier margin, so that they monitor the rock slope conditions at different timings relative to deglaciation (Figure 1), and in slope sectors of different degree of stability and damage. Before the recent deglaciation, the borehole locations were ice covered for approximately 2000-3000 years (Holzhauser et al., 2005). The location of borehole B2 was for the first time uncovered in summer 2017, and the glacier returned in winter 2018 for about 4 months (cf., Hugentobler et al., 2020). According to high resolution aerial images (source: Swiss Federal Office of Topography), the borehole locations B4 and B6 were uncovered from glacial ice in summer 2014 and summer 2016, respectively. Current ice melting rates caused an ice height loss of about 30 m at our study site during the three years of monitoring.
The boreholes are fully grouted, and contain high resolution vertical FBG (fiber Bragg grating) strain sensors and SAA (shape accel array) in-place inclinometers for the measurement of horizontal deformation, a temperature sensor chain, and a pore pressure sensor installed in a 1-2 m long sand filter at the borehole end (see Figure 2). The setup and performance and early data of the monitoring system is described in detail in Hugentobler et al. (2020). Grouting strain sensors into boreholes is the state-of-the-art procedure to couple the sensors with the formation (Choi et al., 2021) and the grouting was conducted by a specialist company for geotechnical monitoring. The grout volume as well as strain and temperature in the borehole was monitored during the grouting process to assure a complete filling of the boreholes.
Air temperature is measured at our study site with a 30 min interval and total precipitation data used in this study is provided by MeteoSchweiz and recorded at the weather station "Bruchji" (canton Valais) located about 6 km away from our study site. Most intense precipitation events recorded at the Bruchji station during our monitoring period could be also confirmed at our study site by a specific seismic noise response measured at a seismic sensor installed at the study site by Manconi et al. (2018). We therefore think that it is reasonable to use the data of this weather station for our analysis.
The rock mass and the hydrogeological properties at our study site have been characterized with surface structural mapping, rock mass classifications (at outcrops surrounding the borehole sites and from borehole optical televiewer images), borehole geophysical and hydro-geophysical investigations, as well as pumping and injection tests in the boreholes. These analyses have identified differences between the three research borehole locations (Hugentobler et al., 2020). The main findings for the individual boreholes related to rock Borehole B2 was drilled directly at the ice margin at the end of a relatively flat plateau and only about 200 m away from the active Moosfluh slope instability. Optical televiewer images of the borehole wall revealed a high fracture density with open joints often filled with fines in the uppermost 20 m (Figure 2). Below this depth, fractures occur less frequently and with lower aperture. The rock mass at outcrops surrounding B2 was classified with a Geological Strength Index (GSI) between 55 and 65. According to the optical televiewer image, these values can be confirmed for the upper part of the borehole, but slightly increase with depth due to the decreasing fracture density (Figure 2). From a borehole water injection test, the transmissivity in the lower part of the borehole (around 15-50 m depth) was estimated around 10 −6 m 2 /s and in the uppermost, fractured part around 10 −3 m 2 /s. Borehole B4 is located about 50 m away from the glacier margin at a relatively evenly dipping slope. Borehole optical televiewer logs show that the rock mass is less fractured compared to B2 (Figure 2), with mainly closed or low aperture structures. A surficial layer of about 5 m with a dense fracturing is present at this location ( Figure 2). GSI values at outcrops around B4 reveal a high-quality rock mass with values between 75 and 85. Similar to B2, the rock mass quality appears to increase with depth with the observed decrease in fracture density. A pumping test in the borehole revealed a transmissivity of about 10 −6 m 2 /s between 9 and 44 m depth and higher assumed values in the uppermost section. Fluid electric conductivity logs show a major transmissive structure at the depth of 29.4 m and minor transmissive structures at 19.5, 25, and 39.8 m.
Borehole B6 was drilled at a cliff position with about 5 m lateral distance to the glacier margin. The area around the borehole is characterized by a disturbed rock mass bounded by a prominent head scarp above. Surface structural mapping indicates that this area is likely an inactive toppling rock slope instability, hosting important groundwater springs (see indicated "Inactive landslide body" in Figure 1). We estimated GSI-values at outcrops around B6 to range between 50 and 60. Because of unstable borehole conditions, a perforated PVC casing was installed in the borehole (cf., Hugentobler et al., 2020). Therefore, only fluid electric conductivity logging and a pumping test could be conducted. These tests resulted in an estimated borehole transmissivity of around 10 −5 m 2 /s with higher assumed values in the top layer and major transmissive structures at a depth of 35.5 and 43 m. A surface creek that originates from springs located a few tens of meters above the borehole location passes a depression next to the borehole and was observed to infiltrate through fractures into the borehole during drilling works. The lower GSI-values, the higher transmissivity, and the borehole stability issues encountered while drilling reveal the lower rock mass quality at this borehole location compared to the other two.

Analysis of Monitoring Data
In this paper, we present a detailed analysis of data from the subsurface temperature sensor chains, as well as from the high resolution vertical FBG strain sensors. Horizontal strain data from the SAA in-place inclinometer system is not discussed in detail due to space reasons but is provided in Section 5.3 ( Figure 12) to support our discussion on the overall observed strain characteristics and its implications for paraglacial rock slope evolution. Subsurface temperatures have been recorded since October 2017 with an hourly measurement interval, and a downhole sensor spacing of 1 m. Some data gaps of a few weeks in total exist in the time series that are related to battery issues. One longer data gap of four months (April-July 2018) in the B2 time series is related to a temporary glacier re-advance that overran the borehole location and disrupted the loggers from the sensor chains.
We use the subsurface temperature data to empirically estimate the in-situ apparent thermal diffusivity D within the annual thermally active layer. D is defined as where λ is the thermal conductivity, ρ the density of the rock mass, and C P the specific heat capacity of the rock. This was done by analytically solving the 1D heat conduction equation for measured periodic surface temperature cycling, and fitting our subsurface temperature data (cf., Koo & Song, 2008;Rajeev & Kodikara, 2016;Schneider et al., 2012).

of 33
In this analysis, annual surface temperature variations are approximated as periodic oscillations and expressed as Fourier series, which are then used as boundary conditions to solve the 1D (vertical to ground surface) heat conduction equation. From the resulting analytical solution with time-periodic surface temperature variations in a semi-infinite half-space, equations for phase delay (linear trend with depth) and amplitude decay (exponential decay with depth) can be derived (Equations 13 and 14 in Rajeev & Kodikara, 2016). These equations were then inverted to calculate the thermal diffusivity from the phase delay and amplitude decay of the annual temperature wave penetrating the subsurface at our borehole sites.
To estimate the phase delay, we applied a cross-correlation analysis between the temperature readings at the uppermost sensor and each other sensor in depth. The amplitude ratio between two time series recorded at different depths was calculated by first aligning the time series using the phase difference from the cross-correlation analysis and then performing a linear regression analysis between the two aligned temperature time series. The slope of the linear regression equals the amplitude ratio. We used these results to interpret which borehole locations are conduction dominated, as well as to parametrize our 2-D continuum heat transfer models.
The FBG sensors are installed as a serial connection of 10 individual sensors with base lengths of 4-5 m and cover the full length of each borehole (see Figure 2). The principles of operation of FBG strain sensors are described in detail in the literature (Morey et al., 1990;Zhu et al., 2017). The individual strain values were calculated by averaging over hourly dynamic measurements of approximately 30 s duration. In Hugentobler et al. (2020) we introduced the strain signals measured with the FBG-system and showed that the system is capable of reliably detecting strains with a precision below 1 µε (or 4 μm).
The individual FBG strain sensors show long term continuous changes in strain (positive strains equal extension and negative strains contraction of the monitored base length), distinct short-term strain events (within hours to days), and rapid strain steps (within one measurement cycle of 1 hr) ( Figure 6, Figures A1-A3). We compared the strain time series of the individual sensors in all three boreholes to variables of potential drivers for deformation (precipitation, pore pressures, and surface temperature). Additionally, we applied an event detection procedure to specifically investigate drivers for short-term and rapid strain variations (hereafter called strain events). An explanation of the different strain signals is provided in the beginning of Section 4.1.3.
To investigate potential drivers for strain events in the data and their distribution over the year, the timing of the most significant strain events in the time series were automatically retrieved using a MATLAB-based function. This function finds the most significant change points in a time series based on statistical properties (in our case mean and slope) and a predefined maximum number of change points (Killick et al., 2012;Lavielle, 2005). The function proceeds in the following steps: (a) Choose a division point in the time series and divide it into two sections; (b) Compute an empirical estimate of the desired statistical property (in our case mean and slope) for both sections; (c) Measure the deviation between the empirical estimate and each point in the section and add the deviations for all points; (d) Add the deviations section-to-section to find the total residual error; (e) Vary the location of the division point until a minimum total residual error is found. A maximum number of 50 changing points per strain time series were defined by trial and error for both the mean and slope criterion and their results were intersected because this yielded the best results. With this approach we did not detect every single strain event but sampled the most significant strain events present in the time series. All automatically detected strain events were then manually checked. A minority of detected events were corrected because they were artifacts related to data gaps, but the majority could be confirmed.
The timing of these most significant strain events was then compared to the potential drivers or causal factors for deformation such as pore pressure variations, precipitation events, and exceptional positive or negative surface temperature peaks or rapid temperature shifts (hereafter called extreme surface temperatures). We expect that the causal factors of strain events may be depth dependent because annual surface temperature cycles mainly affect the shallow subsurface and hydromechanical effects related to pore pressure fluctuations are expected to be more important at greater depth (i.e., below the phreatic groundwater surface of the slope). Therefore, we distinguish deep and shallow strain events. This allows us to analyze the different drivers or causal factors of short term and rapid strain events, which are potentially related to rock mass damage.

Numerical Simulation
We performed 2D, finite-element (FE) numerical modeling in order to understand long term (on the order of 200 years) changes of the slope thermal regime caused by glacial retreat, as well as yearly changes in the thermal regime, and resulting thermo-elastic strain, driven by seasonal temperature cycles. 2D numerical simulations were conducted using the commercial finite element software COMSOL Multiphysics v 5.5. The model geometry was defined using 2D cross-sections that intersect the three borehole sites, and laterally extend from the glacier-filled valley bottom to the ridge of Moosfluh/Hohfluh (see Figure 1). Surface topography above the current glacial ice level was defined using a high-resolution DEM from 2018 of the study area (source: Swiss Federal Office of Topography). Elevations of bedrock sections currently covered by glacial ice were interpolated using data from Rutishauser et al. (2016) in GlaThiDa_Consortium (2019).
To allow a comparison with monitored temperature and strain data recorded in our rock slope at high spatial resolution, we used a triangular FE mesh geometry with two domains (see Figure 3). Domain 1 covers the surrounding of the borehole locations and is composed of a fine mesh with element sizes between 0.03 and 18 m and increasing mesh size with normal distance to the ground surface. This resulted in individual mesh sizes of a few centimeters to maximum a few decimeters along the 50 m deep modeled borehole. Domain 2 contains a coarser mesh with mesh sizes between 0.5 and 95 m and an increase of mesh size with normal distance to ground surface and Domain 1.
Heat conduction is the only heat transfer process occurring in these continuum models and our models solve the heat conduction equation defined as ity D is mainly based on the empirical estimations described in the previous section, and the C P and ρ values used in the 2D model are based on literature values, as summarized in Table 1. We used a depth dependent thermal diffusivity in our models to partially account for the depth dependent structural heterogeneity present at our study site. Between 0 and 10 m depth we assigned the mean estimated in-situ value based on the inversion procedure detailed above. From 20 m downwards a constant literature value for intact rock thermal diffusivity was used (Table 1), and in between these two values the thermal diffusivity was linearly interpolated. The transient subsurface temperature field in the model was calculated by assigning boundary conditions (BC) of a constant geothermal heat flux at the bottom, a transient surface rock temperature at the top (described in more detail below), and no-flow lateral boundaries ( Figure 3). To achieve a realistic geothermal gradient of ∼23°C-km −1 in depth (cf. Rybach & Pfister, 1994), the geothermal heat flux at the bottom boundary of the model was set to 60 mW m −2 which agrees to values used in literature (Wegmann et al., 1998).
The numerical modeling was performed in four steps, with the different boundary conditions visualized in Figure 3 and the thermal boundary conditions detailed in Table 2: 1. The steady-state subsurface temperature field was calculated for the LIA ice extent (see Figure 9a) using 0°C at the top boundary below ice elevation (BC1) and above the ice elevation an elevation dependent mean annual ground temperature (MAGT) function (BC2). 2. Starting from the steady-state temperature field of modeling step 1, a transient model was run to simulate the thermal effect of the lateral ice retreat from the LIA (∼1,860) to the 2017 ice extent. Historical aerial photos show that the lateral ice retreat along our modeled profiles correlates linearly with the documented glacier length change of the Great Aletsch Glacier (GLAMOS, 1881-2019), which has been monitored in detail. Therefore, we used the rate of glacial length change to calculate the lateral ice extent in each time step of the model. With this method, the modeled timing of deglaciation at the borehole locations deviate by a maximum of 2 years compared to the actual year of deglaciation. A constant temperature of 0°C (BC1) was assigned to the part of the top boundary below this lateral ice extent, and an elevation dependent MAGT function (BC2) was assigned above. 3. Following the simulation of ice retreat from the LIA, seasonal temperature variations were modeled using the measured temperature recorded at the uppermost sensor of borehole B4 (BC3), superimposed on the long-term deglaciation trend since LIA. In this modeling step, the temporal resolution was increased  to 1 hr in order to investigate the effects of seasonal temperature cycles on the subsurface temperature field. Similar to modeling step 1 and 2, BC1 was assigned at the top boundary below glacial ice. To achieve the best possible fit of the initial conditions with measured temperatures in borehole B4 we varied the model parameters for modeling step 2 (LIA retreat) by using a constant intact rock thermal diffusivity of 1.4 * 10 −6 m 2 /s instead of the 2-layer thermal diffusivity. Then modeling step 2 was run for 158 years (1860 (LIA) to 2018). Following this, we applied a spin-up period of one annual cycle using again the 2-layer thermal diffusivity and BC3 top condition (Table 2) to allow the models' initial conditions to equilibrate with the boundary conditions. For the simulation of the seasonal temperature cycles we then used the 2-layer thermal diffusivity introduced above. 4. Following the thermal modeling steps 1 to 3, we computed the expected thermoelastic strains resulting from seasonal surface temperature cycles measured at our site using partially coupled thermomechanical numerical simulations. The simulations are partially coupled because temperature changes influence stresses, but stress changes do not influence thermal boundary conditions, thermal properties, or the temperature field. We used a simplified continuum modeling approach and homogenous and isotropic elastic rock mass properties, following Grämiger et al. (2017) (see Table 1). The thermomechanical modeling was conducted in two steps using the same model geometry as for the thermal models. First, we initialized the stresses in the rock slope with a steady-state simulation applying mechanical forces only. Far-field stresses that are of tectonic and exhumation-induced origin were modeled in a simplified way with a horizontal to vertical stress ratio of k = 1 (Kastrup et al., 2004) at the lateral boundaries of the model (Figure 3c). We applied a roller boundary (i.e., constrain the displacement in direction normal to the boundary) at the bottom and a fix point constraint (i.e., no displacement at this point) at the center of this boundary line to keep the model in place. In addition, we added gravity and modeled a constant glacier load at the 2017 ice extent as a hydrostatic stress boundary (cf. Grämiger et al., 2017). After initialization, we switched the lateral mechanical boundaries to roller boundaries (see Figure 3d) and calculated the thermal modeling step 3 including a thermomechanical coupling through the thermal expansion coefficient α. We used the same thermal initial and boundary conditions as described in modeling step 3. Figure 4 shows the subsurface temperature evolution of the glacier adjacent rock slope from 2018 to 2020 at the three borehole locations. The graphs in the three first columns of the figure (a, b, c, e, f, g, i, j, k) show the annual temperature variation in the form of monthly mean temperature profiles, and the graphs in the column to the right of the figure (d, h, l) show the evolution of the annual mean temperatures of the three monitored years. Since we started to monitor the temperature in fall 2017, our provided annual temperature cycles start in November and end in October of the following year (e.g., 2018 contains data from November 2017 until October 2018). Figure 4 shows that the depth of the thermally active layer in B2 and B4, is around 17 m, and annual temperature cycles smoothly diffuse into the subsurface. In contrast, in borehole B6 the annual thermally active layer is around 28 m deep, and the temperature profile is highly irregular. All three boreholes show weak positive and negative temperature anomalies along the profile.

Temperature Data
Borehole B2 is the closest of the three boreholes to the glacier margin, and shows the lowest mean surface temperature in 2018 and a clear warming trend in the two following years. This warming trend is likely BC3 GT above ice elevation T(t,z) = 8.975-0.005*z + T TopB4 (t)  related to the significant change of distance to the glacial ice at this borehole location (from 0 m in 2017 to about 40 m in 2020). The annual variability in this borehole is around of 7-8°C at the uppermost sensor. B4, which is the borehole located furthest from the glacier, shows a mean surface temperature of around 4°C with an annual variability of around 8-9°C. B6 has an intermediate distance to the glacier, and shows an annual mean surface temperature of around 3°C and a variability of around 9°C. All three boreholes display a general warming trend, which can be observed as a shift in positive direction of the mean annual temperature profiles between 2018 and 2020 (Figures 4d, 4h, and 4l). This warming trend can be overprinted by annual temperature variations in the uppermost borehole sections (e.g., boreholes B2 and B6).
The slopes of the mean annual temperature profiles in the three boreholes reveal that the warming front has propagated to different depths in each of the three boreholes. Starting with B2, Figure 4d shows a negative temperature gradient (i.e., cooling downwards) to a depth of about 17 m (in 2018 and 2019). Between 17 and 25 m, the temperature during these years were anomalously low, with measured values around 0°C in years 2018 and 2019, and showing a clear increase in 2020. Below 25 m the temperature increases with depth with a gradient of about 0.2°C/10 m (i.e., similar to the assumed geothermal gradient of 23°/km in the area [Rybach & Pfister, 1994]).
In B4, the temperature profile ( Figure 4h) shows a cooling downward gradient along the whole length of the borehole, although this gradient reduces with depth. This reduction indicates that the trend may switch to warming downward below the borehole depth, similar to the trend observed in B2 below 25 m.
Borehole B6 shows a general temperature decrease with depth below a depth of 10 m, but with a more irreg-  Figure 5 shows the results of the thermal diffusivity inversion based on periodic annual temperature cycles measured during 2019, which include separate estimates based on amplitude decay (Figure 5a), phase delay (Figure 5b), and their theoretical relationship (Figure 5c). These results show a regular amplitude damping and phase delay with depth in boreholes B2 and B4, while B6 shows strong deviations from the one-dimensional heat diffusion model. Figure 5c also shows that the annual temperature signal measured in B2 and B4 plots closely on the theoretical line, suggesting that conductive heat transfer dominates in these two boreholes. Temperature data measured in borehole B6 does not follow this theoretical relationship, hence, non-conductive processes clearly dominate the subsurface thermal regime in the upper 20 m of this borehole.

Thermal Diffusivity
The best fit values for the mean apparent thermal diffusivity are 2.1 * 10 −6 m 2 /s for B2 and 2.2 * 10 −6 m 2 /s for B4, respectively. Around B4, the slope shows a relatively constant dip angle, so the 1-D approach is likely more reliable at this location compared to B2, where 2-D effects are assumed to play a more important role. Therefore, the apparent thermal diffusivity estimate of B4 was used for further investigations. Figure 5d shows the inverted depth dependent apparent thermal diffusivity for B4, based on both the amplitude decay and phase delay. The estimates from both methods show a similar trend, with decreasing apparent thermal diffusivity values with depth. The highest values occur in the uppermost 4 m and below a more gentle decrease with depth is observed. The estimated magnitudes of thermal diffusivity are higher than literature values for conductive heat transfer in intact gneissic rocks would suggest (∼1.4 * 10 −6 m 2 /s) (Wegmann, 1998;Wegmann et al., 1998). We expect the estimated apparent thermal diffusivity values within the thermally active layer to approach a more constant value close to the intact rock thermal diffusivity at greater depth. Hence, the higher estimated value is assumed to be an 'apparent' value that also comprises effects related to the increased fracture density and water/air content in the surficial layer that allow minor non-conductive processes to occur. Such minor non-conductive effects should not be confused with stronger advective heat transfer processes, which are likely important in B6 and can result in significant disturbances of the diffusive profile (e.g., Ge, 1998).

FBG Vertical Strain
We recorded both short-and long-term strain signals with the FBG monitoring system. Two characteristic examples of strain time series are provided in Figure 6, and the complete FBG strain time series of all sensors in the three boreholes including variables of potential drivers for deformation are provided in the Appendix A ( Figures A1-A3).
These strain types, which can be identified in all boreholes are: (a) multi-annual trends in strain that mainly occur in negative direction (e.g., Figures 6a and 6b), (b) annual strain cycles (e.g., Figure 6a), (c) rapid and short term strain events typically showing magnitudes of up to around 20 µε, occur in positive and negative direction (e.g., red arrows in Figure 6), and can be detected with the algorithm explained in Section 3.1.
A comparison of Figures A1 and 3 (Appendix A) shows that the strain response at the three borehole locations varies, but also shows similarities of specific strain signals.
In all three boreholes multi-annual negative strain trends (contraction) are observed at several sensors. In B2 and B6 ( Figures A1 and 3), these negative strain trends are stronger in shallow sensors, while in B4 (Figure A2) they predominate at deeper sensors.
The uppermost sensors that monitor strain from around 0 to 5 m depth show an overall positive (extension) strain trend and some annual cyclic behavior with extension during wintertime and contraction during summer season.
Annual cyclic signals in the opposite direction (i.e., with contraction during the cold wintertime and extension during the warm summertime) are well visible in FBG-2 and 3 of B4 and are further investigated in Section 4.2.3. These signals show a slight phase shift with depth. Similar annual cyclic signals occur in shallow sensors of B2 and B6 but are less clear because they are superimposed by short-term or rapid strain events that occur at similar or greater magnitudes.
These larger strain events with magnitudes between 20 to a few 100 µε occur most frequently in borehole B2 ( Figure A1), and could either be significant deformation events, or an artifact of the FBG signal processing. Most of these steps could either be verified by simultaneous deformation detected with the independent in-place inclinometer system in the borehole (cf., Hugentobler et al., 2020), occurred at several FBG sensors simultaneously (b in Figure A2), or coincided with potential triggers (e.g., rainfall events or pore pressure changes; b in Figure A1). Hence, they are interpreted as significant deformation events. Only one particular large step in sensor FBG-4 of borehole B2 (small form letter a in Figure A1) is likely a signal processing artifact and was removed in the analyses that follow. The most significant large magnitude strain event occurred in B2 at several sensors simultaneously and coincided with an extreme rainstorm event (label b in Figure A1). Figure 7 presents a comparison of automatically detected strain events of the FBG-sensors with potential drivers or causal factors for deformation, including surface temperature, pressure head, and precipitation. Results presented in Figure 7 are summarized in the pie plots of Figure 8, which shows the percentage of the strain events coinciding with the potential drivers or causal factors at the specific boreholes. The data shows that, depending on the borehole, between 60% and 75% of the strain events coincide with precipitation events (Figures 7b, 7f, 7j, and 8). Only about half of these precipitation events caused a clear pressure head increase (Figures 7e, 7i, and 8). On the other hand, 5%-9% of the strain events temporally correlate with pressure head increases not related to direct rainfall infiltration but to snowmelt infiltration. Very few strain events (<2%) coincide with extreme surface temperatures (Figures 7a, 7d, 7h, and 8). These results emphasize the importance of hydromechanical and potentially moisture-related processes for short-term and rapid deformation. However, the interactions are complex, as rainfall also modifies the mechanical and hydraulic conditions of the adjacent glacier. Analyzing these complex interactions is beyond the scope of the present paper and are explored in a future contribution. Further, the conducted empirical attribution of strain events to extreme surface temperatures is strongly simplified, and a more quantitative comparison of monitored strain to modeled thermo-elastic strain is provided in Section 4.2.3.  Between 20% to 30% of the strain events do not coincide with any of the three potential causal factors. To investigate if these events could be related to earthquakes, which are also frequently discussed as destabilizing factors during deglaciation (e.g., McColl & Draebing, 2019), we compared the timing of the strain events with measured earthquakes in the observation period (data provided by the Swiss Seismological Service, SED). Earthquakes that occurred at epicentral distances less than 50 km from the study site and having magnitudes greater than M1 and an earthquake depth of less than 15 km were used for the comparison. However, we did not find any correlation between the around 70 recorded earthquakes during the monitoring period (maximum magnitude M3) and the measured strain events in our boreholes. A characterization of all earthquake events used for the comparison (i.e., magnitude, epicentral distance to study site, and earthquake depth) can be found in the Appendix B ( Figure B1).
The cumulative numbers of detected strain events (Figures 7c, 7g, and 7k) allow us to better visualize the variation of the number of events occurring over the year. Data in these plots is provided as mean number of strain events per sensor interval for shallow sensors (within the annual thermally active layer, i.e., the uppermost ∼18 m) and for deep sensors (∼18 m down to the borehole end). Figure 7 shows that fewer strain events occur during wintertime compared to the spring and summer season and that these strain events occur more frequently in deeper sensors. This difference between shallow and deeper sensors is more strongly pronounced in B4 and B6 compared to B2. In B2, some large magnitude strain events are measured in the shallow layer (see Figure A1, Appendix A). Figure 9 shows simulated subsurface temperatures during LIA deglaciation compared to our measured mean subsurface temperatures at the borehole locations. Figure 9a shows the temperature field of the steady-state solution during the LIA ice extent, which, as mentioned previously, is used as an initial condition for the following transient ice retreat simulation. The modeled temperature field in 2019 is shown in Figure 7. Number of detected strain events during the monitoring period in the three boreholes compared to potential drivers or causal factors for deformation. Panels (a), (d), (h) Daily sum of events versus air temperature at the study site. Panels (e and i) Daily sum of events versus pressure head measured at the borehole end (no functioning pore pressure sensor installed in borehole B2). Panels (b, f and j) Daily sum of events versus cumulative total precipitation per 24 hr measured around 5 km away from the borehole location at the weather station "Bruchji", canton VS (provided by MeteoSchweiz). Panels (c, g, and k) Mean cumulative sum of events per sensor provided for "shallow sensors" (i.e., down to 18 m depth) and "deep sensors" (i.e., from 18 m down to the borehole end).

Figure 8.
Pie plots for the three boreholes (B2, 4, 6) that summarize the percentage of the automatically detected strain events presented in Figure 7 that coincide with potential drivers or causal factors such as precipitation, pore pressure changes, or extreme surface temperatures. Because no functioning pore pressure sensor is installed in B2, a comparison with this potential driver cannot be provided for this borehole. Modeling temperature evolution and deformation during deglaciation. Figure 9b. Figure 9c illustrates the change of the subsurface temperature profile at the specific indicated mid-slope location during glacial retreat.
As can be seen on Figure 9c, the subsurface temperature profile before deglaciation shows 0°C at the ground surface (BC1 below ice) with a nearly linear temperature increase with depth. As soon as the ice retreats below position c, a new mean annual ground temperature is applied to the model (BC2 in Table 2) which results in warmer surface temperatures that penetrate downward with time. The adaption of the subsurface temperature to a new, warmer MAGT causes a transient inverse gradient in the upper subsurface temperature profile (i.e., cooling downwards), which transitions to a positive geothermal gradient at depth. Hence, an inverse thermal profile can be an indication for a transient subsurface temperature adapting to a new, warmer surface temperature (e.g., Pollack et al., 1998). The dashed line in Figure 9c illustrates a newly approached steady-state temperature profile assuming the new MAGT of around 2.6°C to be constant over a very long time at the elevation of the profile. The time required to reach this new steady-state temperature profile is dependent on the thermal diffusivity of the medium.
In order to allow a meaningful comparison of modeled temperatures with temperatures measured in the boreholes, the uncertainty related to the timing of deglaciation at the specific borehole location resulting from the lateral ice retreat function (max. 2 years, cf., Section 3.2) was corrected to get the true year of deglaciation in the models by running the model for maximum two more years. Consistent with the results shown on Figure 9c, all the measured temperature profiles show an inverse temperature gradient at least in the upper part of the borehole, which reflects the preserved cold temperature from the glacial occupation at depth and a more recent warming trend. Figure 9d shows that the measured temperature profile of B2 is not exactly reproduced with our conduction model and the applied boundary conditions. However, the modeled shape of the profile and depth of the warming front is similar to that measured at this borehole location (Figure 9d). The measured temperature profile is shifted by about 0.4°C to colder temperatures compared to the modeled values and has temperatures of around 0°C at a depth of 20-25 m.
Annual mean temperatures measured in borehole B4 can be well reproduced by the conduction model ( Figure 9e). Nevertheless, the model shows a slightly slower warming at a depth below the annual thermally active layer and the borehole end. The measured temperature profiles within the thermally active layer (i.e., the uppermost ∼20 m) deviate from the model because temperatures at these depths are subject to annual temperature variability that can differ from the used MAGT in the model.
Conversely, the measured temperature profile of B6 cannot be reproduced with our thermal conduction models. The measured temperature profile in this borehole shows a very different shape and quicker warming of the subsurface compared to the simulated ground warming after ice retreat by only heat conduction (Figure 9f). The warming front, evident by the turning point to a positive geothermal gradient, is below 50 m in the measured temperature profile and significantly shallower in the results of the conduction model.

Seasonal Temperature Variations
As described above, the results from the annual mean subsurface temperature simulation show good agreement with the monitoring data in borehole B4 (Section 4.2.1). Therefore, we investigated recent seasonal temperature variations (modeling step 3, Section 3.2) at this location, and compared results to monitoring data. Figure 10 shows a comparison of conduction-modeled and measured temperature profiles at different months during the year 2018. The results for other years are similar. The annual surface temperature variability, both in the model and the monitoring data, affects a depth range of approximately 20 m. Measurements and model show a good accordance.
To investigate the mechanisms of seasonal temperature variations in detail, the correlation coefficient between the modeled and measured temperature profile is calculated for each hourly time step and the R-squared value is plotted against time (see Figure 10g). The small form letters in Figure 10g indicate the timing of the different panels (Figures 10a-10f). The correlation coefficient time series shows a cyclic behavior with generally very high values (above 0.95) from midsummer to winter. Lower correlation values are observed during spring and temporally correlate with strong pressure head rise in borehole B4, caused by snowmelt infiltration (see Figure 10g). However, this negative peak in the R-squared value shows a delay compared to the pressure peak. Daily acquired photos from a time lapse camera installed at the opposite valley flank with view to the borehole locations show that this negative R-squared peak coincides with the time when the rock surface surrounding the borehole becomes snow free.

Thermoelastic Strain
Finally, we investigate the magnitude of thermoelastic strains in the subsurface resulting from seasonal temperature cycles (Figure 11). This is based on a comparison between monitored and modeled thermoelastic 10.1029/2021JF006195 19 of 33 strains in borehole B4 which is dominated by conductive heat transport. In addition, the pressure head measured at the borehole end is provided in Figure 11.
The numerical modeling shows that thermoelastic deformation occurs in the thermally active layer, and decreases with depth. The results show that the strain signal of the uppermost sensor cannot be reproduced by the simplified thermoelastic continuum model (Figures 11a). However, strain readings of FBG-sensors number 2 and 3 show a clear annual cyclicity with similar amplitude and phase as the modeled strain (Figures 11b and 11c). In the strain signal of sensor 4 (Figures 11d), a weak cyclic signal with an annual period and an amplitude similar to the modeled one might be superimposed to another signal strongly correlating with the local pressure head. In the deeper FBG-sensors, the modeled thermoelastic strain is an order of magnitude (or more) lower than the actual strain measured in the borehole. Furthermore, the strain variations in these depths often coincide with pressure head fluctuations.

Controls on Subsurface Temperature at a Retreating Glacier Margin
Heat transfer by conduction through the rock mass is an important process explaining many observations of our shallow paraglacial temperature field. Measured temperatures in borehole B4 can be well reproduced with the 2-D heat conduction models and the results for B2 suggest that conduction plays an important role as well. Conversely, the temperature profile measured in borehole B6 strongly deviates from an expected conduction dominated regime and hence cannot be explained with this process. Additionally, all boreholes display a warming trend characterized by an inverse temperature profile (cooling downwards) in the shallow subsurface, which is likely caused by glacier retreat and subsequent exposure of the slope to annual temperature cycles. The depth where the transition to a regular temperature profile (i.e., warming downwards) occurs is observed to increase based on the length of time a borehole has been ice free.
Non-conductive heat transport processes at our monitoring boreholes are potentially driven by groundwater interaction with the adjacent glacier and yearly recharge cycles on the ice free slope. Typical patterns of such advective temperature disturbances have been discussed by Ge (1998) who distinguishes between (a) distributed vertical groundwater flow in semiconfined porous aquifers, (b) horizontal flow in confined aquifers, and (c) localized groundwater flow in preferential groundwater pathways such as fractures in gneiss and granite, all causing substantially different signatures (Bredehoeft & Papaopulos, 1965;Ge, 1998). Such localized groundwater flow through fractures is assumed to be the cause for the observed positive and negative temperatures anomalies along the profiles (cf., Hugentobler et al., 2020). In order to develop localized thermal anomalies, the water flowing in transmissive fractures has to be in thermal disequilibrium with the surrounding rock and must be forced to flow by a natural head gradient, for example, resulting from a stationary or transient groundwater flow field in an inclined slope.
In borehole B2, Figure 9d shows a temperature offset of the modeled profile from the measured data of about 0.4°C between depths of 25 and 50 m. This deviation can be explained by cooling the subsurface with cold water directly connected to glacial ice, which could have caused temperatures of around 0°C down to a borehole depth of 25 m during the time of ice occupation. We suggest that the high fracture density in this shallow layer (see Figure 2) containing many open joints with confirmed high transmissivity from injection tests (cf. Section 2.2), as well as the relatively low normal distance (∼18 m at borehole meter 25) to the rock surface, enabled efficient advective heat transfer during time of ice occupation and during the years when the glacial ice was still close (i.e., 2018 and 2019). The measured temperature increase at around 25 m that is observed between 2019 and 2020 ( Figure 4d) can be related to a decoupling from the advective glacial hydraulic system due to the ongoing ice retreat and the penetration of the conductive warming front from the surface. For the shallow annual thermally active layer in B2, the modeled temperatures show stronger deviations from the measurements compared to B4. We explain this difference with the evident higher variability in MAGT at the B2 location probably related to its proximity to the ice margin and the simplified assumption of a constant MAGT in the model.
Relatively minor non-conductive heat transport processes are likely occurring at the B4 borehole location. This is evident when comparing the seasonal variability of the temperature profile measured in the borehole with the 1-D analytical solution and the 2-D numerical modeling results (Figures 5 and 10). The correlation values reveal an annual repeating decrease during the time of snowmelt in spring, indicating a secondary heat transport process superimposed on the purely conductive transport model (Figure 10g). Lower correlations are observed in the upper 20-30 m of borehole B4. The peak of the reduced correlation values occurs a few weeks after the pore pressure peak (recorded at 43.75 m depth), at a time when the borehole location becomes snow-free. We attribute the observed low correlation to advective heat transfer processes related to water infiltration happening directly at the borehole location.
Borehole B6 temperature measurements cannot be explained using a conduction model, the thermally active layer is 10 m deeper as compared to B2 and B4, and the measurements display larger temperature variations. This borehole is located within an inactive (presumably dormant) landslide body with substantially more disturbed and damaged rock mass compared to the other two borehole locations (cf., Section 2). Additionally, a creek fed by nearby regional springs passes through a morphological depression next to the borehole. During drilling in fall 2017, bubbles from the pressurized air used for drilling extruded all around this depression, revealing a high transmissivity in the shallow subsurface. It is assumed that the availability of creek water to infiltrate into the subsurface at this location and higher transmissivity due to the disturbed rock mass enables strong advective heat transfer processes to occur in the shallow subsurface, as documented in Anderson (2005).

Strains in the Thermally Active Layer
The analyses of Sections 4.1.3 and 4.2.3 have revealed differences in the drivers of the FBG strain signals as a function of depth, as well as between the borehole locations. All boreholes feature an annual cyclic strain signal in shallow sensors that shows consistence with ground temperature and a slight phase shift with depth. Given the correlation to ground temperature, as well as the results of the numerical modeling (Figure 11), it is likely that this component corresponds to thermoelastic deformation.
This can be clearly demonstrated in our monitoring data of borehole B4 with the most intact rock mass conditions. The FBG data from sensor 2 and 3 (Figure 11), which correspond to depth intervals of 4.6-8.9 m and 8.9-13.2 m, show that here the cyclic, thermoelastic deformation signal dominates over other signals. Deeper in the borehole (i.e., from ∼13.2 m downwards) other signals, that show minor or no correlation with the modeled thermoelastic deformation, dominate the strain time series. In borehole B2 and B6, with less intact rock mass conditions, the uppermost few sensors also show signals with an annual cyclicity that is attributed to thermoelastic effects. However, this continuous cyclic deformation is interrupted by abrupt strain events of similar magnitude or superimposed by sometimes even stronger deformation trends (e.g., label a in Figure A3). Additionally, most shallow sensors in all three boreholes show a long term negative strain trend.
An exception is the strain signal measured in the uppermost sensor of all boreholes, which shows an inverse correlation with the model data and hence must be related to a different driving mechanism. Expansion in wintertime and contraction in summertime could be caused by the effect of hygroscopic expansion or moisture induced strain (e.g., Gor et al., 2017) and related to potentially higher relative humidity in the shallow rock layer below the snowpack compared to dry summer conditions. Another possibility to generate extensional strain in wintertime and contraction in summer are frost weathering processes such as volumetric expansion or ice segregation (Draebing et al., 2017). For borehole B4 hygroscopic expansion is the more plausible process, as no negative subsurface temperatures were measured in wintertime in this borehole (Figure 4). In B2 and B6, the temperatures in the shallow 0-3 m were around 0 °C during some winter seasons (Figure 4), and therefore freeze thaw cycles could also be a reasonable explanation for observed cyclic strains in the uppermost sensor.

Strains Below the Thermally Active Layer
Some multi-annual strain time-series measured in deeper sensors clearly correlate with pressure head records in the corresponding borehole (e.g., sensor 4 and 5 in Figure 11) and hence are probably strongly driven by poroelastic hydromechanical effects. Other sensors show rapid or short-term strain events often temporally coinciding with pressure head changes or precipitation events but overall show a low correlation to the pressure head time-history (e.g., FBG-9, 10 in B4, Figure A2). We explain this diversity of reactions with the complex nature of hydromechanical drivers and coupled deformations in our fractured rock slope. The response of an individual FBG monitoring interval to a pressure perturbation strongly depends on the local stress conditions, and the number, orientation and type of fractures intersected by the interval (Krietsch et al., 2020).
In a paraglacial setting the impact of the adjacent glacier on deformations recorded in the slope can be manyfold and are presented in a future contribution. Effective glacial ice load changes, due to variations in the subglacial water pressures related to daily meltwater cycles and rainfall infiltration into the glacier during summertime, might be strong enough to cause measurable strain, pore pressure cycles and additional hydromechanically coupled effects. Also, the higher relative number of strain events detected in deeper, pore pressure controlled sensors compared to the sensors within the thermally active layer (evident in the cumulative sum plots of Figure 7), can be related to hydromechanical effects. During summer and fall a higher activity in the hydraulic systems of the slope and glacier must be expected, that is, groundwater recharge by snowmelt and rainfall infiltration induce recharge variations and higher hydraulic gradients in the rock slope, and daily ice melt in summer strongly increases the dynamics of the subglacial hydrological system (Harper et al., 2005). We interpret these effects to be the reason for the observed increased number of strain events occurring during the summer season.

Irreversible Rock Mass Damage
As discussed in the previous section, strain signals during our 2.5-year monitoring period are partly driven by thermomechanical and hydromechanical processes. The majority of our detected short-term and rapid strain events coincide with precipitation events, of which about 50% were strong enough to cause a clear reaction in the pore pressure signal at 40-50 m depth (Figure 8).
To identify rock mass damage, reversible (elastic) strains and irreversible (plastic) strains have to be distinguished (cf., Hugentobler et al., 2020), and the drivers for the specific strain signals need to be investigated. A strain signal is interpreted as reversible if it returns to its initial value after a given stress perturbation.
In contrast, irreversible strain is the offset that remains after a stress perturbation. However, for longerterm drivers of deformation, such as unloading due to glacial retreat, which is probably a continuously ongoing process at the time scale of our monitoring period, this distinction is not possible. In that case we use strain characteristics (e.g., orientation, spatial differences in strain signals) and additional data such as discontinuity occurrence and orientation to discuss if the signal may be reversible or irreversible. Damage mechanisms in fractured rock masses comprise of tensile or shear fracture propagation, breakage of intact rock bridges, degradation of asperities, smoothing of discontinuity surfaces due to shearing and wear, and formation of new fractures (Pellet & Selvadurai, 2017;Preisig et al., 2016). Measurable irreversible strain signals related to these damage mechanisms are either caused by shear or normal deformations along existing fractures which cross a monitoring interval, or the propagation of new fractures through a monitoring interval. The available data can not clearly distinguish the two processes, but considering the number and size of the mapped fractures and the generally high intact rock strength, we assume that deformation along existing fractures (e.g., slip) is the dominant process observed.
Hence, most rapid events (e.g., label b in Figure A2) are interpreted as irreversible fracture slip driven by pore pressure fluctuations or presumably also glacial load changes from varying englacial water levels. Short-term strain events often show similarities to pressure head signals (e.g., label a in Figure A2) with a strongly reversible strain component and a minor irreversible component. We relate the reversible portion of these strain events to poroelastic reactions and the often present minor irreversible component to plastic deformation and progressive rock mass damage. Therefore, many of the automatically detected rapid and short-term strain events contain at least a minor irreversible strain component, which means that the irreversible deformation is also strongly controlled by hydromechanical effects. This observation agrees with the numerical modeling results of Grämiger et al. (2020), who concluded that hydromechanical effects are more effective in generating long term damage in paraglacial rock slopes compared to thermomechanical effects.
The strain time series also contain 20%-30% short-term strain events that cannot directly be related to any of the investigated drivers or causal factors. We suggest that these signals may be related to fatigue processes occurring in the studied rock slope. As summarized in the introduction, cyclic stresses can promote rock fatigue and drive fracture propagation. In our studied slopes, cyclic loading is mainly driven by annual temperature variations, daily and annual pore pressure fluctuations and presumably also daily fluctuations of the mechanical load placed on the slope by the Aletsch glacier due to changing englacial water levels. In the shallow subsurface, cyclic thermal loading from annual temperature variations dominates the reversible signal causing largest strain amplitudes of ∼60 µε within the uppermost 4 m and strain amplitudes decreasing with depth to ∼5 µε at around 20 m depth. Below around 12 m depth, reversible strain cycles show amplitudes of up to ∼20 µε, are driven by pressure head fluctuations due to snowmelt infiltration or heavy rainfalls and normally occur at frequencies greater than 1/yr. The fact that these strain events often comprise a minor irreversible strain component, shows that many fractures may be close to a critical stress state. This could imply that fractures in our slope normally deform or grow slowly under subcritical conditions supported by cyclic fatigue, but that stronger loading events trigger slip along critically stressed fractures.
The overall negative strain observed in many sensors of all boreholes is composed of (a) a slow, continuous long-term strain signal that can be observed in several sensors ( Figures A1-A3, Appendix A) and (b) the accumulation of deformation from the short-term strain events discussed above. Because of the limited time of monitoring, it is not clear if these longer-term strain trends reflect reversible (elastic) or irreversible (plastic) deformation. Elastic deformation in such deglaciating rock slopes could be caused by rock mass extension due to unloading from ice downwasting or thermal expansion related to the observed ground warming after ice retreat. However, these perturbations would cause an extension of the rock mass, which is the opposite of what we observe in our monitoring data. On the other hand, compressive strain can originate either from slip along toppling fractures steeply dipping into the slope or from slip along sliding fractures dipping in downslope direction (cf., Figure 15 in Hugentobler et al., 2020). Fractures of both orientations frequently occur at our study site (Hugentobler et al., 2020). We therefore interpret these negative strain trends to be irreversible and related to progressive rock mass damage in our deglaciating rock slope.
We find that the magnitude of measured strain responses depends on the local rock mass quality at the given borehole location. This is observed in our research boreholes, as strong pore pressure variations (e.g., due to the heavy rainstorm in early October 2020) cause minor irreversible responses in the more intact rock mass of B4 ( Figure A2), but very clear irreversible deformation in the more strongly disturbed rock mass around B2 (label b in Figure A1). This observation is consistent with the numerical investigations of Gischig et al. (2011a) and Grämiger et al. (2020), who show that the efficiency of thermomechanical and hydromechanical loading cycles to drive progressive rock mass damage in fractured crystalline rocks strongly depend on the magnitude of the drivers, as well as the initial damage state and criticality of the fracture system.

Implications for Paraglacial Rock Slope Evolution
If our interpretations of the underlaying mechanisms hold, near-surface hygroscopic expansion can induce reversible annual strains, which are even larger than the thermo-elastic and poroelastic slope deformations occurring at greater depth. The magnitude of annual vertical strain (50-100 µε) induced by variations in moisture content is at the lower end of what has been derived from lab tests on intact and fractured granite with humidity increase from 20% to 80% (Li et al., 2021). Annual thermoelastic strain magnitudes of about 20-50 µε at 5-10 m depth are consistent with continuum modeling results and an in-situ derived thermal diffusivity ( Figure 11). Annual poroelastic strain cycles are less clearly defined in our borehole pressure records, mainly due to the complex interactions of englacial and slope hydrologic systems.
Modeled and measured longer-term subsurface temperature changes related to the ground warming effect after ice retreat could not be related to monitored strain at our depth of investigation. We explain this with the low expected strain magnitudes resulting from the minor temperature changes of below 1°C in three years, but also to the occurrence of other drivers such as pore pressure or presumably glacial load changes that cause stronger superimposed strain reactions dominating the signal below the thermally active layer. Hence, our findings agree with previous work that showed that thermomechanical effects contribute to progressive rock mass damage in deglaciating environments (Baroni et al., 2014;Grämiger et al., 2018), but that the longer-term effects are minor compared to the stronger hydromechanical and mechanical drivers that occur at greater frequencies during rock slope deglaciation.
As visible in our pore pressure data (borehole B4 vs. B6, Figures A2 and A3), local hydraulic properties of the rock mass control the magnitude of pore pressure fluctuations and therefore the relative importance of hydromechanical versus thermomechanical driven strain at a specific location within the slope. Our study only includes observations close to a retreating glacier margin. In fractured rock slopes the depth of the groundwater water table below surface generally increases toward the crest. In the upper left (southern) slope of the Aletsch Valley, the groundwater table elevation is expected to occur at several 100 m depth (Grämiger et al., 2020). Therefore, the relative importance of thermomechanical loading effects might increase with greater distance to the ice margin.
Our monitoring data shows that numerous rapid and presumably also slow slip events along steeply and moderately dipping fractures ( Figure 2) are superimposed on long-term seasonal strains and multi-annual strain trends. Cumulative slip reaches similar or higher displacement magnitudes than the elastic deformations reported above ( Figures A1-A3). Most of these slip events are irreversible and small (<20 µε) and the total displacement can locally be controlled by a few rare events, such as the 80 mm rainstorm of October 2020 (label b in Figure A1). All these events contribute to very slow irreversible slope movements, even in slopes considered stable. Cumulative displacement magnitudes of these movements as derived from the (horizontal) inclinometer readings range from 8 to 17 mm in two years (see Figure 12) and 2-10 mm for vertical irreversible deformation along our 50 m deep boreholes. Such creep rates are high compared to the irreversible horizontal slope displacements modeled by Grämiger et al. (2020) for a full glacial cycle. In their modeling results, highest displacements occur close to the slope toe (at the location of our monitoring boreholes), but only reach 25 cm after 500 annual cycles. One reason for the lower displacement rates in their numerical models could be related to the simplified assumptions regarding hydraulic and mechanical boundary conditions (especially at sub-annual timescales) and the neglection of subcritical fracture propagation mechanisms. On the other hand, our observed irreversible strain rates cannot occur permanently over thousands of years at the whole slope because this would cause unrealistic high overall slope deformations. This could imply that-as previously suggested-strongest cyclic loading, fatigue and damage occurs directly at the glacier margin and decays with lateral upslope distance. This might also explain the surprisingly high fracture apertures and borehole transmissivities recorded at the borehole locations, in a slope generally showing strong rock with high rock mass quality.
Following this hypothesis, long-term rock slope damage might be driven by accumulated damage occurring at the glacier margin during multiple glacial cycles. In the Great Aletsch Glacier valley the Holocene minimum is located very close (downstream) of our borehole locations (Figure 1). This implies that many Holocene glacier advance and retreat cycles occurred above our research borehole locations, and especially in the valley below the Holocene minimum (cf., Grämiger et al., 2017). This matches very well with the mapped landslide distribution in the Aletsch valley, showing a clear landslide density hotspot below the Holocene minimum. The initiation times of few dated landslides in this hotspot area fall into post-Egesen and pre-LIA time (Glueer et al., 2020;Grämiger et al., 2017), during which multiple glacier retreat and advance cycles have affected these rock slopes. This is also consistent with our proposed mechanism of paraglacial rock slope damage.
In addition to the spatial focusing of damage around the moving glacier margin, our observed high damage rates might also reveal a temporal focus of increased damage during times of strong glacial retreat. This is consistent with observed rock slope failures adjacent to retreating glaciers (Clayton et al., 2017;Evans & Clague, 1994;Fey et al., 2017;Kos et al., 2016;McColl & Davies, 2013;Oppikofer et al., 2008;Strozzi et al., 2010) and increased frequencies of rock slope failures directly after ice retreat revealed by a study of rock slope failure deposits in a Norwegian fiord (Böhme et al., 2015). The many rock slopes where a delayed response of thousands of years to deglaciation had been identified (e.g., Ballantyne, Sandeman, et al., 2014;Ivy-Ochs et al., 2009;McColl, 2012;Prager et al., 2008) could still have experienced increased damage rates during deglaciation. However, because of overall high rock mass quality, these increased damage rates were not strong enough to induce immediate failure.

Summary and Conclusion
We use three full years of hourly measured borehole temperature and high resolution FBG strain monitoring data to investigate the transient subsurface temperature regime and micrometer-scale deformation in a rock slope adjacent to a retreating, temperate valley glacier. The presented data give unique new insights into the thermomechanical and hydromechanical drivers of rock slope deformation adjacent to retreating glaciers over the monitored time scales. Detailed analysis of the subsurface deformation signals allow to identify various potential drivers for reversible (elastic) and irreversible deformation. This knowledge is important for understanding the main processes contributing to short-and long-term progressive rock mass damage that potentially lead to the formation of paraglacial rock slope instabilities. The most important findings of the present work can be summarized as follows: 1. The thermal investigation reveals that the subsurface temperatures at our monitored slope are currently in a transient state adjusting to new surface temperature conditions after glacial ice retreat. We show that heat conduction is a dominant process and explains most of the observed transient subsurface temperatures at sites with intact rock, but locally strong deviations from a conduction dominated regime can occur. These deviations can be explained by non-conductive processes (e.g., heat transfer by water) that become important when the rock mass is fractured (i.e., has an increased transmissivity), and are efficient when water and a hydraulic head gradient is available. The observed non-conductive temperature deviations in our study area are related to groundwater exchange with cold subglacial water, snowmelt infiltration, or creek water infiltration. 2. Our strain monitoring data analysis shows that annual surface temperature changes, precipitation events and pore pressure fluctuations are important drivers for deformation at our study site. However, longerterm trends in the strain time series and a minority of short-term strain events cannot directly be related to any of the investigated drivers or causal factors. We find that thermo-elastic effects due to annual temperature cycles dominate the strain signal in the uppermost 12 m of our stable fractured crystalline rock slope. Below this depth the hydromechanical effects dominate the deformation. Furthermore, we show that both reversible and irreversible short-term deformation more frequently occur during summertime when the hydrological system in the rock slope and in the temperate glacier are more dynamic. 3. Irreversible deformation in our time series, which we relate with progressive rock mass damage, occurs as short-term or rapid strain events within hours or days and as slower, continuous strain trends (mainly in vertical contraction direction). We show that the majority of short-term strain events coincide with precipitation events or pore pressure changes. Other potential causal factors such as extreme surface temperatures or earthquakes occurring in the area do not show temporal correlation with strain events. We propose that cyclic fatigue from thermomechanical and hydromechanical loading cycles, can support long-term subcritical crack grows and micrometer-scale slope deformations even in presumably stable slopes. 4. We show that the magnitude of irreversible deformation depends on the magnitude of the driver but also on the spatially variable pre-existing damage at our three borehole locations. At our most stable borehole location, single events with irreversible strain magnitudes are normally on the order of a few με, whereas at the more damaged borehole location irreversible strain magnitudes up to tens or few hundreds of με have been recorded. 5. We propose that strongest cyclic loading, fatigue, and damage processes occur close to the glacier margin and decay with lateral upslope distance. Following this hypothesis, long-term damage accumulation close to the moving glacier margin during multiple glacial cycles could significantly contribute to the formation of paraglacial rock slope instabilities.

Appendix A: FBG Strain Monitoring Data
Figures A1-A3 contain the temperature corrected hourly mean FBG-strain measurements for each sensor in the three vertical research boreholes. The setup and detailed description of the monitoring system is provided in Hugentobler et al. (2020). All FBG sensors of the three boreholes contain a data gap of about one month in January 2019 which is related to a winter storm that broke the solar panel used for the power supply of system. A second data gap of around two month duration is existing in FBG monitoring data of all sensors in borehole B6 in summer 2019 ( Figure A3). This gap is related to the breakage of the reinforced fiber optical cable connecting the borehole location to the FBG data acquisition system that happened during a debris flow event. Furthermore, sensor one of boreholes B4 and B6 contain incomplete data recordings because of malfunctioning FBG temperature sensors causing traveling reflection peaks that sometimes interfere with the reflected light peaks from the strain sensors because of similar wavelength. Figure A1. Upper plot: The colored lines show a visualization of the temperature corrected strain data of the individual 4.81 m base-length fiber Bragg grating (FBG) sensors installed in borehole B2. The strain data of each sensor is plotted in the center of the specific depth interval, and the vertical axis ticks correspond to anchoring depth of the FBG sensors. Positive strain means elongation and negative strain compaction of the sensor, and the strain scale is provided on the top left in microstrain (μɛ). Label a indicates a strain step in FBG-4 that was identified as a measurement artifact and corrected before the analysis. Label b indicates a significant rapid strain event that was measured in many sensors, coincides with an extreme rainfall event, and is referred to in the text. Lower plots: The vertical blue and cyan bars show the cumulative total precipitation data (per 24 hr) provided by MeteoSchweiz from the weather station "Bruchji" (Valais) located approximately 6 km away from our study site. The cyan colored bars indicate precipitation that presumably occurred as snow (i.e., at surface temperatures below 1°C at our study site [cf., Jennings et al., 2018]). Red time series show surface temperatures measured at the study site during the monitoring period. Figure A2. Upper plot: The colored lines show a visualization of the temperature corrected strain data of the individual 4.31 m base-length fiber Bragg grating (FBG) sensors installed in borehole B4. The strain data of each sensor is plotted in the center of the specific depth interval, and the vertical axis ticks correspond to anchoring depth of the FBG sensors. Positive strain means elongation and negative strain compaction of the sensor, and the strain scale is provided on the top left in microstrain (μɛ). Additionally, the elevation of the pressure head in the borehole is provided as a gray area in the background of the plot. Labels a and b indicate strain events referred to in the text. Lower plots: The vertical blue and cyan bars show the cumulative total precipitation data (per 24 hr) provided by MeteoSchweiz from the weather station "Bruchji" (Valais) located approximately 6 km away from our study site. The cyan colored bars indicate precipitation that presumably occurred as snow (i.e., at surface temperatures below 1°C at our study site [cf., Jennings et al., 2018]). Red time series show surface temperatures measured at the study site during the monitoring period.

Appendix B: Earthquake Events
In Figure B1 all earthquake events used for the comparison with the strain events measured in our research boreholes are characterized by means of earthquake magnitude, epicentral distance to the study site, and earthquake depth. We used all recorded earthquakes that occurred during our monitoring period with magnitudes greater than M1, epicentral distances of less than 50 km to our study site, and with depth between 0 and 15 km. The strain data of each sensor is plotted in the center of the specific depth interval, and the vertical axis ticks correspond to anchoring depth of the FBG sensors. Positive strain means elongation and negative strain compaction of the sensor, and the strain scale is provided on the top left in microstrain (μɛ). Additionally, the elevation of the pressure head in the borehole is provided as a gray area in the background of the plot. Label a indicates a strain event referred to in the text. Lower plots: The vertical blue and cyan bars show the cumulative total precipitation data (per 24 hr) provided by MeteoSchweiz from the weather station "Bruchji" (Valais) located approximately 6 km away from our study site. The cyan colored bars indicate precipitation that presumably occurred as snow (i.e., at surface temperatures below 1°C at our study site [cf., Jennings et al., 2018]). Red time series show surface temperatures measured at the study site during the monitoring period.

Data Availability Statement
Data are accessible under: https://doi.org/10.3929/ethz-b-000505871. Figure B1. Earthquake catalog used for the comparison with strain events detected in our research boreholes (data provided by the Swiss Seismological Service [SED]). The earthquake magnitude is provided versus the epicentral distance of the individual earthquake event from our study site. Depths information for the individual earthquake events are shown with the color scale.