Combining Site Characterization, Monitoring and Hydromechanical Modeling for Assessing Slope Stability

Abstract

Rainfall-induced landslides are a disastrous natural hazard causing loss of life and significant damage to infrastructure, farmland and housing. Hydromechanical models are one way to assess the slope stability and to predict critical combinations of groundwater levels, soil water content and precipitation. However, hydromechanical models for slope stability evaluation require knowledge about mechanical and hydraulic parameters of the soils, lithostratigraphy and morphology. In this work, we present a multi-method approach of site characterization and investigation in combination with a hydromechanical model for a landslide-prone hillslope near Bonn, Germany. The field investigation was used to construct a three-dimensional slope model with major geological units derived from drilling and refraction seismic surveys. Mechanical and hydraulic soil parameters were obtained from previously published values for the study site based on laboratory analysis. Water dynamics were monitored through geoelectrical monitoring, a soil water content sensor network and groundwater stations. Historical data were used for calibration and validation of the hydromechanical model. The well-constrained model was then used to calculate potentially hazardous precipitation events to derive critical thresholds for monitored variables, such as soil water content and precipitation. This work introduces a potential workflow to improve numerical slope stability analysis through multiple data sources from field investigations and outlines the usage of such a system with respect to a site-specific early-warning system.

1. Introduction

Landslides are common natural hazards in many parts of the world, often triggered by increased pore pressure after heavy rainfall. Due to climate change, more intense rainfall events are expected and the frequency of destructive landslides may increase [1,2,3]. In order to evaluate hazards associated with landslide-prone hillslopes, several modeling approaches can be used to determine critical precipitation events that may lead to slope failure. Common limit-equilibrium models tend to overestimate the factor of safety as a measure of slope stability in complex geometrical setups [4]. The concept of a local factor of safety in Coulomb stress-field based finite element models is one approach used to overcome this limitation [5]. Studies applying limit-equilibrium and continuous finite element models have partly found good agreement between the results for stability analysis e.g., [6]. However, variation between models depends strongly on the methods used and the application scenarios [7,8].

The incorporation of spatial variability of soil and rock types into models for slope stability evaluation is crucial to assess the structural and hydrological state of a hillslope e.g., [9]. Slope morphology and spatial distribution of the material properties affect the slope failure potential e.g., [10,11]. Especially bedrock topography and soil depth were identified as important factors with respect to slope hydrology and stability [9,12,13]. Geophysical characterization e.g., [14,15,16] and monitoring e.g., [17,18] are increasingly used for structural and hydrological assessment of landslide-prone hillslopes. Such geophysical studies often combine multiple methods to study the subsurface structures and seismics and electric resistivity tomography (ERT) are the most commonly used e.g., [19,20]. Seismic refraction is particularly useful to identify lithological layers and slip surfaces e.g., [21,22], whereas ERT is able to provide information on the water content distribution in the subsurface e.g., [23,24] by using the correlation between bulk electrical resistivity and saturation e.g., [25,26]. Compared to drilling methods or point sensors, geophysical measurements usually provide information with a higher spatial resolution at lower cost and are only minimally invasive [27]. The combination of geophysical methods allows collection of complementary information and can be analyzed using data fusion methods [28]. Typically, supporting laboratory studies are used to improve geophysical monitoring concepts and to evaluate innovative geophysical methods, such as self-potential measurements, for detecting critical hydrological conditions [29,30]. An extensive review of geophysical monitoring methods for failure-prone hillslopes is given in Whiteley et al. [31]. In addition, the development of cost-effective sensor networks for monitoring soil water content and slope movement has gained momentum in recent years, which allows us to bridge the gap between costly boreholes and extensive geophysical monitoring and surveying e.g., [32,33].

A few case studies have combined hydrogeological and geomorphological site characterization, geophysical monitoring and hydromechanical modeling. For shallow landslides in pyroclastic soils, water content, derived from electric resistivity profiles, were combined with statistical modeling using a cellular automaton to derive a relationship between the factor of safety of a hillslope with in situ measurable quantities [34]. For the La Clapiere landslide in France, vertical electrical sounding was used to obtain the underground structure used in a geomechanical model [35]. An extensive ERT survey was combined with groundwater measurements and meteorological data to study the groundwater dynamics of a translational landslide and to develop a conceptual model [36]. To assess the slope stability of the Brzozowka landslide near Cracow (Poland), ERT monitoring was combined with drilling and laboratory tests to improve a stability analysis in which critical conditions for slope stability were derived from simulations for extreme precipitation events [37].

In this study, we present a slope stability analysis of a failure-prone hillslope in Germany based on the combination of thorough site investigation and monitoring with a hydromechanical finite element model. The site investigation includes a geological, geomorphological, hydrological, and geophysical characterization and soil water content and groundwater monitoring using a sensor network as well as geoelectric measurements. We introduce a work flow incorporating the different data sources into the model setup and for validating the model results. The data used in this study are partly obtained from former studies at the study site, and include (i) groundwater level [38], (ii) borehole logs [38], (iii) laboratory tests for soil hydraulic and mechanical properties [38,39], (iv) slope movement and (v) precipitation [40], as well as recent (vi) seismic refraction surveys, (vii) electric resistivity tomography, (viii) soil water content monitoring network and (ix) a digital elevation model. The parameterized, constrained and verified hydromechanical model was subsequently used to study potential hazardous precipitation events. We compared the results to the observations in the field and discuss how the developed approach could be extended towards a site-specific early warning system.

2. Study Area

The study site is a failure-prone hillslope at the Dollendorfer Hardt located 16 km south-east of the city of Bonn in the Siebengebirge, Germany (Figure 1). This site has been investigated in several previous studies since the 1990s [39,40,41]. At least two major landslides occurred at the study site in the past 100 years (1958 and 1972), and minor movements have been observed since the last failure event [39,40]. There is anecdotal evidence that attributes the initiation of the first major landslide to the construction of a trail at the upper part of the current scarp area. However, intense rainfall is suspected to be the triggering factor of both major events [38].

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].

The scarps of the two landslides at the study site cover an area of almost 30,000 m${}^2$ 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 $\pm 3$ cm year${}^{-1}$ 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 suggest that the seismic wave velocity range associated with the three layers are <300 m s${}^{-1}$, 400–600 m s${}^{-1}$, and >800 m s${}^{-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:

• a top layer consisting of clayey sediments, trachyte tuff and loess, transported and mixed by the landslides;
• an intermediate layer of Tertiary sediments mainly consisting of silt and clay;
• a base layer of Devonian bedrock, strongly weathered at the top.

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.

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.

The net infiltration, ${\mathit{I}}_{\mathit{n}\mathit{e}\mathit{t}}$ [mm month${}^{-1}$], is the difference between precipitation P [mm month${}^{-1}$] and actual evapotranspiration $\mathit{E}\mathit{T}$ [mm month${}^{-1}$]:

$${\mathit{I}}_{\mathit{n}\mathit{e}\mathit{t}}=\mathit{P}-\mathit{E}\mathit{T}$$

(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

$$\mathit{P}\mathit{E}\mathit{T}=16\left(\frac{{\mathit{L}}_{\mathit{d}}}{12}\right)\left(\frac{\mathit{N}}{30}\right){\left(10\frac{{\mathit{T}}_{\mathit{a}}}{\mathit{l}}\right)}^{\mathit{a}}$$

(2)

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

(3)

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

(4)

where $\mathit{P}\mathit{E}\mathit{T}$ [mm month${}^{-1}$] is the estimated potential evapotranspiration, ${\mathit{L}}_{\mathit{d}}$ [h] is the average daylength of the month, N [-] is the number of days of the month, ${\mathit{T}}_{\mathit{a}}$ [$°$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 $\mathit{P}\mathit{E}\mathit{T}$ was estimated to be 644 mm for the Dollendorfer Hardt. This value is in good agreement with $\mathit{P}\mathit{E}\mathit{T}$ 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 $\mathit{P}\mathit{E}\mathit{T}$ for the study area is shown in Figure 4. It can be seen that monthly $\mathit{P}\mathit{E}\mathit{T}$ 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, $\mathit{P}\mathit{E}\mathit{T}$, 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${}^{-1}$ 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${}^{-1}$ 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 agreement 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 agreement 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 $\mathit{T}(\mathit{t},\mathit{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 $\Delta \mathit{R}$ [$\Omega$] for the measured resistance R [$\Omega$] [54]. The absolute error was set to a = 0.001 $\Omega$ 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 }_{\mathit{b}}=\frac{{\mathit{S}}^{\mathit{i}}}{{\mathit{F}}_{\mathit{f}}}\left({\sigma }_{\mathit{w}}+\frac{{\sigma }_{\mathit{s}}}{\mathit{S}}\right)$$

(5)

with saturation, S [-], water conductivity, ${\sigma }_{\mathit{w}}$ [S m${}^{-1}$], surface conductivity, ${\sigma }_{\mathit{s}}$ [S m${}^{-1}$], formation factor, ${\mathit{F}}_{\mathit{f}}$ [-] and saturation exponent, i [-]. Water conductivity, ${\sigma }_{\mathit{w}}$ = 0.1 [S m${}^{-1}$], was gained from in situ measurements in the boreholes. As surface conductivity depends on the clay content, the empirical relationship between clay content, $\mathit{c}\mathit{c}$ [%] and the surface conductivity, ${\sigma }_{\mathit{s}}$ [S m${}^{-1}$], in the work of Rhoades et al. [58] was used:

$${\sigma }_{\mathit{s}}=(2.3\times \mathit{c}\mathit{c}-2.1)\times {10}^{-3}$$

(6)

where $cc$ is the clay content ($\mathit{c}{\mathit{c}}_1$ = 57% and $\mathit{c}{\mathit{c}}_2$ = 40%, where the index 1 denotes the upper layer and 2 the lower layer) obtained from Schmidt [38]. The formation factor ${\mathit{F}}_{\mathit{f}}={\phi }^{-\mathit{j}}$ was calculated based on the porosity $\phi$ [-] and the cementation exponent $\mathit{j}=2$, where the porosity values (${\phi }_1$ = 55%, ${\phi }_2$ = 40%) were also taken from Schmidt [38] and the cementation exponent as well as the saturation exponent $\mathit{i}=2$ were based on literature values e.g., [59].

The derived water content values from ERT (Figure 5) were in agreement 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${}^3$ cm${}^{-3}$ 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${}^{-1}$ over several days resulted in increasing saturation in the upper 0.50 m of the soil.

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 \mathit{K}\left(\mathit{h}\right)\nabla \mathit{H}+\mathit{W}=\frac{\partial \theta \left(\mathit{h}\right)}{\partial t},$$

(7)

with the hydraulic head, H [m], pressure head, h [m], hydraulic conductivity, $\mathit{K}\left(\mathit{h}\right)$ [m s${}^{-1}$], volumetric water content, $\theta \left(\mathit{h}\right)$ [m${}^3$ m${}^{-3}$] and possible source/sink terms, W [s${}^{-1}$]. Following the model of van Genuchten [60], the water retention curve is given by

$${\mathit{S}}_{\mathit{e}}=\frac{\theta \left(\mathit{h}\right)-{\theta }_{\mathit{r}}}{{\theta }_{\mathit{s}}-{\theta }_{\mathit{r}}}={\left(1+{\left(\alpha \right|\mathit{h}\left|\right)}^{\mathit{n}}\right)}^{-\mathit{m}}$$

(8)

with effective saturation, ${\mathit{S}}_{\mathit{e}}$ [-], residual and saturated water content, ${\theta }_{\mathit{r}/\mathit{s}}$ [m${}^3$ m${}^{-3}$], and soil specific parameters, $\alpha$ [m${}^{-1}$], n and $\mathit{m}=1-1/\mathit{n}$ [-]. The hydraulic conductivity is based on the effective saturation

$$\mathit{K}\left(\mathit{h}\right)={\mathit{K}}_{\mathit{s}}\sqrt{{\mathit{S}}_{\mathit{e}}}{\left(1-{(1-{\mathit{S}}_{\mathit{e}}^{1/\mathit{m}})}^{\mathit{m}}\right)}^2$$

(9)

with saturated hydraulic conductivity, ${\mathit{K}}_{\mathit{s}}$ [m s${}^{-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, $\nu$ [-]. In this case, momentum equilibrium is given by the total stress tensor, $\mathbf{\sigma }$ [Pa], the unit weight, $\gamma$ [N m${}^{-3}$] and the body force vector, ${\mathit{F}}_{\mathit{v}}$[-], as

$$\nabla .\mathbf{\sigma }+\gamma {\mathit{F}}_{\mathit{v}}=0$$

(10)

From the total stress tensor, the effective stress tensor, ${\mathbf{\sigma }}^{\prime }$[Pa], can be calculated based on the pore water pressure, ${\mathit{P}}_{\mathit{w}}$ [Pa], [61,62]:

$${\mathbf{\sigma }}^{\prime }=\mathbf{\sigma }-{\mathit{S}}_{\mathit{e}}{\mathit{P}}_{\mathit{w}}\mathbf{I}$$

(11)

with the unit vector, I[-]. To assess the stability of a hillslope, the local factor of safety, $\mathit{L}\mathit{F}\mathit{S}$ [-], is calculated based on the Mohr–Coulomb failure criteria

$$\mathit{L}\mathit{F}\mathit{S}={\displaystyle \frac{cos{\varphi }^{\prime }({\mathit{c}}^{\prime }+{\sigma }_{\mathit{m}}tan{\varphi }^{\prime })}{\mathit{R}}}$$

(12)

with effective friction angle, ${\varphi }^{\prime }$ [$°$], effective cohesion, $c^{\prime }$ [Pa], mean effective stress, ${\sigma }_{\mathit{m}}=0.5\times ({\sigma }_1+{\sigma }_3)-{\mathit{S}}_{\mathit{e}}{\mathit{P}}_{\mathit{w}}$ [Pa] and the radius of the Mohr circle, $\mathit{R}=0.5\times ({\sigma }_1-{\sigma }_3)$ [Pa]. Here, ${\sigma }_{1,3}$ [Pa] are the major and minor principal stresses. A slope is considered to be prone to failure for $\mathit{L}\mathit{F}\mathit{S}<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.

4.1. Model Setup

Surface topography was taken from a digital terrain model with a spatial resolution of 1 m${}^{-2}$ 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${}^{-1}$, 600 ms${}^{-1}$ and 800 ms ${}^{-1}$ were selected as limit values for the lithological layers in agreement 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].

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).

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

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

Description	Unit	Trachyte Tuff	Tertiary Clay	Devonian Clay/Silt
Sand	%	26	11	3
Silt	%	40	41	64
Clay	%	34	48	33
${\rho }_{\mathit{b}}$	kg m${}^{-3}$	1900	2000	1900

and

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

Symbol	Unit	Trachyte Tuff	Tertiary Clay	Devonian Clay/Silt
${\theta }_{\mathit{s}}$	‒	0.40	0.35	0.40
${\theta }_{\mathit{r}}$ *	‒	0.06	0.07	0.065
$\alpha$ *	m${}^{-1}$	1.9	2.0	1.1
n *	‒	1.22	1.18	1.31
${\mathit{K}}_{\mathit{s}}$	m s${}^{-1}$	${10}^{-6}$	${10}^{-9}$–${10}^{-7}$	${10}^{-10}$–${10}^{-6}$ **
${\mathit{K}}_{\mathit{s}}$ ***	m s${}^{-1}$	${10}^{-6}$	${10}^{-7}$	$8.0\times {10}^{-9}$
${\rho }_{\mathit{b}}$	kg m${}^{-3}$	1900	2000	1900
E **	MPa	15	15	30
$\nu$ **	‒	0.35	0.35	0.35
${\varphi }^{\prime }$	$°$	34	32	30
${\mathit{c}}^{\prime }$	kPa	20	10	30

* 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⁻¹.

. 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${}^{-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, $\nu$) were estimated from typical ranges provided in the literature [68,69] (

Table 3. Typical values of Young’s modulus for cohesive material (MPa) [68,69].

Soil Type	Very Soft to Soft	Medium	Stiff to Very Stiff	Hard
Silts with slight plasticity	2.5–8	10–15	15–40	40–80
Silts with low plasticity	1.5–6	6–10	10–30	30–60
Clays with low-medium plasticity	0.5–5	5–8	8–30	30–70
Clays with high plasticity	0.35–4	4–7	7–20	20–32

). 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, $\nu$ [-], 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.

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${}^2$ > 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 agreement 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.

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 (R${}^2$-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 agreement 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 low-intensity 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 ${\mathit{I}}_{\mathit{n}\mathit{e}\mathit{t}}$ [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 cm${}^3$cm${}^{-3}$ as the maximum water content because a full saturation at a water content of 0.40 cm${}^3$ cm${}^{-3}$ 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 cm${}^3$ cm${}^{-3}$ 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.

In the final step, simulated slope stability in response to an event rainfall of 20 mm h${}^{-1}$ 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.

Taking these two sets of initial conditions, the response to an event rainfall of 20 mm h${}^{-1}$ 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 mid- and 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 agreement with the findings of Moradi et al. [13] that specifically highlighted the importance of bedrock topography and soil layering for slope stability evaluation.

6. Discussion

The coupled hydromechanical model of Lu et al. [5] was applied for slope stability assessment of the hillslope at the Dollendorfer Hardt (Bonn, Germany) with a relatively complex geometry and heterogeneity in material properties. The test site has a long history of investigations and slope analysis. Drilling, groundwater monitoring, laboratory tests, geophysical, and geomorphological surveys as well as a soil water content sensor network were used to design and parameterize the hydromechanical model. The soil water content obtained from the sensor network was successfully combined with the hydromechanical model to provide realistic initial conditions for the near-surface water content distribution. The ERT monitoring was only applied on a segment of the hillslope and was not used to define initial conditions for the simulations. As the ERT revealed very little changes in the water content distribution at greater depth even during summer time or heavy rainfall events, the sensor network covers all relevant shallow water dynamics. The results of the simulations show that both the lithological layers and the initial conditions play an important role in the redistribution of the pore water pressure and thus determine the position of the potentially unstable locations. The initially wetter conditions require significantly less rainfall than the drier initial conditions to reach potentially unstable conditions. Instabilities develop at locations that facilitate the accumulation of water due to the subsurface topography of less permeable layers. The obtained model allows us to study the influence of initial conditions and precipitation events on the slope stability.

For the test site itself, a significant mass movement seems unlikely as only potentially unstable locations close to the surface were identified in the model. However, the modeling results clearly identify the regions of most significant movement during long-lasting rainfall events at reasonable precipitation rates for the region in agreement with previous studies [38,40]. The hillslope is possibly moving in response to heavy precipitation but at a low rate of approximately 1–3 mm year${}^{-1}$ [40]. Modeling of more severe precipitation events would be required for a complete hazard assessment, also incorporating expected future changes in precipitation patterns due to climate change. Further, the model results indicate that the scarp is continuously unstable due to its high inclination. This is consistent with field observations of ruptures of small sizes that occur regularly. A slow extension of the scarp towards the top of the hill seems likely, potentially compromising hiking tracks.

In principle, the presented workflow for slope stability assessment can be extended towards a model-based early-warning system as almost all available data sources were incorporated into the model. So far, site-specific early warning systems are predominantly sensor-based with warning thresholds derived from conventional stability analysis e.g., [73,74,75]. Incorporating precipitation forecasts, a potentially critical state of a hillslope could be calculated based on the current water content distribution using the current state as initial conditions in the numerical model. However, computational demands are often the bottleneck to include near-real time complex numerical simulations into early warning systems. Therefore, it is challenging to calculate slope stability for a predicted rainfall based on the most recent soil water content measurements. A well-constrained model as presented in this work can be used to derive thresholds for implemented sensor networks or geophysical monitoring based on pre-calculated scenarios. Besides such a model-centered early-warning system, multiple data sources could be used for a more robust early-warning system by acknowledging the value of the in situ measurements itself. Data robustness could be improved through the combination of multiple different types of sensors [70]. In this work, soil water content sensors and ERT were combined to monitor the water content distribution in the subsurface. Through redundancy, erroneous measurements of single sensors or electrodes can be detected in order to reduce the risk of false alarms. Additionally, extensometers or GPS sensors could be used to capture slope movement [75]. By combining precipitation forecasts, soil water content sensor networks and geophysical monitoring to capture water content dynamics and redistribution, and model-derived thresholds, it may be possible to establish a multi-level early warning system with high accuracy.

This study showed the rich possibilities that arise through the combination of various survey and monitoring methods with a hydromechanical model for slope stability. The complementary data sources allow for constraining the model in most aspects and reduce uncertainty in the model design. However, as seen in the comparison with the soil water content sensor network, there is a small scale heterogeneity that can not be resolved using geophysics as the contrast in geophysical parameters, e.g., resistivity or wave velocity, is too small. The high clay content in the soil made the application of ground penetrating radar, as another good methodology to detect near surface heterogeneity, impossible. On other hillslopes, this may be a useful addition to reveal near-surface heterogeneity. With additional near-surface soil sampling and extensive laboratory testing, it may be possible to describe the variability in parameters of the upper soil layers and to include the variability in the numerical model using geostatistical methods. For the studied test site, the assumption of homogeneous soil layers was considered to be adequate as slope topography and layering were the dominant factors controlling slope stability. The added value of ERT monitoring was not high in this study, but this method is expected to be more useful for slopes with stronger water content dynamics at greater depth. This would require an elaborate ERT setup with multiple electrodes and various electrode spacings. In more dynamic scenarios, material parameters and also soil layering may also change within an observation period. In the introduced framework, this could be included through an updated model.

In the future, an approach considering data assimilation seems desirable, especially in the context of a continuous monitoring of the test site. This should preferably be explored with a more sophisticated numerical model for slope stability analysis. In particular, a more advanced coupling of surface flow and infiltration into the soil would significantly improve the model. Field observations showed that strong precipitation events with more than 20 mm h${}^{-1}$ only resulted in small increases in soil water content, since a large amount of the water ran off as surface water. Surface morphology supported this observation with visible runoff channels at the surface. Modeling surface runoff and ponding of water based on the surface morphology might alter infiltration dynamics along the slope during those rare heavy precipitation events. Daily water content dynamics along the hillslope are also altered by root-water uptake. While in the current study, net infiltration was considered, a more complex model for root-water uptake would change the water content distribution with depth and increase the heterogeneity of the water content distribution in the subsurface. If this process is adequately considered, the hydromechanical model could be operated with a daily or even higher resolution for the meteorological boundary conditions. In addition, differences between overgrown and bare parts of the slope could be incorporated into the model and possibly linked to the effective cohesion, as roots are known to contribute to the slope stability. With growing computational power, the full three-dimensional geometry of the hillslope and adjacent slopes could be modeled to eliminate simplifying assumptions and effects resulting from the geometrical reduction. An extension towards the calculation of plastic deformation would allow to determine slope movement and water dynamics for situations with a factor of safety larger than one, which not necessarily result in immediate rapid slope failure. Additional model verification could be achieved through measurements of slope movement.

7. Conclusions

A workflow to design, parameterize and validate a hydromechanical model using various data sources from site characterization, monitoring and laboratory tests was presented. Through the combination of complementary data sources, a detailed numerical model was implemented despite the site heterogeneity. Monitoring with point-scale sensors revealed small scale heterogeneity in the hydraulic properties of the upper soil layer, which was not resolved during site characterization. This heterogeneity had minor influence on the slope stability analysis as it is dominated by slope and bedrock topography as well as initial conditions. The multiple data sources and the presented model can be extended towards a model-centered early warning system as well as towards an early warning system based on multiple data streams with thresholds defined through the model.