Forecasting landslides by taking their temperature: Combining mathematical and experimental modeling with field data on the case of the El Forn landslide (Andorra)

In this study we are suggesting a temperature-based modeling approach for deep-seated landslides, validated through combined field monitoring and experimental testing. The Silurian shales of the shear band of El Forn landslide (Andorra) have been characterized through thermal and rate controlled triaxial tests, thereby calibrating a mathematical model that is used to monitor the behavior of deep-seated landslides. We show that by measuring the temperature inside the shear band of the active landslide, we are able to quantify and reduce the uncertainty of the model’s parameters, to adequately monitor and forecast the response of the selected deep-seated landslide. In particular, it is shown that by calibrating the model’s parameters through the experimental and field data of a traditionally unobservable quantity, basal temperature, the model is able to reproduce the evolution of the observable displacement of the slide.


Introduction
Deep-seated landslides are typically involving slow motion of earth over heavily deforming zones of intense shear (shear bands) at their base, before collapsing catastrophically. The shear bands are usually formed by clays or clay-like materials that can be very sensitive when the material is sheared and experiencing changes in pressure and temperature (Seguí, Rattez, & Veveakis, 2020). In earlier works, the authors (Seguí, Rattez, & Veveakis, 2020;Veveakis, Vardoulakis, & Di Toro, 2007), presented a mathematical model that is able to reproduce the behavior of a deep-seated landslide from its creeping phase (secondary creep (Intrieri, Carlà, & Gigli, 2019)) to its catastrophic collapse (tertiary creep), considering the changes in temperature of the shear band material because of frictional heating (Goren & Aharonov, 2007;Goren, Aharonov, & Anders, 2010). In this mathematical model, the constitutive equations used were following the work of Vardoulakis (2002), considering that the clay material located inside the shear band behaves as rate hardening and thermal softening, when the material is at critical state (neglecting any volumetric changes when the material deforms). However, these theoretical considerations have never been tested in a controlled case.
In the present study, we aim at validating the constitutive equations of the model mentioned previously, presented by Veveakis et al. (2007) and Seguí, Rattez, and Veveakis (2020), and the theoretical considerations of the material inside the shear band made by Vardoulakis (2002). To this end, we have instrumented a deep-seated landslide located in Andorra ( Figure 1A), called El Forn landslide (Figure1B), installing a thermometer inside its shear band ( Figure 1C). In addition, we have tested the core samples of the shear band material, to assess its rate and thermal sensitivity, and we are hereby discussing the results.
The El Forn landslide is a large deep-seated landslide located in Andorra, in SW Europe ( Figure 1A). This landslide has approximately a sliding mass of 300M m 3 (Figure 1B)that currently creeps with an average velocity of 2cm/year (Seguí, Tauler, Planas, Moya, & Veveakis, 2020). Inside the large landslide, there is an lobe that slides faster (with a maximum velocity of 4cm/year) which is the Cal Ponet-Cal Borronet lobe, with approximately 1M m 3 of sliding mass ( Figure 1B). This lobe moves as a deep-seated landslide itself ( Figure 1C inlet) and is the one that the present study focuses on. The shear band of this lobe is located at 29m depth ( Figure 1C) and is formed by Silurian shales very rich in phyllosilicates (muscovite, paragonite and chlorite) containing about 80% of the shales, and about 20% of quartz (Seguí, Tauler, et al. (2020)).
The instrumentation installed in this lobe is inside the borehole S10 ( Figure 1A), and consists of an extensometer (measuring the horizontal displacement), three piezometers (measuring the water pressure), and a thermometer inside the shear band (measuring the changes in temperature of the material, seen in Figure 1C). We observe from the displacement profile in the inlet of Figure 1C, that the landslide is indeed creeping as -2- a a rigid block on top of a deforming shear band, located at approximately 29m depth. The field data that we obtain therefore include the displacement of the slide ( Figure 1D inlet), the groundwater pressure ( Figure 1D), and the temperature of the shear band (Figure 1D). Regarding the groundwater pressure, only the piezometer that is located below the shear band ( Figure 1C) has recorded variations of the water pressure ( Figure 1D), because the landslide movement is triggered by the lower aquifer due to seasonal snow melting from the top of the mountain, and seasonal rainfall. The piezometer is 6m below the shear band ( Figure 1C inlet) and has high water pressure readings (in the order of 90-100kPa) ( Figure 1D), thus, applying high pressures to the shear band. The temperature variations that have been recorded are between 6.34 and 6.39 • C ( Figure 1D), which can traditionally be considered negligible. However, as shown in Figure 1D, the temperature signal is very sensitive to changes in pressure, essentially echoing any pressure variations within a day.
2 Rate and Thermal Sensitivity of the shear band material Apart from the monitoring data, the models are requiring knowledge of the material parameters of the shear band, such as thermal diffusivity, rate sensitivity, thermal sensitivity, as well as the thickness of the shear band. To obtain the values of the thermal and rate sensitivities, we have performed tests on the core samples of the shear band material (Seguí, Tauler, et al., 2020). The experimental tests have been performed in a triaxial machine with temperature control, on remolded core samples from the shear band of the Cal Ponet-Cal Borronet lobe of the El Forn landslide, located between 29 and 29.5m depth. These samples are Silurian shales and have been previously characterized mineralogically (Seguí, Tauler, et al. (2020)), showing that the fabric of the samples of the shear band is completely aligned and parallel to the shearing direction.
For the tests performed on all the samples in the triaxial machine, we have followed always the same steps in the same order to have consistency in the results and repeatability. As the cylindrical samples (38mm diameter, 65mm height) are not consolidated after being remolded, we have first compressed isotropically the samples at 200kPa in undrained conditions, and calculated a static friction angle of 30 • . Once the pressure of the triaxial has stabilized, and maintaining a constant isotropic compression at 200kPa, we have increased the axial load at different rates, with loading-unloading cycles until the sample reaches a critical state, at which the deviatoric (differential) stress (q), the confining stress, the pore pressure, and the temperature remain constant (Figure 2 A).
While at critical state, velocity stepping is performed (Figure 2 B) at 5 different velocities (from 0.0001−1mm/min), allowing the sample to reach a new critical state before performing the next velocity step. Through this exercise, the rate sensitivity of the material's shearing resistance at critical state, q cs , is evaluated to be ( Figure 2B inlet): where V is the velocity [mm/min], N is the rate sensitivity of the material [−], q ref and V 0 are the reference values for the shear resistance (deviatoric stress at critical state [kP a]) and the velocity [mm/min], respectively. As shown in the inlet of Figure 2B, for this ma- After the velocity steps, the axial load was kept constant at a compression of 200kPa, keeping the sample at critical-state (Figure 2A). At this point, the thermal tests start by increasing the temperature of the sample to obtain the thermal sensitivity of the material ( Figure 2C). For the thermal tests, we have increased approximately 3 degrees every hour, so that the temperature of the sample has time to stabilize and reach successive critical states. Once the temperature of the sample is constant at each temperature variation, we mark the deviatoric stress values to obtain the thermal sensitivity of the -4-manuscript submitted to EarthArXiv Following the rate and thermal sensitivities tests, we may now combine the two effects on the shearing resistance of the material (Equations 1 and 2), by accepting the multiplicative decomposition suggested by Vardoulakis (2002): whereγ ≈ V /H is the deviatoric (differential) strain rate calculated in the lab under negligible radial deformation rate (H being the height of the sample, in this case 65mm). From the experimental results of the El Forn shear band material, we obtain M = 0.04 Equation 3 can be solved for the strain rateγ, to obtain the following viscoplastic flow law:γ

Mathematical model of the shear band
The mathematical model and the constitutive equations used to forecast the behavior of a deep-seated landslide were first described by Vardoulakis (2002), and then implemented by Veveakis et al. (2007) and Seguí, Rattez, and Veveakis (2020) for the Vaiont landslide (Italy) and the Shuping landslide (in Three Gorges Dam, China). The equations used in the mathematical model focus on the behavior of the material located inside the shear band, and assume that the material is at critical state (deforming under constant volume), fully saturated in water, visco-plastic, and its mechanical properties vary along the vertical axis, z, of the shear band (inlet of Figure 1C).
Using the arguments presented in details by Rice (2006); Seguí, Rattez, and Veveakis (2020); Veveakis et al. (2007), the momentum balance equations inside the shear band ( ∂σ xz ∂z = ∂σ zz ∂z = 0) yield constant profiles of the effective stresses inside the shear band and equal to their external values: σ xz = τ d (t) for the shear stress and σ zz = σ n (t) for the normal stress. Correspondingly, since the material is at critical state (i.e. deforming under constant volume), the mass balance yields the incompressibility condition for zero volumetric strain rate˙ V =˙ zz = 0. Therefore, the main equation describing the response of the basal material is the energy equation (Rice, 2006;Seguí, Rattez, & Veveakis, 2020;Vardoulakis, 2002;Veveakis et al., 2007)), reading: with boundary conditions T = T boundary at the boundaries of the shear band, z = − ds 2 , ds 2 (ds is the thickness of the shear band). In this equation, ρC m is the heat capacity of the shear band material and c th = jk m /ρC m is the thermal diffusivity, and jk m being the thermal conductivity.
A necessary step in using the laboratory derived law of Equation 3 in the mathematical model, is to convert the invariant expression of the deviatoric stress, q, into the corresponding loading stresses in the field. This requires knowledge of the stress conditions in the field, which is not the case in the El Forn landslide. However, given that the slope of the slide (shown in Figure 1C) is approximately 30 • , which is the same as the static friction angle of the material, we may use Rankine's theory of slopes (see Seguí, Rattez, and Veveakis (2020) and Collin (2001)) to deduce that the normal stresses acting on the shear band should be equal (σ zz = σ xx = σ yy ), and the only shear stress applied on the shear band should be in the down-slope direction, σ xz . Based on that, we can use the definition of the deviatoric stress in the conditions considered, q cs = √ 3J 2 = 3σ 2 xz , and the definition of the deviatoric strain rate at critical state,γ =˙ xz √ 3. Through these expressions, the dissipation (source) term of Equation 5 becomes, in light of Equation 4: To fully characterize the dependence of the shear stress, τ d , on the groundwater pressure and its variations, a regional hydro-geomechanical model is required, as described by Seguí, Rattez, and Veveakis (2020). In the specific case under consideration, the El Forn landslide, such an analysis cannot be easily performed since the landslide is fed/loaded from the pressure changes of the groundwater below the shear band. These changes are not communicated to the overburden, since the shear band is acting as a flow barrier. This, in turn, may suggest that the overburden is subjected to a constant gravitational load, and water pressure variations below the shear band directly affect the shearing resistance of the landslide, with the rim of the shear band always remaining at critical state.
To assess the validity of this suggestion, we will hereby use the outcomes of the hydrogeomechanical analysis of Seguí, Rattez, and Veveakis (2020), suggesting that the shear stress, τ d , varies linearly with the water pressure, p f , recorded below the shear band. We will therefore set τ value of the shear stress applied on the shear band, when the fluid pressure is at the reference value, p f 0 . The validity of this assumption will be re-assessed and discussed after the modeling and forecasting process, presented in the next sections.
Following all these considerations, Equation 5 can be combined with Equation 6 and Equation 4, and further be reduced to a single parameter dimensionless equation where the following dimensionless parameters have been used: The dimensionless group, Gr, is the so-called Gruntfest number (Gruntfest (1963)), defined as follows: with The Gruntfest number, Gr, expresses the ratio of the mechanical work converted into heat over the heat diffusion capabilities of the material. This parameter includes all the material properties at hand (thermal conductivity, rate and thermal sensitivities, and reference rate), as well as the thickness of the shear band and the shear stress, τ d , applied on the shear band from the external loading sources (gravity and groundwater). Following the considerations thus far, Gr evolves with the groundwater pressure (p f ), and therefore, with time, acts as a link between the external loading conditions with the internal response of the material. Since this is the only parameter of the mathematical equations, the performance of the model depends on accurately constraining the value of Gr. This will be done next, by combining the experimental and field data, together with probabilistic considerations for the least constrained parameters in Gr.

Modeling and forecasting the landslide's response
In this section, we will use the experimental results presented in Figure 2, together with the field data of Figure 1D, to calibrate the parameters of the mathematical model of Equations 4 and 7. The process that will be followed includes the subsequent steps: 1. Input parameters: The first step is to input as many parameters as possible, so that the Gr number of Equation 7 can be constrained in terms of both possible range of values, as well as their probability distribution. 2. Uncertainty minimization through the field temperature: Once this step is performed, the temperature equation, Equation 7, will be solved using a Monte-Carlo scheme (Raychaudhuri (2008) and Hastings (1970)), and the field temperature will be used to obtain the optimal value of Gr, thereby minimizing its uncertainty. 3. Forecasting and parameter constraint: With a calibrated temperature equation, the next step is to forecast the field displacement shown in the inlet of Figure 1D. This will be done using Equation 4, adjusted for the problem at hand.
These steps are detailed in the following paragraphs.

Input parameters
As it becomes obvious by the definition of the only dimensionless group of the problem, the Gr of Equation 9, incorporates all the material, environmental and loading parameters of the problem at hand, and needs to be constrained so that the model can be tested. Therefore, as a first step, we rely on the results obtained from the experimental tests, providing values for the rate sensitivity, N = 0.0136 and the thermal sensitivity, M = 0.04 • C −1 of the shear band material. In addition, we obtained the reference laboratory values for q ref = 719.4 kP a andγ ref = 1.5 × 10 −2 s −1 . The former two parameters have been evaluated at laboratory conditions where T lab = 20 • C, while the field features at an ambient temperature of T boundary = 6 • C (see Figure 1D). Finally, from the field data of the groundwater pressure, shown as the blue line in Figure  1D, we identify as reference value for the pore fluid pressure the initial value of the recordings, p f 0 = 90 kP a.
In addition to the above, since the shear band material consists of phylosilicates (Seguí, Tauler, et al., 2020), we accept standard literature values (Vardoulakis, 2002) for their thermal conductivity, jk m = 0.45 J/ • Cms and thermal diffusivity, c th = 1.6 × 10 −7 m 2 /s. Note that these values have small deviations across all clay materials (Ghuman & Lal, 1985) and are hereby considered to have no uncertainty.
Having constrained these values, the remaining parameters in the expression of Gr are the reference values of the loading stress τ d,ref , and the shear band thickness, d s . Both quantities cannot be easily determined from field data of a single borehole, mainly because of the lack of representativeness of the values. As discussed by Seguí, Tauler, et al. (2020), the material retrieved from the borehole S10 of the field (see Figure 1B) suggests that the shear band thickness, d s , should be between a few centimeters, up to half a meter. However, this is information received from a single borehole, and cannot be constrained or validated further. Correspondingly, the loading stress, τ d,ref , at the beginning of the recording cannot be easily determined, since the landslide is mainly driven by the over-pressurized aquifer below the shear band. To constrain the loading stress one should develop a regional hydro-geological model of the water flow, to calculate the volume of the aquifer, the induced seepage force, and its inclination with respect to the slope of the landslide.
In the absence of further information for both, the shear band thickness (d s ) and the loading stress (τ d,ref ), ranges of values were given to these parameters (d s = 1 ÷ 50 cm and τ d,ref = 100 ÷ 1000 kP a). We then assessed, based on the mineralogical study of Seguí, Tauler, et al. (2020) on the shear band material, that d s should obey a normal distribution with its mean (most probable) value being d s = 10 cm, and τ d,ref should obey a uniform distribution between its end member values. Under these considerations, the base value of the Gruntfest number, G 0 (Equation 10), was found to obey a log-normal distribution.

Uncertainty minimization through the field temperature
The next step in the modeling process is to solve the stochastic equation for the basal temperature, Equation 7. To do so we have used a Monte-Carlo sampling approach (Hastings, 1970;Raychaudhuri, 2008), by sampling 500 random points in the log-normal distribution of G 0 , shown in the inlet of Figure 3A.
Following this sampling, Equation 7 was solved 500 times, and the solutions have been compared with the field temperature, as shown in Figure 3A. By calculating the least square error of the Monte-Carlo simulations with the field temperature, in the training window highlighted in Figure 3A, we then obtained the simulation with the least error. This solution allowed us to obtain the value of G 0 that minimizes the error and fits the field temperature best, G 0 = 0.71 [−] ( Figure 3A). Through this value, and the definition of G 0 (Equation 10), we obtain a first scaling equation between the loading stress and the shear band thickness ( Figure 3B):

Forecasting and parameter constraints
Once the temperature is calculated, the next step is to determine the velocity, V , and cumulative displacement, u, of the landslide. This is achieved in our model by integrating in time and space the strain rate (Equation 4). For the velocity: where V 0 is the initial velocity of the field, which should also satisfy the following equation of the model (from Equation 4): The initial velocity of the field can be calculated from the displacement data to be approximately V 0 ∼ 2 − 4 cm/year, a value that is in accordance to literature values as well (Corominas et al., 2014). Using the parameter values previously discussed, Equation (13) can be solved for τ d,ref to yield: By requiring Equations 11 and 14 to be equal, we may solve for V 0 (in cm/yr) as a function of d s (in cm), to obtain:   By assuming one of the above inversion set of parameters, for example the average value of V 0 = 3 cm/yr (hence setting d s = 25 cm and τ d,ref = 540 kP a in the remaining parameters), the system's parameters are fully constrained, and the mathematical model should be able to forecast the displacement of the landslide. Indeed, upon numerical time integration of the velocity of the model, Equation 12, the model forecasts satisfactorily the landslide's displacement, as shown in Figure 3C.

Discussion and Conclusions
Following this parameter constrain and forecasting exercise, what remains to be discussed is the validity of the parameter inversion and the assumptions used. For the former discussion, the model has received input from experimental and field data and provided output for the shear band thickness and the loading stress of the landslide. It was deduced that, for the background velocity range derived from the data, but also re-ported in the literature (Corominas et al., 2014), the shear band thickness was inverted to be between 15−35 cm and the loading stress between 530−550 kP a. This is a reasonable range for the shear band thickness, which can be confirmed through the observations from the core logs, and in accordance with shear band thicknesses inverted for other deep-seated landslides by the authors. Indeed, Veveakis et al. (2007) calculated a shear band thickness of approximately 16 cm for the Vaiont landslide in Italy, and Seguí, Rattez, and Veveakis (2020) derived a shear band thickness of approximately 70 cm for the Shuping landslide in China. The median value of 25 cm for the El Forn landslide is therefore in line with these results.
Regarding the inverted range of values for the loading shear stress below the shear band, τ d,ref = 530−550 kP a, its validation requires in principle detailed hydro-geological data like the seepage flow-lines under the shear band and the size of the over-presurized aquifer. In the absence of such information, we may attempt to assess the validity of this range of values, by recalling that near impermeable barriers like the shear band, the flow is parallel to the barrier. Since the aquifer is over-pressurized, the driving stress would consist of a gravitational component and a seepage component (see the analysis of Seguí, Rattez, and Veveakis (2020)): τ d,ref = ρgHsin(δ) + ρ w gh aq sin((δ + α)/2). In this expression ρ = 2.5 g/cm 3 is the density of the overburden at depth H = 29 m, ρ w = 1 g/cm 3 the density of water, g = 9.81 m 2 /s the acceleration of gravity, δ = 30 • the slope of the landslide, α = δ the slope of the seepage force and h aq the height (depth) of the parallel flow inside the over-pressurized aquifer below the shear band. For these parameter values, the range of τ d,ref = 530−550 kP a is achieved for h aq = 30−35 m, which is a reasonable size for parallel flow in an over-pressurized aquifer in a mountainous area.
After assessing the validity of the inverted values for the parameters, the next step is to discuss the validity of the assumptions of this study. The main assumptions were that: the thermal sensitivity of the material is sufficient to drive the creeping process; the loading stress is scaling linearly with the groundwater pressure. The latter assumption was discussed previously in the context of the seepage forces exerted on the shear band by the aquifer. The former assumption, which was this work's main hypothesis has been illustrated through a combination of field data, experimental results and modeling approaches. Through this combination of techniques, the present study has shown the importance of the temperature in the shear band of a landslide. The case study that has been used is the Cal Ponet-Cal Borronet lobe, in the El Forn landslide (Andorra). This lobe creeps as a deep-seated landslide, and has been instrumented with a thermometer in the shear band, piezometers, and an extensometer.
The field data of this lobe shows that the movement of the slide is triggered, first, by the increase of water pressure in the shear band (by the lower aquifer), and, second, by an increase of the temperature in the shear band, resulting in the acceleration of the sliding mass. In this letter, we have presented the experimental results of the rate and thermal sensitivities of the shear band's material of the lobe. The results of the laboratory tests have shown that the material inside the shear band is indeed thermal and rate sensitive. The mathematical model implemented and the constitutive equations assumed and redefined by the laboratory experiments, have demonstrated the causes of the acceleration of the lobe, as well as forecasted the behavior of the landslide. The field data of the basal temperature and the laboratory experiments have reduced the uncertainty of the model, allowing us to decrease the parameters needed to be inverted, and forecast the behavior of the landslide. The present model can, therefore, become a useful tool to forecast the landslide motion through field and laboratory data.