*2.1. Site Characterization*

The Dollendorfer Hardt is a horst structure with a height of roughly 250 m above sea level, representing the most northern peak of the Siebengebirge of volcanic origin. Towards the north, there is the lower Rhine Bay. The lithology of the area is characterized by Devonian shales, consisting of interchanging bedding of slate and graywacke, on which Tertiary sediments, mainly clay and quartz sand, are deposited. These lithological layers are partly overlain by basaltic and trachytic tuffs as a result of volcanic eruptions in this region [39].

In order to obtain subsurface information for the study site, drillings were performed between March and August of 1998 [38]. The techniques included closed core percussion drilling, open core percussion drilling and manual auger and percussion drilling (Figure 2). The soil of the drill cores was used to determine hydraulic and mechanical properties of the different identified layers according to the guidelines of the German Industrial Standard of that time (DIN 1993) [38]. In particular, the following parameters were determined: particle size distribution, soil water content, consistency limits, particle density and maximum soil water content but were not completed for all soil types due to the limited availability of samples. Further, shear strength was determined by direct shear tests and triaxial tests. Samples from outcrops were used to determine the hydraulic conductivity [38].

**Figure 2.** (**a**) Digital elevation model of the study area provided by the Geobasisdatendienst NRW with geoelectric monitoring profile line, seismic profile line, borehole locations and water content sensors, (**b**) as well as a simplified geomorphological map.

The scarps of the two landslides at the study site cover an area of almost 30,000 m<sup>2</sup> located at a steep (up to 35–40◦) south-facing hillslope with a river at the valley bottom. The area is approximately 330 m long with a variable width ranging from 17 m at the narrowest passage, up to 65 m close to the scarp (Figure 2). The well defined scarp is the northern boundary of the landslide zone and is a result of a rotational movement. Relatively undisturbed rotational blocks can be identified in the landslide mass in the upper part of the slope [38], which then convert to a mass flow in the middle part of the slope (transport zone). The transport zone is the narrowest part of the landslide and over 140 m long. The earthflow is constrained by one, sometimes two sets of levees at both sides, originating from the two landslide events. The translational mass has been deposited in the debris zone in the lower part of the slope up to a small river with little to no inclination [39,42]. The landslide mass consists mostly of trachytic tuff and clayey sediments from the Tertiary. The first landslide mass additionally contains loess from the Quaternary loess cover. The slope is naturally covered by a forest consisting mainly of beech trees.

Previous studies provided evidence of regular slope movement in the middle part of the scarp area using inclinometers and tiltmeters. Movement of ±3 cm year<sup>−</sup><sup>1</sup> was observed for the transport zone in 'extraordinary wet conditions' in spring when heavy rainfall coincided with high groundwater levels [38]. In addition, continuous soil creep of small magnitude on the order of a few millimeters per year was observed in the transport zone [40]. Schmidt [38] attributed the vulnerability of the site for slope instability to the specific geological setting of the area with an abundance of clay-rich soil layers. This also explains the reported elastic swelling/shrinking of the lower rotational block of the scarp zone and the elastic movement associated with groundwater level changes in the debris zone [38]. Ruptures and breakups of tuff blocks of ≈0.3 m can be observed regularly along the scarp.

To investigate the extent of the lithological layers in the study area, ten seismic refraction profiles were acquired between September 2015 and January 2017. The measurement profiles had lengths between 50 to 150 m with a geophone spacing varying between 1 and 2.5 m for the different profiles. The refracted seismic waves were generated by strikes with

a 5 kg hammer on a metal plate close to each geophone. For a good signal to noise ratio, ten strikes per geophone were stacked. Signals were recorded with a SUMMIT II Compact (DMT-Group, Essen, Germany). The maximum recording time was 350 ms. The data were processed and analyzed using the software ReflexW (J. Sandmeier, Karlsruhe, Germany). A bandpass filter was applied to cut frequencies below 10 Hz and above 150 Hz. First arrivals were picked and an inversion was performed using regular grids with a grid spacing of one quarter or less of the geophone distance.

The results for the first arrival of the refracted seismic waves indicated a three layer case (Figure 3). The inversion results sugges<sup>t</sup> that the seismic wave velocity range associated with the three layers are <300 m s<sup>−</sup>1, 400–600 m s<sup>−</sup>1, and >800 m s<sup>−</sup>1. These layer velocities were consistently found in all measured profiles. Based on the seismic data and additional core drilling data from [38] and literature values for the rock and mineral types [43,44], the three different layers are interpreted as follows:


The seismic measurements also indicated isolated sand structures within the Tertiary sediments. Furthermore, the thickness of the Tertiary sediments was found to decrease downhill and to disappear completely within the transport zone. The abrasive character of the former landslides eroded the originally deposited sediments and replaced them with the deposited landslide mass. In the transport zone, the sliding surface corresponds to a lithological boundary and was identified in the seismic refraction. In the upper part of the rotational landslide, the sliding surface lies within the Tertiary sediments and therefore cannot be determined by seismic refraction [38,39]. Based on previous studies, the upper layer is assumed to be temporally unstable and thus prone to landslides [41]. A similar layering was found on the adjacent slopes. This supports the assumption that only the upper layer of tuff and parts of the Tertiary sediments were transported and mixed by the landslides.

**Figure 3.** Tomogram of p-wave velocity derived for seismic profile 1 (see Figure 2a) and its subsequent interpretation.

## *2.2. Meteorological Data*

Precipitation was continuously measured at the meteorological station at Königswinter-Heiderhof as well as at the station of the Department of Meteorology of the University of Bonn (MIUB). The Königswinter-Heiderhof station is located 3 km southwest of the test site, whereas the MIUB station is approximately 10 km northwest of the study site (Figure 1). At the Königswinter-Heiderhof station, precipitation and temperature data are recorded at a daily resolution since 1990 and 2000 onward, respectively. For the MIUB station, daily precipitation data are available between 1999 and 2001, and 10-min resolution precipitation data are available since 2010. Temperature is recorded with a 10 min resolution since 1995. The mean annual precipitation for the time period of 1995–2017 varied between 650–950 mm with an average of 769 mm. For the same period, the mean annual temperature varied between 8 ◦C and 13 ◦C with an average of 9 ◦C (Figure 4). Since temperature data were missing at the Königswinter-Heiderhof station for the period 1995– 1999, temperature data from the MIUB station were used in this period. As common in the temperate climate zone, the strongest rainfall events were associated with thunderstorms in summer. The maximum monthly precipitation was 235 mm in July 2014.

**Figure 4.** (**a**) Monthly precipitation, (**b**) temperature, (**c**) potential evapotranspiration, and (**d**) net infiltration for the Dollendorfer Hardt test site based on the data of Königswinter-Heiderhof station from 1995–2017.

The net infiltration, *Inet* [mm month−1], is the difference between precipitation *P* [mm month−1] and actual evapotranspiration *ET* [mm month−1]:

$$I\_{\rm nt} = P - ET \tag{1}$$

To obtain net infiltration from precipitation, an estimate of the actual evapotranspiration is required. To minimize the number of required parameters, the temperature-based approach of Thornthwaite [45] was used here. In this approach, the monthly potential evapotranspiration was calculated using

$$PET = 16 \left( \frac{L\_d}{12} \right) \left( \frac{N}{30} \right) \left( 10 \frac{T\_d}{l} \right)^a \tag{2}$$

$$a = (6.75 \times 10^{-7})I^3 - (7.71 \times 10^{-5})I^2 + (1.792 \times 10^{-2})I + 0.49239 \tag{3}$$

$$I = \sum\_{i=1}^{12} \left(\frac{T\_{ai}}{5}\right)^{1.514} \tag{4}$$

where *PET* [mm month−1] is the estimated potential evapotranspiration, *Ld* [h] is the average daylength of the month, *N* [-] is the number of days of the month, *Ta* [◦C] is the mean daily temperature of the month and *a* [-] is an exponent that is a function of the annual Thornthwaite heat index, *I* [-]. Based on this approach, the mean annual *PET* was estimated to be 644 mm for the Dollendorfer Hardt. This value is in good agreemen<sup>t</sup> with *PET* values provided for the Roleber station (8 km north of the study site) by the Deutscher Wetterdienst (DWD, German Weather Service), which were calculated based on the Haude method [46].

The resulting monthly *PET* for the study area is shown in Figure 4. It can be seen that monthly *PET* occasionally exceeded precipitation in the summer months when temperature is relatively high. In these cases, net infiltration was considered to be negative, which was modeled as outflow. The mean annual precipitation, *PET*, and net infiltration for the simulation period of 1995–2017 were 769, 644, and 125 mm, respectively. It should be noted that it was assumed here that PET was equal to actual evapotranspiration. This assumption is certainly not valid in all conditions [47]. However, the Dollendorfer Hardt is experiencing precipitation all over the year and has a relatively high groundwater level. Therefore, it was assumed that the trees did not experience water stress and were able to transpire with rates dictated by the atmospheric conditions.

Slope instability at the Dollendorfer Hardt site is suspected to be triggered by intensive rainfalls. In order to identify the potential magnitude of high-intensity rainfall events both the daily rainfall data from the Königswinter station and the high-resolution data from the MIUB station were analyzed. The maximum daily rainfall amount at the Königswinter and MIUB stations are 88 mm and 45 mm, respectively. The maximum hourly precipitation per year was derived from the MIUB station with 35 mm h−<sup>1</sup> in 2013. The mean value of maximum hourly precipitation over the years 2010 to 2017 was 15.7 mm h−1. Based on this analysis, a rainfall intensity of 20 mm h−<sup>1</sup> was selected to represent an intensive rainfall event at the study site which is also not too rare.

#### *2.3. Monitoring of Ground Water Level and Soil Water Content*

Following the drilling in 1998, 26 standpipe groundwater gauges were installed [38]. For long-term groundwater monitoring, an electronic water level indicator was used. Twelve tubes showed strong variations in groundwater level and were monitored with hourly resolution with pressure transducers [38]. Groundwater monitoring data are available from April 1999 to May 2001 [38].

To characterize spatial variability of near-surface soil water content, the soil water content sensor network SoilNet [32] was installed at the Dollendorfer Hardt test site in August 2016. The installed network consisted of 20 SoilNet nodes. The locations of the nodes were chosen to achieve a homogeneous distribution of sensors across the slope. At each location, six SMT100 sensors (Truebner GmbH, Neustadt, Germany) were installed horizontally at depths of −0.05, −0.20 and −0.50 m. Two sensors at each depth were used to increase the measurement volume and to provide a control of the measurement quality. The operating principle of the SMT100 sensor and the calibration approach for

determining the relative dielectric permittivity from the sensor response are described in Bogena et al. [48]. To link the measured soil permittivity to soil water content using petrophysical relationships, 12 undisturbed soil samples were taken at 4 locations along the slope at the depths of −0.10, −0.30 and −0.50 m. The soil samples were saturated in the laboratory and the soil permittivity as well as the weight was measured daily during a drying period of 42 days at room temperature. Subtracting the dry weight of the samples from the measured weight, the water content was determined and linked to the soil permittivity. The samples were roughly categorized into two classes based on their grain size distribution. The coarse-grained samples were located in the rotated blocks while the fine-grained material was found at all other locations. No difference was found for the petrophysical relationship for clay and tuff samples. For the coarse-grained soil samples, the petrophysical relationship of Roth et al. [49] showed the best agreemen<sup>t</sup> and was used for the conversion of the soil permittivity to water content for the SoilNet sensors 11–14 and 16–18. For the fine-grained samples, the petrophysical relationship of Robinson et al. [50] showed the best agreemen<sup>t</sup> with the laboratory data and was subsequently used for the remaining SoilNet sensors. Soil water content measurements were taken from August 2016 to July 2018.

To assess the water dynamics at greater depth and to achieve a higher spatial resolution than the SoilNet sensor network, an electric resistivity monitoring system was installed in the transport zone of the landslide. Measurements were conducted between March 2016 and May 2018 with various time intervals from daily to monthly measurements. An ABEM Terrameter LS (Guideline Geo AB, Sundbyberg, Sweden) with 96 electrodes was used with an electrode spacing of 0.5 m. A dipole–dipole configuration with skips of 0, 2, 4 and 6 electrodes was combined with multiple gradient measurements. In the data processing prior to the inversion, data points were removed due to systematic errors, such as bad electrode connections, problems with power supply or high current strength (>1 A). To account for temperature effects, the electric conductivities were corrected to the mean subsurface temperature of 10 ◦C following the procedure of Hayley et al. [51]. Due to the lack of temperature measurements at depth, the model of Brunet et al. [52] was used to calculate the required temperature information *<sup>T</sup>*(*<sup>t</sup>*, *z*)[ ◦C] over the given time and depth range. The preprocessed electrical data were inverted using the finite element based inversion code CRTomo [53]. We used a resistance error model with parameters *a* and *b* for absolute and relative resistance error contributions, resulting in a resistance error Δ *R* [Ω] for the measured resistance *R* [Ω] [54]. The absolute error was set to *a* = 0.001 Ω and the relative one to *b* = 3%. To improve the resistivity estimation, the inversion was performed with the seismic layer boundary between bedrock and landslide as a structural constraint e.g., [18,55,56]. To better uncover temporal changes, a difference inversion was also performed [57].

In the ERT inversion results, the three layers found with seismic refraction could not be identified, as the difference in electric conductivity of the two upper layers was too small. The inversion results only show the electrically conductive and heterogeneous landslide mass and the more resistive bedrock. From these results, water content was estimated using the relationship as described in Waxman and Smits [26]:

$$
\sigma\_b = \frac{S^l}{F\_f} \left( \sigma\_w + \frac{\sigma\_s}{S} \right) \tag{5}
$$

with saturation, *S* [-], water conductivity, *σw* [S m<sup>−</sup>1], surface conductivity, *σs* [S m<sup>−</sup>1], formation factor, *Ff* [-] and saturation exponent, *i* [-]. Water conductivity, *σw* = 0.1 [S m<sup>−</sup>1], was gained from in situ measurements in the boreholes. As surface conductivity depends on the clay content, the empirical relationship between clay content, *cc* [%] and the surface conductivity, *σs* [S m<sup>−</sup>1], in the work of Rhoades et al. [58] was used:

$$
\sigma\_s = (2.3 \times \text{cc} - 2.1) \times 10^{-3} \tag{6}
$$

where *cc* is the clay content (*cc*1 = 57% and *cc*2 = 40%, where the index 1 denotes the upper layer and 2 the lower layer) obtained from Schmidt [38]. The formation factor *Ff* = *ϕ*<sup>−</sup>*<sup>j</sup>* was calculated based on the porosity *ϕ* [-] and the cementation exponent *j* = 2, where the porosity values (*ϕ*1 = 55%, *ϕ*2 = 40%) were also taken from Schmidt [38] and the cementation exponent as well as the saturation exponent *i* = 2 were based on literature values e.g., [59].

The derived water content values from ERT (Figure 5) were in agreemen<sup>t</sup> with values measured by the SoilNet sensor network, which were limited to a depth of −0.50 m of the soil. In general, observed water dynamics were low during the monitoring period, as only the first two meters showed a response to precipitation events and a decrease in soil water content during dry periods. However, the change in volumetric water content below the top layer was less than 2 cm<sup>3</sup> cm<sup>−</sup><sup>3</sup> and declined strongly with depth over a 7 week dry period from 5 April to 24 May 2017. Continuous low-intensity precipitation with less than 10 mm d−<sup>1</sup> over several days resulted in increasing saturation in the upper 0.50 m of the soil.

**Figure 5.** Water content derived from ERT measurements taken at 31 May 2017.

## **3. Hydromechanical Model**

The applied hydromechanical model is based on the approach of Lu et al. [5] solving the Richard's equation to describe the transient water flow in the subsurface. In the mixed form, the Richard's equation is given as

$$
\nabla K(h)\nabla H + \mathcal{W} = \frac{\partial \theta(h)}{\partial t},
\tag{7}
$$

with the hydraulic head, *H* [m], pressure head, *h* [m], hydraulic conductivity, *K*(*h*)[m s<sup>−</sup>1], volumetric water content, *θ*(*h*)[m<sup>3</sup> m<sup>−</sup>3] and possible source/sink terms, *W* [s−1]. Following the model of van Genuchten [60], the water retention curve is given by

$$S\_{\varepsilon} = \frac{\theta(h) - \theta\_r}{\theta\_s - \theta\_r} = \left(1 + (a|h|)^n\right)^{-m} \tag{8}$$

with effective saturation, *Se* [-], residual and saturated water content, *θr*/*s* [m<sup>3</sup> m<sup>−</sup>3], and soil specific parameters, *α* [m−1], *n* and *m* = 1 − 1/*n* [-]. The hydraulic conductivity is based on the effective saturation

$$K(h) = K\_s \sqrt{S\_c} \left( 1 - (1 - S\_c^{1/m})^m \right)^2 \tag{9}$$

with saturated hydraulic conductivity, *Ks* [m s<sup>−</sup>1].

The mechanical response of the soil is considered to be linear elastic with no poroelasticity and described by the Young's modulus, *E* [Pa] and Poisson's ratio, *ν* [-]. In this case, momentum equilibrium is given by the total stress tensor, *σ* [Pa], the unit weight, *γ* [N m<sup>−</sup>3] and the body force vector, *Fv*[-], as

$$
\nabla \mathcal{L} \sigma + \gamma \mathcal{F}\_{\overline{\nu}} = 0 \tag{10}
$$

From the total stress tensor, the effective stress tensor, *σ-*[Pa], can be calculated based on the pore water pressure, *Pw* [Pa], [61,62]:

$$
\sigma' = \sigma - S\_\text{t} P\_\text{w} I \tag{11}
$$

with the unit vector, *I*[-]. To assess the stability of a hillslope, the local factor of safety, *LFS* [-], is calculated based on the Mohr–Coulomb failure criteria

$$LFS = \frac{\cos \phi'(c' + \sigma\_m \tan \phi')}{R} \tag{12}$$

with effective friction angle, *φ* [◦], effective cohesion, *c* [Pa], mean effective stress, *σm* = 0.5 × (*<sup>σ</sup>*1 + *<sup>σ</sup>*3) − *SePw* [Pa] and the radius of the Mohr circle, *R* = 0.5 × (*<sup>σ</sup>*1 − *<sup>σ</sup>*3)[Pa]. Here, *<sup>σ</sup>*1,3 [Pa] are the major and minor principal stresses. A slope is considered to be prone to failure for *LFS* < 1 as the restricting forces are smaller than the downhill forces in this case.

#### **4. Combination of Field Observations and Hydromechanical Modeling**

In this study, the hydromechanical model was used to assess slope stability through the calculation of the local factor of safety. Site characterization, laboratory tests and continuous monitoring of soil water content and precipitation are input for the model to increase the quality of the assessment (Figure 6). The input can be separated into one-time constraints and dynamic information. The base for the numerical model is the physical/mathematical model as described in Section 3. To solve the governing equations, the parameters, the geometry of the domain as well as boundary and initial conditions need to be specified.

**Figure 6.** Schematic sketch of the different input streams into the hydromechanical model for stability analysis.

## *4.1. Model Setup*

Surface topography was taken from a digital terrain model with a spatial resolution of 1 m<sup>−</sup><sup>2</sup> derived from remote sensing and provided by the Geodatenbasisdienst NRW. Bedrock topography as well as the varying thickness of soil layers was obtained from the seismic refraction and the drilling logs. For this purpose, the wave velocities of 300 ms<sup>−</sup>1, 600 ms<sup>−</sup><sup>1</sup> and 800 ms −1 were selected as limit values for the lithological layers in agreemen<sup>t</sup> with the lithological units identified in the borehole logs. The position of the interfaces

was digitized in the seismic tomograms, georeferenced and converted to absolute depths using the digital terrain model. The information was entered into the open source software GrassGis (GrassGis Development Team, 2017) and interpolated to surfaces using inverse distance weighting. Hence, a volume model of the investigated area was created based on the generated surfaces (Figure 7). The geophysically-derived volume model was used to generate a 3D geometrical model of the study site in Comsol Multiphysics [COMSOL Inc, Stockholm, Sweden].

**Figure 7.** Derived 3D geological model based on refraction seismic survey and borehole logs. Digital elevation model provided by Geobasisdatendienst NRW.

Despite the complex geometry of the study site that can influence subsurface flow and the resulting stability, the hydromechanical model was applied to the mid-cross-section of the landslide area only for computational reasons. The transition from the 3D to 2D geometrical model and the hydromechanical simulations were performed by COMSOL Multiphysics. The 2D domain was discretized using an unstructured triangular mesh with an increasing mesh size from surface to bottom, so that the highly dynamic hydrological conditions in the top soil are captured with a reasonable computational efficiency and accuracy. The mesh size near the slope surface was ≈0.05 m. The maximum mesh size in the top layer was ≈0.2 m. The maximum mesh size increased to 0.3 m in the mid layer and increased further towards the deeper part of the bottom layer. The modeling domain was defined to be substantially larger than the region of interest to reduce the impact of the boundary conditions.

The groundwater monitoring data from the boreholes was used to define the hydraulic boundary conditions (Figure 8). It has been observed that the groundwater level can occasionally reach the surface in the lower part of the slope. Accordingly, the slope surface is defined as a mixed boundary. The inflow and outflow rate depend on pressure, net infiltration, storage capacity of the soil and also its hydraulic conductivity. The bottom boundary was defined as a no-flow boundary. A fixed pressure head boundary at the lower right lateral boundary of the domain was defined using the river level at the slope toe. Here, it was assumed that the river level was constant and at a height of 100 m above the bottom of the modeling domain. The top 25 m of this boundary was a seepage boundary in which water is free to exit from the saturated subsurface. On the left side of the domain, a pressure head boundary was defined using the minimum level of the measured groundwater level

at the closest borehole to this boundary. The upper 38 m of the left lateral boundary was a no-flow boundary. From the mechanical point of view, the ground surface is a free boundary with no external loads and constraints, whereas the lateral boundaries were defined as Roller boundaries and the bottom boundary was fixed (Figure 8).

**Figure 8.** The 2D cross-section used for hydromechanical modeling including boundary conditions for the hydrological and mechanical model components.

In order to parameterize the three layers of the subsurface model, each layer was considered as a homogeneous and isotropic medium. Estimates for bulk, dry and saturated soil density, particle size distribution, porosity, soil cohesion, friction angle and saturated hydraulic conductivity were given by Schmidt [38] and are summarized in Tables 1 and 2. The bulk and dry density of the soil were determined by oven drying following the German Industrial Standard (DIN) of that time, DIN [63]. The maximum moisture content was determined using the method described in DIN [64]. The saturated conductivity was determined using a constant head as well as a falling head permeameter based on DIN [65]. The soil particle distribution were defined by analyzing the samples taken from the drilling cores along the landslide zone and by lithological interpretations of the borehole logs [38]. Soil cohesion and friction angle were determined by a shear box and a triaxial test based on Schmidt [38], DIN [66]. The laboratory tests showed good repeatability with a statistical variance between repeated measurements of less than 10% [38]. The Tertiary clay is the material with the highest density but the lowest cohesion (Table 2). The Devonian clay/silt has the highest cohesion of around 30 kPa but the Trachyte tuff has the highest friction angle of all tested materials. The derived maximum water content is around 35 to 40% for all materials, while the residual water content is around a 6%. The hydraulic conductivity for the Tertiary clay and the Devonian clay/silt layers were provided as a range of possible values and constrained through a calibration of the model with regard to the measured mean groundwater levels of the period 1999 to 2001 by assuming a net infiltration of 125 mm year<sup>−</sup>1. Using the available information on bulk density and particle size distribution (Table 1), estimates of the soil hydraulic parameters of each layer were obtained using the Rosetta pedotransfer function (Rosetta Lite v 1.1) [67] (Table 2). The soil elastic moduli (*E*, *ν*) were estimated from typical ranges provided in the literature [68,69] (Table 3). Here, we have considered a stiff clay-rich Trachyte tuff and Tertiary clay layers with medium to high plasticity (*E* = 15 MPa) and a harder Devonian clay/silt layer with relatively lower plasticity (*E* = 30 MPa). The typical Poisson's ratio, *ν* [-], for silty soils is given as 0.3–0.35 and for unsaturated to saturated clay as 0.1–0.5. Accordingly, for the unsaturated Trachyte tuff and Tertiary clay layers and saturated Devonian clay/silt layer, a Poisson' ratio of 0.35 is considered. All parameters required for the simulations are listed in Table 2. While all parameters may be subject to changes over time in principle, e.g., due to ongoing compaction or root growth [70], they were considered to be constant in our simulations because the test site showed little dynamics in recent monitoring periods.


**Table 1.** Particle size distribution and bulk density of each soil layer of the Dollendorfer Hardt test site (derived from Schmidt [38]).

**Table 2.** Soil model parameters for the Dollendorfer Hardt test site.


\* Calculated value by the Rosetta pedotransfer function. \*\* Values derived from literature. \*\*\* Obtained by model calibration using mean groundwater level in the period of 1999–2001 using a net infiltration of 125 mm year<sup>−</sup>1.



#### *4.2. Model Calibration and Validation*

The model was initialized by simulating the long-term average state of the slope using a spin-up period of 300 years with a constant mean annual net infiltration of 125 mm. This long spin-up period was used to ensure that the model reaches steady-state flow condition. The hydraulic conductivity of the three soil layers was calibrated with regard to the measured mean groundwater level for the period of 1999–2001 using data from three boreholes along the mid cross-section of the test site within the value range given above. The mean groundwater depth from the surface at the boreholes B8, B14 and B21 in the measurement period of 1999–2001 are −7.4 m, −2.5 m, and −2.8 m, respectively (Figure 2). For every set of hydraulic conductivity values, the 300-year model initialization was repeated. For the final choice of parameters, the measured and simulated values were highly correlated (R<sup>2</sup> > 0.9) but the simulated mean groundwater level was 0.66 m lower than the measured level. After calibration and model spin-up for 300 years, the groundwater level is lower in the upper part of the slope and closer to the surface in the middle part of the slope, which is in agreemen<sup>t</sup> with the measured values at the site.

The calculated steady-state pressure distribution was considered as the initial condition for the next simulation step. The mean monthly conditions for the test site were simulated using the mean monthly net infiltration for the period of 1995–2017. These simulation results with the calibrated model were verified using two series of measured data. First, the simulated groundwater level was compared to the mean monthly groundwater level (Figure 9) by averaging daily measured groundwater levels. Second, simulated soil water content for the top soil was compared to mean monthly measured soil water content obtained from the SoilNet data at the site. Here, the measured soil water content from the

twelve sensors nodes that are located along the mid-cross-section were used (sensor nodes 1, 2, 3, 4, 5, 8, 9, 10, 12, 14, 18 and 20 in Figure 2). It should be noted that unreliable parts of the measured data were not considered.

**Figure 9.** Time series of measured and simulated mean monthly groundwater level for the three boreholes used for model calibration and validation.

This simplified model could capture some of the key features that are deemed important for slope stability assessment. The seasonal pattern of the simulated groundwater levels correlated reasonably well with the measured values (R2-values ranging from 0.49 to 0.59). However, the simulation shows less pronounced groundwater variations. This is attributed to the use of the mean monthly net infiltration that ignores daily variations associated with high-intensity rainfall, root water uptake and actual evapotranspiration at the site. Accordingly, there are few fluctuations in the precipitation and little variation between wet and dry conditions, which is directly reflected in the lower dynamics of the simulated groundwater level. Moreover, the use of low-intensity mean monthly net infiltration results in the absence of infiltration fronts as well as a lack of water perching due to soil layering, which also contributes to the reduced fluctuation of the simulated groundwater level. In addition to the overall decreasing groundwater levels from top to mid-part of the slope, the simulation captures the remarkably different dynamics of the groundwater level for nearby cross-sections as a result of variation in depth and topography of soil layers. This is consistent with the measured values at the site. Further, the maximum mean groundwater level was observed in January, February and March. This is in agreemen<sup>t</sup> with the findings of Schmidt [38], who found that most slope movements occurred during the rainy spring.

In the next step, the simulated soil water content of the top 0.5 m for the period from September 2016 until May 2017 was compared to the measured values from SoilNet. The results show that the seasonal variation as well as the variation with depth is much smaller for the simulated soil water content. The low variation in the simulated water content with depth is mainly attributed to the implementation of the low-intensity mean monthly net infiltration in which the fluctuation in precipitation and the extremely rainy and dry periods are moderated. Notably, the use of a mean monthly low-intensity net infiltration results in the absence of infiltration fronts. In combination with the homogeneous soil hydraulic properties within each layer and the lack of depth-dependent root water uptake, this results in highly simplified and incorrect water content distributions with depth. Therefore, it was concluded that the simulation results do not seem to capture relevant features of the measured near-surface water content distribution, which is the area of interest for slope stability evaluation. Thus, it is evident that the evaluation of slope stability in response to intensive rainfall events cannot be based on monthly net infiltration. Therefore, the measured soil water content data will be considered directly for slope stability evaluation. As the ERT monitoring showed very little dynamics in the deeper layers, soil water content below the top 0.5 m was taken from the simulations using lowintensity mean monthly net infiltration.

#### **5. Model Results for Precipitation Events**

To study the slope stability in response to precipitation, two days were selected from the available SoilNet data as the initial conditions for the simulation of a hypothetical rainfall event. These two days were chosen so that both dry and wet initial conditions are considered. According to the values of *Inet* [mm month−1] within the SoilNet measurement period, the driest day occurred at the end of September 2016 after a three-month period with negative net infiltration. The wettest initial conditions occurred at the end of March 2017 after an extended period with positive net infiltration starting in October 2016. Accordingly, the data from 30 September 2016 and 31 March 2017 were used to define dry and wet initial conditions, respectively. The measured soil water content was horizontally quite variable due to the heterogeneity of the soil hydraulic properties within each soil layer. Since the soil layers were assumed to be homogeneous in the simulations, some data points exceeded the saturated water content for the wet conditions. Therefore, the measured soil water content was normalized with 0.395 cm3cm−<sup>3</sup> as the maximum water content because a full saturation at a water content of 0.40 cm3cm−<sup>3</sup> resulted in numerical issues. In order to ensure numerical stability, it was also required that the water content distribution with depth varied by at least 0.02 cm3cm−<sup>3</sup> between depth. Therefore, some normalized soil water content values were manually adjusted. The measurements in dry conditions (30 September 2016) were not normalized and used as is. The measured values were linearly interpolated between the sensor locations by Comsol Multiphysics.

The saturation and simulated LFS distributions before the start of the event rainfall on 30 September 2016 and 31 March 2017 are shown in Figure 10. As expected, the most failure-prone areas besides the scarp area appeared in the middle part of the slope (i.e., locations B and C), where the groundwater was close to the surface and the water content was higher. The LFS at position A in the scarp area was low and close to 1 due to the geometry of the slope with a high inclination in this area. In the following, the simulation results for the two most failure-prone locations B and C will be discussed, because the slope instability at those spots is attributed to water level changes. The initial water content distribution along the depth for locations B and C of the mentioned days is shown in Figure 11. For both days, the groundwater level and the soil water content are higher at location B than at location C. In case of the wet conditions, the groundwater level is similar at both locations. The sharp changes in water content below −2 m are related to soil layering.

**Figure 10.** The initial saturation (**<sup>a</sup>**,**b**) and LFS distribution (**<sup>c</sup>**,**d**) for the rainfall simulations on 30 September 2016 and 31 March 2017.

**Figure 11.** Initial conditions for water content for the event rainfall on 30 September 2016 and 31 March 2017 at the two most vulnerable locations B and C indicated in Figure 10.

In the final step, simulated slope stability in response to an event rainfall of 20 mm h−<sup>1</sup> is presented for the two selected days. The simulated LFS profiles at location B and C for the dry initial conditions have a maximum difference at the soil surface because of the difference in water content (Figure 12). Below a depth of −0.5 m, the differences in LFS are associated with the differences in the mean monthly conditions of the soil. Accordingly, water content and pressure decreased uniformly from the groundwater table upwards below the depth of −0.5 m. As discussed before, this is attributed to the homogeneous soil layers with no root water uptake and the implemented low-intensity mean monthly net infiltration to obtain the initial conditions below the depth of −0.5 m. The groundwater level and the soil water content below the depth of −0.5 m in March 2017 was higher than in September 2016, which resulted in a relatively lower LFS in March 2017 for the equivalent locations.

**Figure 12.** Variation in LFS during the event rainfall for two days (30 September 2016 and 31 March 2017) with different initial conditions, comparably dry and wet soil, for the two most vulnerable locations B and C indicated in Figure 10.

Taking these two sets of initial conditions, the response to an event rainfall of 20 mm h−<sup>1</sup> was simulated. Figure 12 shows the temporal development of the LFS for the two locations during the rainfall until a LFS of 1 is reached. The critical amounts of rainfall needed are 60 mm and 230 mm for the wet and dry initial conditions, respectively. It can be seen that significantly less rainfall is required on 31 March 2017 with the overall wetter initial condition compared to 16 September 2016. These results are consistent with previous studies [71,72] that have shown that initial hydrological conditions play an important role in the timing of failure initiation. It is interesting to note that the instability threshold was reached first in location C for both wet and dry initial conditions, although the top 0.5 m of the soil at the beginning of the rainfall was drier at location C for the both dry and wet initial conditions. Within a short time (<2 h) after the start of the rainfall, the near surface water content at location C becomes higher than that at location B. Consequently, location C reaches the failure threshold earlier. This is attributed to the difference in bedrock topography and soil depth at these two locations. In particular, the depth of the less permeable midand base layer is only 0.4 m from the surface at location C compared to 2.2 m at location B. Accordingly, the local pore pressure and saturation increase faster at this location, which

means that a potential failure state is reached after less rainfall. This is in agreemen<sup>t</sup> with the findings of Moradi et al. [13] that specifically highlighted the importance of bedrock topography and soil layering for slope stability evaluation.
