Recharge response and kinematics of an unusual earthflow in Liechtenstein

Understanding landslide behavior over medium and long timescales is crucial for predicting landslide hazard and constructing accurate landscape evolution models. The behavior of landslides in soil that undergo periodic displacements, termed earthflows or compound soil slides, is especially difficult to forecast at these timescales. This is because velocities can increase by orders of magnitude over annual to decadal timescales, due to processes such as changing recharge conditions, erosion of the landslide toe, and retrogression of the landslide head. In this paper, we provide a detailed analysis of the Schlucher landslide, an unusual earthflow that is perched above the village of Malbun, Liechtenstein. This landslide had been displacing by 10 to 20 cm/year until 2015, when displacements on the order of 2 m/year occurred from 2016 to 2018. These large displacements damaged landslide mitigation measures, caused numerous surface deformation features, and threatened the local population downstream of the earthflow. This landslide has an unusually long monitoring record, with accurate displacement and climatic data available since 1983. We analyze this nearly 40-year monitoring time series to estimate recharge from snowmelt and rainfall, and its correlation with displacement. We also analyze recently collected, high-resolution surface and subsurface data in order to understand landslide response to recharge, landslide kinematics through time, and catastrophic failure potential. We find that interannual displacements can be explained with variations in recharge; however, periodic surges with recurrence times of tens of years must be explained by other mechanisms. In particular, recharge into the landslide during the recent acceleration (2016 to 2018) was not anomalously high. Instead, we argue that loss of internal strength is responsible for this recent acceleration period, and that this mechanism should be considered when forecasting the surge potential for certain earthflows and soil slides.


Introduction
Landslides in soil that undergo periodic, slow to rapid displacements are among the most difficult types of landslides to understand and predict. These landslides, termed earthflows and compound slides by Hungr et al. (2014), are expected to undergo periodic accelerations, during which the velocity of movement increases from cm/year to m/day or m/hour (Mackey and Roering 2011;Hungr et al. 2014). Therefore, the typical risk posed is to infrastructure, and risk to life is expected to be low, although catastrophic failure can occur in some rare cases (Fletcher et al. 2002;Aaron et al. 2017Aaron et al. , 2019. These landslides also represent an important sediment source in some landscapes (Keefer and Jonhnson 1983;Fletcher et al. 2002;Mackey and Roering 2011).
Understanding and predicting the motion of earthflows and compound soil slides require detailed knowledge of landslide response to recharge, as well as landslide kinematics. When analyzing these factors, it is important to select a timescale of analysis, due to the fact that time-dependent changes in landslide velocity must be accounted for (Mackey et al. 2009;Mackey and Roering 2011;Nereson and Finnegan 2018). In the present work, we consider the daily to decadal timescales, over which monitoring data is occasionally available (Nereson and Finnegan 2018).
At the annual timescale, many researchers have established relationships between landslide recharge, from rain and snowmelt, and displacement (Iverson and Major 1987;Matsuura et al. 2008;Schulz et al. 2009a, b;Coe 2012;Nereson and Finnegan 2018). This work has found that periods of acceleration are often related to times of high recharge, such as during spring snowmelt periods. There can be a substantial delay between recharge and displacement, on the order of weeks to months (Iverson and Major 1987). Based on these correlations, relationships have been established between historical records of recharge and displacement to forecast earthflow displacements for various climate scenarios (e.g., Coe 2012). Mechanisms other than variations in recharge have been hypothesized to control long-term velocity variations in earthflows and soil slides, including fluvial incision and availability of regolith in the source zone (e.g., Mackey and Roering 2011;Nereson and Finnegan 2018).
Earthflows and compound soil slides are distinct landslide types. This is mainly caused by differences in the percentage of plastic, clayey material, which leads to different deposit morphologies. As will be summarized below, both landslide types deform along discrete basal, internal, and lateral shear surfaces in response to pore-pressure changes, and can display similar kinematic modes. However, earthflows can undergo episodes of distributed internal shear strain, or failure along multiple shear surfaces, to produce a flow-like morphology that is absent in other types of soil slides (Hungr et al. 2001. Earthflow and compound soil slide kinematics can vary through time, and have been interpreted based on surface displacement measurements such as GPS, InSAR, and air photo analysis (Coe et al. 2003;Squarzoni et al. 2005;Baldi et al. 2008;Dewitte et al. 2008;Guerriero et al. 2014Guerriero et al. , 2017Clapuyt et al. 2017;Schulz et al. 2017). These analyses have revealed that these landslides are often subdivided into distinct kinematic elements, defined as sectors with consistent spatial and temporal kinematics, and/or kinematic zones, when bounding extensional and compressional features can be identified (Guerriero et al. 2014(Guerriero et al. , 2017Schulz et al. 2017). Local variations in the velocity of kinematic elements/ zones are controlled by factors such as spatial variations in snowmelt and slope angle (e.g., Delbridge et al. 2016;Guerriero et al. 2017). During acceleration periods, direct subsurface measurements are difficult to obtain, as shearing of boreholes can result in their becoming unusable. For earthflows, indirect subsurface Landslides measurements, as well as earthflow morphology, suggest that internal straining and flow-like motion can occur during surges (Keefer and Jonhnson 1983;Mainsant et al. 2012;Bertello et al. 2018).
The existence of distinct kinematic zones in earthflows, which undergo extension, compression, and/or shearing (Delbridge et al. 2016;Guerriero et al. 2017), implies that internal strains develop within the landslide body. These internal strains would activate the internal strength of the landslide material, which provides additional resistance to movement. In soil slides, internal strength can result from both frictional and cohesive strength of the landslide body, which can be reduced through time by processes such as grain crushing, loss of apparent cohesion through wetting, loss of true cohesion through breakage of grain cementation, and shearing of fine-grained soils from peak to residual strength (e.g., Fletcher et al. 2002;Augustesen et al. 2004;Sassa and Wang 2005;Stark et al. 2005).
This paper presents and analyzes the exceptional, long-term monitoring data of the Schlucher landslide, an unusual earthflow located in Liechtenstein. This landslide displayed anomalously high displacements between 2016 and 2018, prompting a detailed site investigation. Accurate displacement and climatic data are available since 1983, and a detailed site investigation, which included subsurface drilling, was conducted in 2017. We analyze this unique dataset to derive continuous recharge estimates, and examine the spatial displacement field and internal strain accumulation through time. Additionally, we perform finite element modelling to clarify the mechanisms governing the motion of the Schlucher landslide. These analyses enable the examination of whether interannual variations in recharge can explain the observed displacement trends, provide a basis for delineating the kinematic zones that compose the landslide and their persistence through time, and allow us to investigate the role of the internal strength of the landslide body in regulating landslide response to recharge.

Site overview
The Schlucher landslide is an unstable soil mass located in Liechtenstein. An overview of the study area is shown in Figs. 1 and 2. As can be seen in Fig. 1, the landslide is perched above the town of Malbun, which is a popular ski destination. Figure 2 shows that the landslide extends from rocky cliffs at the upstream end, to a bulging toe which is fronted by a large boulder and confined in a channel. A talus cone overprints the landside near the headscarp. Figure 3 shows cumulative landslide displacements since 1983, and demonstrates that the landslide has accumulated significant displacements since 2015. Figure 3 also shows that Pt. 13 has accumulated much less displacement than the downstream GPS points, indicating that this upslope portion is relatively inactive. These displacements have resulted in noticeable deformation features at the landslide surface, shown in Figs. 2 and 4. Starting at the top of the active portion of the landslide, a new backscarp can be seen on Fig. 4 a and d. This feature, which became active after 2015, separates the relatively inactive upper portion from the recently active lower portion. Alternating compressive and extensional features can be seen downstream of this backscarp (Fig. 4c, d, e). These are linear features which span the width of the earthflow, and the sense of movement is apparent from the deformed mitigation measures which have been pulled apart (Fig. 4d, c) or compressed. Within the lower portion of the landslide, an extensional feature can be seen which is perpendicular to the upstream features (Figs. 2 and 4b). The toe of the landslide is over-steepened, and features a large boulder (Fig. 2). The lateral margins of the landslides are apparent from shear deformation of a wood lined channel (Fig. 4f).
The Schlucher landslide was documented for the first time on an aerial image in 1946, where the landslide body is vegetated and shaped by fluvial erosion. In the same catchment, relatively frequent debris flows have been documented since 1944, as well as at least one small debris slide in 1981. Mitigation work to prevent infiltration into the landslide was undertaken from 1983 to 1985 and 2006 to 2008. As described above and shown in Fig. 4 c, d, and f, these surface mitigation measures have been heavily damaged due to landslide displacements which occurred from 2015 to 2018.
Four boreholes, the locations of which are shown in Fig. 2, were drilled into the landslide mass in October 2017. These boreholes were instrumented with a variety of sensors, described in the "Data" section, which include inclinometers and a coaxial cable (used to collect time domain reflectometry measurements). The borehole tubes were sheared off in spring 2018, after approximately 120 cm of landslide displacement. In borehole BL1, subsurface displacements were measured at depths of 7.8 meters below ground surface (mbgs) and 12.3 mbgs. In BL2, a primary rupture surface was measured at 15.3 mbgs, but numerous internal deformations were noted throughout the borehole depth. In BL3, displacements were measured at 12.4 mbgs depth, and in BL4, measured displacements occurred at 6.6 mbgs and 14.1 mbgs. In all boreholes, measured displacements occurred at discrete horizons, and were thus not distributed with depth (which would imply creep-like motion).
These subsurface measurements show that the landslide is approximately 15 m thick, and features internal sliding surfaces at some locations. The zone of the landslide that has accumulated significant displacements since 2015, which is shown in Fig. 2, has a planar area of approximately 10,700 m 2 . This leads to an estimated volume of 160,000 m 3 . Including the relatively inactive upslope areas of the landslide body, the volume of the landslide is approximately 210,000 m 3 .
The Schlucher landslide is located in a catchment composed of Hauptdolomite, Raibler Formation, and Arosa Zone of the Austroalpine nappes. The Raibler and Arosa Zones contain geotechnically problematic rock types, such as gypsum, rauhwacke, and red or green clay shales/sandstones. The landslide body itself comprises a heterogenous soil deposit composed of remobilized bedrock weathering products, forming layers dominated by different bedrock lithologies. During drilling, continuous core was extracted and visually described, and eight samples were taken and tested in the laboratory. BTG (2018) used this data to develop an interpreted stratigraphic cross section, a simplified version of which is shown in Fig. 5. The top most layer is primarily composed of blocky gypsum deposits, followed by layers dominated by red and green clay/siltstones layer which contain all sliding surfaces. Below these clay/siltstone-rich layers, a dolomite rich soil layer (Hauptdolomite) was encountered in BL 1 (14 m) and BL2 (15.4 m). The base above the bedrock contact (assumed to be Raibler Formation) is composed of dense glacial till with a block of gneiss detected in BL3 (18.3 m). The lithological layering, a thin organic layer observed in borehole BL1 (7.5 m) and a tree branch in BL2 (15.8 m) (both corresponding to measured displacement depths) suggest that the landslide material may be composed of colluvial deposits formed during multiple postglacial events.
The visual description of the core in BL1 showed that the landslide mass at this location is primarily composed of silty sandy gravel, although thin (0.3 to 0.7 m) layers of silt were noted. The depth of measured displacements in BL1 corresponds to two of these silt layers. In BL2, silt, clay, and fine sand units were observed from 0 to 7 mbgs, and 10.2 to 13.0 mbgs, and a clayey, silty sand and gravel unit is in between these fine-grained layers. The measured displacements occurred in a silt layer. In BL3, the core is mapped as a sand, with varying proportions of silt, clay, and gravel, and the measured displacements occur in a silty sand layer. In BL4, most units were mapped as sand, with varying proportions of silt and gravel, although a sandy, silty gravel layer was noted from 8.5 to 11.8 mbgs. Measured displacements occurred in fine sand units with a high proportion of silt. A fine-grained colluvial layer was mapped below the depth of displacements in BL1 and BL2. A gravely sand layer, interpreted to be a moraine deposit, was mapped below the locations of the deep rupture surface in BL3 and BL4.
These visual interpretations are consistent with laboratory material testing, which showed that the material sampled from the boreholes ranged from clayey sand to clayey gravel. The laboratory testing also showed that the percentage of fines ranged from 17 to 47%, which, based on Atterberg limit tests and the USCS classification system, is classified as low plastic silt or low plastic clay (Table 1). Standard penetration tests (SPT) were conducted at 5 m intervals in BL2, BL3, and BL4. SPT N values, shown in Fig. 6, indicate that the material is medium dense to dense.
Measured surface and subsurface displacements, described in detail below, show that the landslide displays intermittent movement, along multiple discrete rupture surfaces. Based on this, the Schlucher landslide can be classified as an earthflow (using the terminology of Hungr et al. (2014)). It should be noted that earthflows are typically composed of soils with a higher proportion of plastic, fine material than that present in the Schlucher landslide. Therefore, this landslide does not display every characteristic of this class of mass movement and is transitional in nature to a compound debris slide.

Data
During the site investigation, both surface and subsurface data were collected. Surface data includes multiple orthophotos, highresolution digital surface models of the site, oblique timelapse camera imagery, measurements of surface displacements, and climatic measurements from a nearby weather station. Orthophotos were taken in 2015, 2017, and 2018, with 50 cm, 1.5 cm, and 1.5 cm resolution, respectively. The 2017 and 2018 photos were captured using an autonomous drone, and were used to derive high-resolution digital surface models (25 cm pixel spacing). The timelapse camera was installed in March 2018, and overlooks the entire Schlucher catchment. Twelve images of the catchment were recorded every day.
Periodic manual geodetic and GPS measurements of marked points on the landslide have been taken annually or semi-anually since 1983. During 2016, the interval between point measurements was decreased to once a week. Continuous GPS measurements have been taken from three permanent stations since October 2017. Interruptions in this data have been caused by avalanches in the winter of 2019, as well as relocations of the base station.
The location of the weather station used in the present analysis is shown in Fig. 1. This station has a heated rain gauge, and has collected precipitation data since 1972. Continuous temperature data is available from this station since 1998. Since 1998, periodic measurements of snow water equivalent are made during the winter at the location of the weather station.
The subsurface data, a subset of which was presented in the "Site overview" section, includes the following information gathered from four boreholes which were drilled in October 2017 ( Fig. 2): & Visual description of the core recovered during drilling. & Photos of the core material recovered during drilling. & Pore-water pressure, fluid electrical conductivity, and groundwater temperature.  Fig. 4 were taken (green dots). The sense of displacement of the deformation features is interpreted based on the deformed surface mitigation structure. Bottom: detailed overview of the most active portion of the landslide, showing the locations of the four boreholes (black dots), three permanent GPS stations (magenta dots), and four periodically measured GPS points (red dots). The section line A-A′, which is used for the stratigraphic section and finite element analysis, is shown on the top Figure   Fig. 3 Cumulative displacement of the Schlucher landslide, measured at five GPS points since 1983. The locations of the GPS points are shown in Fig. 2 Pore-water pressure, fluid electrical conductivity, and groundwater temperature were initially measured using three downhole sensors/loggers installed in boreholes BL2, BL3, and BL4 (Fig. 2). Each of these boreholes contained one or two 1 m long screened inclinometer tube intervals (with geotextile), located at the bottom of the borehole below the depth of the deepest discrete shear surface (Fig. 5). It can be assumed that initially, the porepressure sensors recorded pressure levels from these screened intervals.
Landslide displacements, which accumulated in April and May 2018, resulted in these sensors becoming unusable. After this time, new pressure sensors were installed at shallower depths, and are likely sampling the pore-water pressures from the upper shear surfaces, i.e., the location where the inclinometer tubes were destroyed. BTG (2018) investigated these sheared inclinometer tubes using a downhole camera. This analysis revealed that the sheared borehole tubes appear to be blocked at varying depths within the landslide mass, and that the inclinometer tubes are broken at discrete depths. BL1 appears to be sealed by 3 stones at a depth of 7.2 mbgs, and no breaks in the tube were noted. Therefore, the measured water level in BL1 may be representative of that at 7.2 mbgs. In BL2, the camera could be lowered to 14.6 m depth, where the tube is filled with fine-grained material. Breaks in the tube were noted at 4.7 to 4.9 m, 5.0 m, and 12.1 to 12.4 m. It is difficult to interpret the depth where the measurement is representative; however, water levels in BL2 appear to stay below 5.0 m. Therefore, the measured pressures may be representative of 12 to 14 m, although water may enter the tube at 4.7 m during infiltration periods. Pumping tests revealed that measurements in BL3 were not reliable. In BL4, the camera could only be lowered to a depth of 6.1 m, where the inclinometer pipe was blocked. Pressure measurements in this tube are likely representative of pressures at 6.1 m depth. Table 2 shows the time interval of the three pore-water pressure sensors.
Due to the differing temporal resolution of the available monitoring data, three separate time periods are analyzed in the present work. These are: & October 2017 to August 2019, when timelapse camera imagery during snowmelt periods, continuous temperature data, subsurface pore-pressure data, and surface and subsurface displacement measurements are available.
& January 1998 to August 2019, when continuous temperature data and measured snow water equivalent are available.
& September 1983 to June 2015, when annual or semi-annual surface displacement data is available. During this time period, no continuous temperature data is available, and measured snow water equivalent is only available from 1998.

Quantification of recharge
Understanding recharge conditions into the Schlucher landslide is important to understand the drivers of motion and their variance over time. We aim to test whether variations in recharge can account for measured variations in displacement. Thus, we estimate recharge from both rainfall and snowmelt based on temperature and precipitation data from a weather station that contains a heated rain gauge and is located at the catchment outlet ( Fig. 1).
Liquid precipitation is measured directly by the rain gauge; however, time series of snowmelt must be estimated. We use the degree day method (e.g., Rango and Martinec 1995) to estimate daily snowmelt depth, using Eq. 1.

Original Paper
Landslides where M is the snowmelt depth (cm), a is the degree day factor cm°C day , and T is the number of degree days (°C day). A degree day is defined as a daily mean temperature of one degree above a threshold temperature, taken as 1°C in the present analysis. Therefore, a daily mean temperature of 4°C would result in a value of 3 . The degree day factor (in Eq. 1) is known to vary over the course of the year, with low values during the winter and high values in April, May, and June. We parameterized the time varying degree day factor based on the values given in Rango and Martinec (1995) and a simple time varying function. During winter months until March 1st, a value of 0.07 cm°C −1 day −1 was used, followed by a linear increase in the value of the degree day factor to 0.4 cm°C −1 day −1 by April 1. We verified our degree day factor function by comparing measured and modelled snow water equivalents for the period from 1998 to present.
To estimate total snow water equivalent available for snowmelt infiltration, precipitation (in mm) that occurred between December and March when measured temperatures were below 2°C was summed. Thus, the procedure estimates a continuous time series of snow water equivalent accumulation and melting, which is used for recharge estimation.
Precipitation data is available throughout the entire monitoring period; however, available temperature data is different for different time periods. Additionally, a timelapse camera was installed in 2017, which allows further refinement of the snowmelt estimates in the catchment of the Schlucher landslide. Therefore, three different time periods were considered, and recharge estimated differently for each. These are: 1. 2018 to 2019, when timelapse camera data is available. 2. 1998 to 2019, when continuous temperature data is available. 3. 1983 to 2019, when daily temperature measurements, taken between 07:00 and 08:00 are available.
The methods used for each time period are summarized in the following sections.

to present
For the analysis of the two snowmelt periods from April to May 2018 and March to June 2019, timelapse camera imagery is available. The catchment area was subdivided into three zones, and the mean daily temperature at the median elevation in each zone was determined by linearly extrapolating the temperature from the weather station (Fig. 1). The slope of this linear extrapolation was determined by fitting a linear temperature function to three nearby weather stations installed at different elevations. The three zones used to estimate snowmelt are shown in Fig. 7. The spatial area covered by snow in each of the three zones was assessed individually based on the timelapse camera imagery, by estimating the ratio of snow covered pixels to snow free pixels. This only provides an approximate estimate of the true area covered by snow, as each pixel represents a different spatial area due to topographic effects. To partially address this limitation, each of the snowmelt zones were subdivided into facets of approximately uniform slope angle, and a weighted sum of these facets was used to estimate the average snowmelt in each zone, with weights based on the true area of each facet determined with a high accuracy digital elevation model.
Finally, a modified degree day method that accounts for the area covered by snow was applied to estimate snowmelt in each of the three zones: and c (%) is the percent of the area covered by snow, estimated based on the timelapse camera images and corrected based on a digital elevation model.  For this period, a continuously evolving snow water equivalent profile was estimated. Precipitation that occurred when temperatures were below 2°C was added to available snow water equivalent, and snowmelt estimated from the degree day method (Eq. 1) reduced available snow water equivalent and was treated as recharge. During this time period, evapotranspiration was accounted for using the Turc method (Turc, 1961).

to present
The only available temperature measurements that extend back to 1983 are periodic measurements taken daily between 07:00 and 08:00 during the months of January to April. Therefore, to estimate total snow water equivalent from 1983 to present, precipitation (in mm) that occurred between December and March when measured morning temperatures were below 2°C was summed. The degree day method was used to estimate snowmelt during the winter. The difference between the cumulative precipitation and estimated snowmelt was then used as the snow water equivalent available for recharge during the snowmelt period. Snowmelt is assumed to have occurred by March 31 st each year; however, the results of the analysis are not sensitive to this assumption due to the frequency of displacement measurements, summarized below.

Displacement and strain calculation
A number of different methods, including periodic GPS measurements (semi-yearly to yearly), continuous GPS measurements, and digital image correlation, were used to quantify landslide displacements and kinematics. For the period from 1983 to 2019, measurements from five periodic GPS points, shown in Fig. 2, were used to analyze landslide displacements.
To assess the relationship between displacement and recharge using the periodic GPS measurements, the displacement that occurs between successive measurements (termed "incremental displacement") was compared to estimated recharge from liquid precipitation and snowmelt that occurs in the same time period (termed "incremental recharge"). Thus, if 160 cm of recharge (combination of rainfall and snowmelt) is estimated to have occurred between successive displacement measurements, during which 25 cm of movement occurred, incremental recharge would be 160 cm, and incremental displacement would be 25 cm. In October 2017, three permanent GPS stations, which provide continuous displacement monitoring, were installed on the landslide body (GPS1, GPS2, and GPS3).
To quantitatively investigate landslide kinematics, strains within the landslide body were estimated based on two methods. For the period spanning from 1984 to 2019, annual 1D strains were estimated between each of the GPS points shown in Fig. 2. This was done by taking the ratio of the change in the distance between the two points to the original distance. A similar procedure was used by Guerriero et al. (2017) to estimate 2D strains from GPS measurements.
A more detailed strain estimate was made for the period from 2017 to 2018 by performing digital image correlation (DIC) on two high-resolution orthophotos, using the algorithm described in Bickel et al. (2018). Briefly, this algorithm tracks displacements of groups of pixels that can be identified in successive images, a process which results in a 2D horizontal displacement vector field with a high spatial resolution.
This displacement field was then further processed in order to derive a 2D horizontal strain field for this time period. This was done by dividing the digital image correlation results into quadrilateral elements, with the initial state corresponding to the grid of points derived from the 2017 imagery, and the deformed state corresponding to the locations of the same points in 2018 (see inset on Fig. 12 for a visual representation of this procedure). We then derive continuous strain fields, using the Green definition of the strain tensor (e.g., Bechtel and Lowe 2015). The Green strain tensor (denoted as E below) eliminates rigid body rotations, and is given as: where F is the deformation gradient tensor, and I is the identity tensor.
We parameterize F using the standard finite element interpolation algorithm and quadrilateral elements with linear shape functions (Reddy 2019). We then estimate the Green strain tensor at the centroid of each of the quadrilateral elements. This tensor corresponds to a coordinate system with the x-axis aligned east-west, and the y-axis aligned north-south. In order to assist with the interpretation of the strain tensor, we further transform the estimated Green strain tensor into the principle strain tensor at the centroid of each element, using standard formulas (e.g., Bechtel and Lowe 2015).

Finite element modelling
We perform simplified, 2D finite element modelling using RS2 (Rocscience 2015) in order to provide mechanical insights into the behavior of the Schlucher landslide. In particular, the numerical modelling was used to test whether internal deformation patterns interpreted from surface displacement data are compatible with the interpreted rupture surface, as well as the implications of these internal deformation features for landslide displacement magnitudes. We use the geometry shown in Fig. 5, and simplified the material types of the landslide into one material that represents the colluvium, and one higher strength unit for all material below the rupture surface. The model setup is shown in Fig. 7 Discrete snowmelt zones. Zone 1 has a median elevation of 1740, zone 2 has a median elevation of 1840, and zone 3 has a median elevation of 1950 Fig. 15a. We specified the location of the rupture surface based on the interpretation of the borehole data (Fig. 5), and treated it as a low strength joint.
We selected material properties that were compatible with the results of the site investigation, and these properties are shown in Table 3. Elastic properties (Young's modulus and Poisson's ratio) were selected to be representative of coarse material with a significant fines component. Additionally, a friction angle of 30°, typical for medium dense coarse-grained soil, was selected. As mentioned previously, the goal of the finite element modelling is to test the influence of the internal strength of the earthflow body. We therefore tested two scenarios, one with high cohesion in the earthflow body, and one with low cohesion. For the high internal strength scenario, a value of 50 KPa was used, and for the low internal strength scenario, a value of 10 KPa was used. We note that mechanical testing of the material that composes the earthflow body is not available at present, so model results are only able to show the results of various material strength assumptions, which are compatible with available core descriptions and index properties. Preliminary analysis indicated that the limit equilibrium apparent basal friction angle for this landslide is 18°, and displacement monitoring (described in more detail below) indicates that the landslide is continuously creeping. Thus, we hypothesize that the apparent basal friction angle (which includes the influence of pore pressure) remains near 18°, but may be lowered during periods of high infiltration. We therefore used an apparent basal friction angle of 16°in our simulations to simulate periods of high recharge. This apparent friction angle corresponds to an approximately 3 m rise in the groundwater table in the landslide, which is consistent with, but somewhat higher than, pressures measured within the earthflow body (described below). Thus, our simulations likely represent conservatively high pore-pressure conditions.
We validate the results of our finite element modelling by comparing simulated to measured displacement patterns, as well as simulated to estimated kinematic zones. In particular, we compare the 3D inclination of displacement vectors, measured during the recent acceleration period (2015 to 2018) with the inclination of simulated displacement vectors. We also compare simulated zones of internal yielding to those determined from the displacement and strain calculations. Figure 8 shows the estimated snow water equivalent compared to measured values since 1998. Modelled and measured values are similar, suggesting that the time varying degree day factor adequately captures snowmelt behavior at this site. Additionally, Fig. 8 shows measured and modelled maximum snow water equivalent for each year since 1983. The two are again similar, indicating that snow water equivalent estimates since 1983 capture the variance of this quantity from year to year.  Figure 9a shows that the GPS stations displaced between approximately 2 and 2.5 m, and that there is a delay between snowmelt and displacement. Figure 9b shows that the displace-  (Fig. 10a), although a small displacement trend can be observed. In both 2018 and 2019, snowmelt occurs at all elevations simultaneously (Figs. 9a and 10a). This suggests that the snowmelt estimated at the weather station shown in Fig. 1 is a good proxy for landslide recharge during snowmelt periods, and this is used to understand landslide response to recharge events from time periods before timelapse camera imagery was available. Figure 10a and b show that displacements in 2019 were about two orders of magnitude less than in 2018, and peak snowmelt values in 2019 were approximately half of those measured in 2018 (this observation is discussed further below). Figure 10b shows water levels above the rupture plane in BL2 compared to measured GPS displacements. This plot shows that the minor displacements observed in 2019 correlate with peaks in pore-water pressure, although these pressures only correspond to an increase of less than 1 m of water in BL2. Displacements in September/October of 2018 correspond with water level fluctuations in BL4 (compare Fig. 10b and c); however, the displacement time series show a constant displacement trend, with no velocity fluctuations which correlate with pressure fluctuations. This may indicate that BL4 is measuring pressures in a perched aquifer within the landslide body, and does not represent pressures at the rupture surface. Figure 10c shows a direct temporal correlation between the porepressure levels recorded in the sheared-off boreholes BL1, BL2, and BL4 and rainfall and snowmelt events, with a time shift between the onset of rainfall and pore pressure of only a few days. The measurement depth in BL1, BL2, and BL4 is 7.2, 14.6, and 6.1 m, respectively, and the highest amplitude pore-pressure responses occur at the shallowest measurement depth (BL4), and appear to get more diffuse with depth. It should be noted that these measurement points are offset in both depth and location on the slope (Fig. 2), so 2D effects may also influence the relative amplitude of the measurements.

to present
1998 to present Figure 11a shows incremental recharge, adjusted for evapotranspiration, compared to incremental displacement between 1998 and 2019. The highest incremental recharge occurred in 1999, and during this year, the highest pre-2015 displacements are measured (~1 m). During the following 2 years, recharge remains high and displacements are relatively low. Figure 11a also shows that the displacements measured from 2016 to 2018 are orders of magnitude greater than in the preceding years, despite similar quantities of incremental recharge, and that the displacements measured in 2019 were an order of magnitude less than those measured in 2017 and 2018, despite similar recharge conditions. Figure 11b and c show incremental displacement as a function of incremental recharge for Pt. 1, resampled so that the spacing between measurements is approximately 1 year (thus incremental displacements approximately correspond to yearly landslide velocities). Four outliers can be seen in Fig. 11b, and these outliers are removed in Fig. 11c. Figure 11c shows that landslide sensitivity to recharge has varied on decadal timescales, and is lowest from 2006 to 2014, and highest from 1996 to 2006. This observation is further discussed in the "Discussion" section. Figure 11d shows incremental recharge compared to incremental displacement for 1984 to 2015 for five locations. In 1996 and 1999, high incremental recharge is estimated, and correspondingly high incremental displacements were measured. However, in 2000 and 2001, recharge is still high but corresponding incremental displacements are low. These observations are discussed further once the results of the landslide strain analysis and finite element modelling are presented.  Fig. 13. The raw vector field (Fig. 12) shows that the upper part of the landslide has accumulated more displacement than the toe between 2017 and 2018. The orientations of the vectors also show that the landslide is bending during displacement, as the upper vectors are orientated WSW, whereas the vectors at the toe are oriented SW. The inset of Fig. 12 provides an example visualization of the results of the strain analysis, and shows that the landslide undergoes both a rigid body translation and internal straining in order to move from the original (blue rectangles on Fig. 12) to deformed state (red polygons on Fig. 12). By comparing this inset to the results of the strain analysis shown in Fig. 13a, it can be seen that the rectangles that underwent the most strain are highlighted as a zone of high negative strain. Figure 13a and b show that strains within the landslide body are accommodated along a number of prominent linear features. The toe of the landslide is undergoing compression (Fig. 13a), with high compressive strains noted at the location of the large boulder (Fig. 2). Upslope of the boulder, a prominent compression feature is noted (Fig. 13a), corresponding to the location where the fast moving upper portion of the landslide is compressing the slower moving lower portion (Fig. 12). The upper portion is then further subdivided into three compartments by prominent extensional and compressional feature (Fig. 13a). Figure 13 b shows that the toe of the landslide is undergoing extension along a curved linear feature, which corresponds to the vectors in Fig. 12 that are splitting around the large boulder at the toe.

September 1983 to June 2015
The results of the strain estimates since 1983 are shown in Fig. 14, and are plotted as cumulative strain though time. As can be seen in Fig. 14, the strain rate stays at a low constant value for most years in the monitoring period. This pattern is interrupted in 1999, and between 2015 and 2018, when significant strains are estimated between GPS points 1, 7, and 13. These years correspond to years where the landslide underwent significant displacement (Fig. 11). Figure 14 shows that strains between Pt. 7 and Pt. 8 are positive (extension) during years where displacements are limited, but switch to negative (compression) during years with large displacements, potentially indicating a change in kinematics during years with high displacement.
We interpret these results to develop the kinematic model shown in Fig. 13c. We subdivide the landslide into two major zones, and each zone has separate sub-compartments. Zone 1 is bounded by the main scarp at the head of the active landside compartment, and by a zone of compression (Fig. 13a) at the front. Zone 1 has three subzones, separated by extensional and compressive features (Fig. 13a). Zone 2 is bounded by the zone of compression with zone 1, and the toe of the landslide, which is being compressed into and against the sides of the channel (Fig. 13a). There is a subzone within zone 2, apparent in Fig. 13b, where a portion of zone 2 is moving down channel. Fig. 10 a Estimated snowmelt and cumulative 3D displacement of the landslide since June 24/2018. b 3D GPS displacement for GPS 3 compared to pore pressure in BL2. Note that the gaps in the pore-pressure data are due to the water level dropping below the level of the pore-pressure sensor. c Water level measured in BL1, BL2, and BL4 compared to estimated recharge. The locations of the boreholes and GPS stations are shown in Fig. 2  Fig. 11 a Incremental recharge vs incremental displacement for the period from 1998 to 2019. The legend for these points is given in panel d. b Incremental recharge vs incremental displacement for the period from 1984 to 2019 for Pt 1. c Data from plot B, with outliers removed. d Incremental recharge and incremental displacement for GPS measurements from 1984 to 2015. Evapotranspiration is not accounted for in the estimate of recharge. In all panels, the GPS measurements have been resampled so that there is approximately 1 year between each measurement Figure 14 shows that the kinematic model determined from the digital image correlation is consistent with measured strains since 1983. In particular, Fig. 14 shows that strains between points 1 and 3, and between points 7 and 8 are near zero, indicating that significant internal strains have not been accumulating within zones 1 and 2 during the nearly 40-year monitoring period. Conversely, compressional strains have been developing between points 1 and 7, and extensional strains between points 7 and 13, consistent with the kinematic zonation shown in Fig. 13c and indicating strain accumulation between zones 1 and 2. Other researchers have found constant earthflow kinematics through time (Guerriero et al. 2014;Schulz et al. 2017).
All of the extension and compression features have been verified in the field, as shown in Figs. 2 and 4. Figure 2 shows a number of lineaments that span the width of the landslide, and correspond to the area of compressive strain that separates zone 1 and zone 2 (Fig. 13c). Figure 4 a and d show the main scarp at the head of the landslide, which is highlighted as a zone of extension in Fig. 13a. Figure 4c and d show a surface drainage system undergoing extension, at the location corresponding to the boundary between zone 1a and zone1b in Fig. 13c. Figure 4b shows a crack undergoing extension, located at the boundary between zone 2a and zone 2b in Fig. 13c. The shear strains accumulating at the landslide boundary are highlighted in Fig. 13a, and have damaged a surface drainage system (Fig. 4f).
A remarkable feature of the displacement field is that the large boulder located at the edge of the deposit (labelled "boulder" in Fig. 13c) appears to strongly influence landslide displacements and kinematics. Figure 12, as well as DIC performed on oblique imagery (using two images taken from the same angle as Fig. 4a), shows that the landslide is splitting around the boulder, and the location of the scarp that bounds subzone 2a is controlled by the boulder location. This boulder appears to be acting as a high strength unit that forces the surrounding lower strength material to deform around it.

Finite element modelling
The results of the finite element modelling are shown in Fig. 15. Figure 15a shows the simulated deformation vectors compared to the measured deformation vectors along the same cross section.  Additionally, Figure 15b shows the boundaries of interpreted kinematic zones compared to modelled locations of internal yielding. As can be seen, the correspondence between measured and modelled results is high.
When high internal strength was used, limited displacements of the landslide were simulated, even though the entire basal rupture plane yields (Fig. 15e). As shown in Fig. 15c, limited internal yielding of the earthflow body occurred. Conversely, when low internal strength was used (Fig. 15b, d), significant internal yielding occurred in areas where the inclination of the basal rupture plane changes. These simulations never converged, indicating that acceleration of the mass was ongoing in these simulations. As shown in Fig. 15b, the location of zones that experience significant internal yielding of elements corresponds well to the kinematic zones interpreted from the DIC analysis (Fig. 13).
Thus, the finite element analysis indicates that the Schlucher landslide exhibits two displacement regimes. When the earthflow is destabilized (i.e., the limit equilibrium factor of safety drops below 1), and the internal strength of the earthflow body is not overcome, then landslide displacements are simulated to be minimal. Conversely, when the internal strength of the landslide is exceeded, the landslide undergoes large displacements, accumulating significant internal strains. Another interesting observation from these model results is that, after the landslide has displaced some distance, the landslide will stabilize unless a new zone of internal failure forms. This can be seen in Fig. 15b, where the two white triangles indicate zones of internal yielding that only formed after the landslide had accumulated greater than 2 m's of displacement. Had these zones of internal yielding not formed (due to strength heterogeneity within the earthflow body), the landslide would have stabilized.

Relationship between recharge and displacement
The long-term record of recharge and displacement available for the Schlucher landslide demonstrates a complex and time varying response of the landslide to recharge. Starting with the recent, high-resolution data, Fig. 9a and b show that, in 2018, snowmelt is simulated to occur in two pulses, and two corresponding pulses in landslide acceleration are noted (Fig. 9b). There is an approximately 25-day delay between the onset of snowmelt and landslide acceleration, as shown in Fig. 9b. When considering a time lag of 25 days, the first acceleration period does not match details of the recharge function but the second acceleration period coincides very well with the simulated recharge from snowmelt (Fig. 9b). This could indicate that during the first recharge period, the heterogeneous landslide body re-saturates and during the second recharge period, the landslide velocity increases linearly with the amount of recharge. Therefore, we hypothesize that in 2018, it took 25 days for recharge to result in a significant pore-pressure increase at the depth of the main sliding surface, mainly delayed by the time necessary for infiltration in the unsaturated zone and pore-pressure diffusion to the deep shear surface. Response times of weeks to months have been noted in other studies of earthflows (Iverson and Major 1987;Handwerger et al. 2013).
In 2019, measured displacements were over one order of magnitude less than those measured from 2016 to 2018, an observation discussed in more detail below. Two minor displacement events occurred from April 22 to May 01 and from May 24 to June 29, and correlate with peaks in the measured pore-pressure levels of BL2 ( Fig. 10a and b). As discussed above, the pressure levels in BL2 are likely representative of groundwater pressures at the depth of the main shear surface (12 to 14 m). Contrary to the snowmelt period in 2018 (Fig. 10), no clear correlation between these minor accelerations and the snowmelt offset by 25 days is apparent.
Based on the observations detailed above, we assume that the landslide displacements that occurred in April and May 2018, as well as the minor displacements from April to June 2019, were caused by snowmelt infiltration increasing pore pressures in the main rupture planes of the landslide within a period of about a month. This detailed understanding of landslide response to recharge will be used to understand landslide dynamics since 1983, as shown in Fig. 11.
The recharge estimates shown in Fig. 11a were made using continuous temperature measurements, and show that recharge from 2016 to 2018 was not anomalously high despite significant displacements. Additionally, comparing Fig. 11b and c show that, if the four outlier years are removed, landslide displacements Our analysis of recharge and displacement has revealed a complex relationship between the two that changes over time. The landslide appears to have a background displacement rate, on the order of 5 to 20 cm per year (Fig. 11c). This suggests that, similar to many earthflows (e.g. Mackey and Roering 2011;Nereson and Finnegan 2018), the Schlucher landslide is in a critical state, and therefore sensitive to variations in degree of saturation and pore-water pressure. Superimposed on this background displacement rate are annual and decadal movement cycles. Annual cycles appear to be well explained by variation in recharge; however the decadal movement cycles must be controlled by a different driver. Another crucial observation from the analysis of this detailed time series is that the exceptional displacements observed in 1999, and from 2016 to 2018, labelled as outliers in Fig. 11b, cannot be explained by recharge alone, as Fig. 11a and d show that recharge conditions from these years are not exceptional compared to other years. This observation is further discussed in the following section.
Relationship between internal deformation and displacement As discussed above, large displacements that occurred in 1999 and from 2015 to 2018 cannot be completely explained by variations in Fig. 15 Results of finite element modelling of the Schlucher landslide. a Model setup and calculated displacement vectors compared to measured for the low internal strength simulations. b Yielded elements for low internal strength simulations. The lines and arrows indicate the location and direction (compression vs. extension) of the kinematic zone boundaries determined from the DIC and strain analysis (Fig. 13). The two white triangles denote zones of internal yielding that only formed after significant displacement (> 2 m). c Yielded elements for high internal strength. d Total displacement for low internal strength. Note that these simulations never converged, and displacements are ongoing. e Total displacements for high internal strength simulations, following model convergence recharge. Figures 13 and 14 show that, when large displacements occur, significant internal strains are accumulated within the landslide body. This is supported by the finite element analysis, which shows that, for significant landslide displacements to occur, the internal strength of the landslide body must be overcome, and when this occurs, zones of yielding result in places where boundaries between kinematic zones are observed (Fig. 15). We interpret this to indicate that the Schlucher landslide has two displacement regimes, which govern landslide sensitivity to recharge, and are regulated by the internal strength of the landslide body at locations where the inclination of the rupture surface changes, analogous to the strength loss behavior exhibited by many compound rockslides (Glastonbury and Fell 2010). As noted by Guerriero et al. (2014), these locations are determined by channel geometry (for areas of the earthflow located outside the source area), and their location remains stable through time.
When the internal strength of the landslide body is not overcome, as appears to occur in the majority of years in our monitoring record, the landslide is much less sensitive to recharge than when internal yielding can occur. This is because internal strength of the landslide body prevents internal deformation from occurring, and internal deformation is necessary for accumulating significant displacements (Fig. 15). During these years, we hypothesize the internal strength of the landslide body could be degraded by mechanisms including grain crushing, shearing of dense granular material, and shearing of low plastic fine-grained layers. These mechanisms would preferentially occur at areas of high stress concentration, which are the same areas where internal yielding later develops, leading to a heterogenous distribution of internal strength within the landslide body. We note that this behavior would not be exhibited by most earthflows, which tend to be composed of ductile, high plasticity fine-grained material.
Once the internal strength has been sufficiently degraded, snowmelt and/or rainfall can trigger significant displacements, such as those observed from 2015 to 2018. This would lead to internal yielding of the weakened material located in zones of the earthflow where the basal rupture surface changes, as shown in Fig. 15b. Additionally, during years where significant displacements occur, strains between Pt. 7 and Pt. 8 switch from extension to compression (Fig. 14). Once enough displacement has accumulated, a new zone of internal failure must develop in the landslide body (this is shown by the white triangles in Fig. 15b), which would switch the landslide displacement regime back to that of low sensitivity, until the internal strength of the earthflow body can once again be overcome. This switch of displacement regimes may have occurred in 2019, when snowmelt was approximately half that observed in 2018, and displacements were limited. This preliminary finding should be further tested with longer time series.
Thus, our results suggest a further mechanism that can lead to velocity variations in earthflows and soil slides, in addition to fluvial incision, availability of regolith, and changes in surface moisture balance (Nereson and Finnegan 2018). Given that none of these three factors could have resulted in the altered state of the landslide body observed from 2015 to 2018, we suggest loss of internal strength as an additional mechanism that can lead to periodic landslide acceleration, a mechanism also highlighted by Fletcher et al. (2002). This would only effect earthflows with a variably inclined rupture surfaces, and brittle internal strength, provided either by continuous lithological layers of low plastic fines and/or dense granular materials.
The data collected during the site investigation provides evidence that internal deformation has occurred from 2015 to 2018, and that both loosening of material, and periods of high landslide sensitivity to recharge have occurred in the past. Evidence of internal deformation is available from surface displacement measurements, summarized in Figs. 13 and 14, and can also be interpreted from subsurface measurements. BL2 is located on the compressional feature that separates kinematic zones 1 and 2 in Fig. 13, and numerous displacements of the borehole casing were noted at various depths overlying the primary rupture plane. Additionally, multiple sliding planes were noted in BL1 and BL4, which could indicate internal shearing of the earthflow material. Further, with the exception of two measurements, the SPT N values systematically decrease towards the toe of the landslide (Fig. 6), indicating loosening of the material with increasing displacement. Finally, aerial imagery captured in 1946 and 1984 indicates that the boulder at the front of the landslide (Fig. 2) has displaced approximately 43 m during this 38-year period. If a conservative assumption of 0.3 m of displacement per year during low sensitivity years is assumed (Fig. 11c), displacement of the boulder would be underpredicted by about 32 m, giving an indication that the high sensitivity landslide regime has occurred in the past.
The analyses presented here have implications for forecasting the behavior of soil landslides that undergo periodic displacements over decadal timescales based on future climate scenarios. For the present case, fitting a recharge-displacement function to the high quality, 32-year long record of recharge and displacement prior to 2015 would have underpredicted the displacements of the Schlucher landslide by one to two orders of magnitude during the period from 2015 to 2018. As summarized above, we argue that a change in the internal strength conditions of the landslide, caused by a loss of internal strength, changed the force balance of the Schlucher landslide, leading to higher sensitivity to pore-pressure changes. Considering this mechanism when forecasting landslide motion requires mechanistic models that consider the force balance within the landslide.

Conclusions
Earthflows and compound soil slides represent a major sediment source in many landscapes, and displacements during surges can damage infrastructure. Thus, predicting velocities at the annual and decadal timescales is crucial for understanding landscape evolution and future landslide risk. We used the exceptionally detailed monitoring data available for the Schlucher landslide to assess the drivers of motion during the 36-year monitoring history. We show that the displacements of 5-10 meters per year measured in 1999 and from 2015 to 2018 are unique when compared to the entire available displacement history, and were likely driven by a loss of internal strength. Our results suggest that loss of internal strength may be a mechanism that can lead to surging in certain earthflows, and that this mechanism is difficult to account for based on empirical fitting of recharge-displacement time series.
Liechtensteinische Landesverwaltung. We are grateful to Stephan Wohlwend and Emanuel Banzer of Abteilung Naturgefahren, Amt f ü r B e v ö l k e r u n g s s c h u t z ( A B S ) , L i e c h t e n s t e i n i s c h e Landesverwaltung as well as Reto Thöny, Daniel Figi and Matthias Vogler of Büro für Technische Geologie AG for collecting and providing the data used in this study. Clément Roques kindly provided the script used to estimate evapotranspiration. Two anonymous reviewers provided constructive and thorough reviews, which greatly improved the manuscript.

Funding
Open Access funding provided by ETH Zurich.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.