Evaluating landslide response in a seismic and rainfall regime: a case study from the SE Carpathians, Romania

There have been many studies exploring rainfall-induced slope failures in earthquake-affected terrain. However, studies evaluating the potential effects of both landslide-triggering factors – rainfall and earthquakes – have been infrequent despite rising global landslide mortality risk. The SE Carpathians, which have been subjected to many large historical earthquakes and changing climate thus resulting in frequent landslides, comprise one such region that has been little explored in this context. Therefore, a massive (∼9.1 Mm²) landslide, situated along the river Bâsca Rozilei, in the Vrancea seismic zone, SE Carpathians, is chosen as a case study area to achieve the aforesaid objective (evaluating the effects of both rainfall and earthquakes on landslides) using slope stability evaluation and runout simulation. The present state of the slope reveals a factor of safety in a range of 1.17–1.32 with a static condition displacement of 0.4–4 m that reaches up to 8–60 m under dynamic (earthquake) conditions. The groundwater (GW) effect further decreases the factor of safety and increases the displacement. Ground motion amplification enhances the possibility of slope surface deformation and displacements. The debris flow prediction, implying the excessive rainfall effect, reveals a flow having a 9.0–26.0 m height and 2.1–3.0 m s⁻¹ velocity along the river channel. The predicted extent of potential debris flow is found to follow the trails possibly created by previous debris flow and/or slide events.

1 Introduction

Landslides, though a normal process of hillslope erosion, pose socio-economic risk to human life and infrastructure (Gupta et al., 2017; Froude and Petley, 2018; Pollock and Wartman, 2020; Kumar et al., 2021). Despite the rising global landslide mortality risk, effective evaluation of disastrous influences of landslides has been infrequent (Sassa, 2015; Haque et al., 2019; Klimeš et al., 2019). Such evaluation approaches could be regional (susceptibility/hazard/risk/vulnerability) or local (slope stability, runout prediction, monitoring/change-detection mapping) (Fell and Hartford, 1997; Van Westen et al., 2006; Margottini et al., 2013; Hungr, 2018). However, effectiveness in such approaches cannot be justified until the main landslide-triggering factors – rainfall and earthquakes – are evaluated together. Despite the numerous case studies of rainfall-induced slope failures in earthquake-affected terrain (Lin et al., 2006; Helmstetter and Garambois, 2010; Tang et al., 2011; Durand et al., 2018; Bontemps et al., 2020), studies predicting the potential effects of both factors have been relatively rare. The necessity of such studies becomes more critical in view of an annual average of >4000 landslide-related deaths worldwide in the last decade (Pollock and Wartman, 2020).

Figure 1 Study area. Inset (a) (source: NOAA NCEI, USA) and (b) (after Ustaszewski et al., 2008) highlight the position of the study area. Geological setting is based on Murgeanu et al. (1965), Tischler et al. (2008), and Pospíšil et al. (2012).

Owing to the capability to represent the progressive deformation in the slope under various loading conditions, numerical-modeling-based analysis can be considered one of the few approaches for effective evaluation of slope instability and associated risk (Jing, 2003; Fenton and Griffiths, 2008). Though the continuum-modeling-based approaches have been common for local-scale evaluation of hillslope response (Griffiths and Lane, 1999; Jamir et al., 2017; Kumar et al., 2018, 2021), their limitations in estimating large strain, particularly during dynamic analysis, make the discontinuum modeling a better option (Havenith et al., 2003; Bhasin and Kaynia, 2004). Apart from the stability evaluation, prediction of potential runout during the slope failure constitutes a principal risk evaluation approach (Hungr et al., 1984; Hutter et al., 1994; Rickenmann and Scheidl, 2013). Among different types of landslide, debris flows have been shown to have the maximum outreach, cause more fatalities, and cause secondary effects like river damming and subsequent outburst flooding (Jakob et al., 2005; Ding et al., 2020; Kumar et al., 2021). Among different runout prediction approaches, dynamic-model-based Rapid Mass Movement Simulation (RAMMS) (Christen et al., 2010), FLO-2D (O'Brien et al., 1993), and MassMov2D (Beguería et al., 2009) have been the more useful ones (Rickenmann and Scheidl, 2013; Kumar et al., 2021).

In view of these understandings, the present study aimed to infer the potential response of a landslide slope under seismic and extreme rainfall conditions using stability evaluation and runout simulation. Such simulations/modeling outputs depend upon certain input parameters and criteria, the values of which might be affected by uncertainties due to nonlinear behavior of material. Therefore, a parametric analysis is also performed to evaluate uncertainty. In order to achieve the aforementioned objectives, a massive (∼9.1 Mm²) landslide in the Vrancea seismic zone, SE Carpathians, is chosen as a case study area. The region has been subjected to frequent earthquakes and relatively wet climatic conditions that induce frequent landslides and related socio-economic losses (Micu et al., 2013, 2016; Micu, 2019; Mreyen et al., 2021).

Figure 2 Landslide features. (a) Landslide marked with different features. (b) Crown portion. (c) Right flank. (d) Left flank. (e) Signs of failure in the flow deposits. Image source: © Google Earth.

2 Study area

2.1 Geological setting and geomorphology

The landslide is situated at 45^∘30^′23^′′ N, 26^∘25^′05^′′ E along the river Bâsca Rozilei in the SE Carpathians, Romania (Fig. 1). The earliest record of this landslide is mentioned in the geological map by Murgeanu et al. (1965). Unfortunately, no previous record and/or dating is available at present. The slope is composed of shale belonging to the Miocene thrust belt that separates the external foredeep in the north, east, and southeast from the inner Carpathian mountain ranges. Thrust faults, strike-slip faults, and folds traverse the region in and around the vicinity of landslide slope. The origin of these structural features has been related to the Eocene–Miocene collision of the ALCAPA and Tisza–Dacia plates against the Bohemian and Moesian promontories that gave rise to the Carpathians (Tischler et al., 2008). The SE part of the Carpathians, however, is still uplifting at a rate of 3–8 mm yr⁻¹ due to the foreland coupling of the converging plates (Pospíšil and Hipmanova, 2012; Maţenco, 2017).

The landslide toe along the river hosts the village of Varlaam (Figs. 1 and 2a). The landslide has a slope gradient ranging between 15–20^∘ and encompasses an area of ∼9.1 Mm². The landslide-affected area is covered by shrubs and scattered trees towards its flanks and with grasslands in the inner parts, mainly used as pastures and hayfields. The landslide crown region has a depression that might be a surficial imprint of the paleo-detachment (or depletion zone) (Fig. 2b). Near the right (or southern) flank, a seasonal flow channel (or gully) emerges near the paleo-detachment depression and finally merges at the river channel (Fig. 2c). Near the left (or northern) flank, the slope surface comprises flow relics, possibly of paleo-debris flow and/or slide events (Fig. 2d), as also inferred from loose/unconsolidated deposit at the slope toe (Fig. 2e). This flow deposit is noted to develop 100–150 m wide minor scarps (Fig. 2e). Such scarps may further grow and result in the debris flows during extreme rainfall and/or earthquake events and hence pose a risk to the nearby human settlement.

2.2 Rainfall and earthquake regime

The study area is subjected to an increasing rainfall trend (Fig. 3a). The average monthly rainfall was 50±1.6 (SE) mm during the years 1982–2019 and has increased in recent decades (2000–2019) to 53±2.3 (SE) mm (Fig. 3b). Monthly rainfall patterns further reveal higher rainfall in the months of May, June, and July (Fig. 3c). Such enhanced summer (June–July) rainfall has been related to the existing positive phase of the North Atlantic Oscillation (NAO) index that allows the strengthening of continental climate, Mediterranean retrogressive cyclones, and the Siberian High in central and southern Europe (Constantin et al., 2007; Magyari et al., 2013; Obreht et al., 2016).

Figure 3 Rainfall (RF), soil moisture (SM), and surface runoff (RF) pattern. The region from which these parameters are extracted is highlighted in Fig. 1. Dataset (except j) source: FLDAS_NOAH01_C_GL_M model (McNally, 2018). Spatial resolution of dataset: 0.1^∘ (∼10 km). Data source of (j): GPM IMERG Final Precipitation (Huffman et al., 2019). Spatial resolution of dataset: 0.1^∘. The rainfall data in (j) are available only from 1 June 2000. The SM and SR data are not available at a daily scale for the study area. The blue line (in a, d, g) indicates linear regression, and the shaded region around it refers to the 95% confidence interval. Dots in box plots refer to outliers. The red line in (j) refers to extreme rainfall (30 mm d⁻¹). The grey-shaded region in (k) refers to those months that witnessed above-average values of RF, SM, and SR.

Apart from the rainfall, soil moisture and surface runoff also show increasing trends during the years 1982–2019 (Fig. 3d and g). Though the average monthly soil moisture and surface runoff have also increased in recent decades, their trends do not follow rainfall entirely (Fig. 3e and h). This occurs due to the fact that the temporal coexistence of rainfall, surface runoff, and soil moisture depends upon the rainfall threshold and soil conditions (antecedent soil moisture). Further, the surface runoff (water, from precipitation that flows over the land surface) correlates relatively well with the rainfall, unlike the soil moisture, which retains part of the rainfall before achieving saturation and hence does not correlate well (Fig. S1 in the Supplement). This difference in correlation is further visible in the monthly pattern (Fig. 3f and i). It is of note that the surface runoff and soil moisture data are based on the FLDAS (Famine Early Warning Systems Network Land Data Assimilation System) model (McNally, 2018). This utilizes precipitation datasets and analyses like CHIRPS (Climate Hazards Group InfraRed Precipitation with Station data) and MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2) along with land cover data to derive variables like soil moisture and surface runoff.

Figure 4 Earthquake pattern. (a, b) Position of study area. (c) Depth and earthquake magnitude. (d) Distance of epicenters from landslide and earthquake magnitude. The blue line (in d) indicates linear regression, and the shaded region around it refers to the 95 % confidence interval. Data source: National Institute for Earth Physics, Romania.

Further, the daily rainfall data of the years 2000–2019 revealed 48 extreme rainfall events (Fig. 3j). “Extreme” rainfall pertains to >30 mm d⁻¹ in the region on the basis of previous studies exploring the rainfall variability (Apostol, 2008; Croitoru et al., 2016). Out of these 48 events, 28 events occurred in the last decade, particularly in the years 2005, 2007, 2010, and 2016. The debris flows and flash floods in the region in the years 2005 and 2010 (Micu et al., 2013; Grecu et al., 2017) can be related to these extreme rainfall events in the region. The years 2005 and 2010 had relatively high precipitation due to synoptic conditions that involved pressure lows and front systems moving along a SE–NW trajectory from the Mediterranean Sea and Black Sea towards central Europe and in a west-to-east direction from the Atlantic Ocean to eastern Europe. These trajectories led to severe flood and slope failure events in different parts of central and eastern Europe (Mihailovici et al., 2006; Micu et al., 2013; Grecu et al., 2017). The influence of these trajectories is also visible in the regional rainfall pattern (Fig. S2), where years 2005 and 2010 have higher rainfalls. Though the years 2007 and 2016 also had many extreme rainfall events, these years did not have as high surface runoff as years 2005 and 2010 (Fig. 3h). Notably, many conceptual and physically based models have been proposed relating the initiation of debris flow to surface runoff conditions (Simoni et al., 2020).

Further, the temporal pattern of higher values (above average) of rainfall, surface runoff, and soil moisture revealed that May–September months dominate the trend, having the majority of the events when all three variables had extremes (i.e., above average) (Fig. 3k). These above-average values refer to the monthly scale. The temporal overlapping of these variables further justifies the occurrence of debris flows and flash floods in this region in the last decade and possibility of more such events in the near future (Micu et al., 2013; Ilinca, 2014; Grecu et al., 2017; Micu, 2019).

Apart from the temporally enhanced rainfall, surface runoff, and soil moisture, the study area is also subjected to frequent earthquakes owing to its position in the Vrancea seismic zone, which is one of the most active seismic zones in Europe (Fig. 4a and b). This region has received ∼469 earthquakes (M_w≥4) during the years 1960–2019. The earthquake event cluster represents a NE–SW trend (Fig. 4b). About 75 % of the total earthquake events occurred in a depth range of 60–180 km (sub-crustal depth), and four out of five events having a magnitude ≥ 6 occurred within 60–100 km depth (Fig. 4c). The relative dominance of M≥6 earthquakes in this depth range has been related to the reverse faulting mechanism in this depth range (Radulian et al., 2007; Petrescu et al., 2019). The possible explanation of the pattern of earthquakes has been divided into the following two categories: (1) it might be associated with a descending relic ocean lithosphere beneath the bending zone of the SE Carpathians or (2) it might be associated with a continental lithosphere that was delaminated after the collision (Bokelmann and Rodler, 2014; Petrescu et al., 2019).

Figure 5 Landslide model construction. (a) Measured peak frequency distribution. Based on Cauchie et al. (2019). (b) Digital elevation model. (c) Soil (or debris) thickness pattern in the landslide. (d) Cross-sectional view of landslide model.

Though the majority of earthquakes have their magnitude smaller than 5 and quite deep (mostly between 60 and 180 km), their epicenters are situated within 50 km of the study area (Fig. 4d). Such intermediate to deep earthquakes in the Vrancea region (study area) have triggered landslides as far as 250–300 km from their epicenters (Havenith et al., 2016). Further, any major future earthquake might have ground effects in a much larger area (150 000 km²), possibly causing more landslides (Havenith et al., 2016). The regional distribution of the annual rainfall and earthquakes around the study area is also shown in Fig. S3.

3 Methodology

In order to evaluate the landslide response under seismic and extreme rainfall conditions, our approach involved dynamic slope stability analysis and runout simulation, respectively. Both the techniques required a landslide model that was constructed using field-based ambient noise analysis, empirical equations/values, a digital elevation model, and geological modeling software. Details are as follows.

3.1 Debris (or loose material) depth estimation

We analyzed seismic ambient noise at 56 measurement points to estimate the depth of impedance contrasts. The equipment was composed of seven Güralp CMG-6TD 30 s velocimeters, one Lennartz 5 s velocimeter, and one CityShark II velocimeter. The technique aims at estimating the site resonance frequency by computing the spectral ratio between horizontal (NS, EW) and vertical components (Nakamura, 2009). Under particular geological conditions where impedance contrast exists at depth, as representative of a loose/soft material overlying bedrock, the resulting horizontal-to-vertical spectral ratio (HVSR) curve presents a peak in correspondence of the site resonance frequency (f_o). Figure 5a represents the location of the inferred f_o in a range of <1.5–4.5 Hz. Lower frequencies, generally implying higher thickness of loose material, are noted in the central part and near the right flank.

The thickness (h) of the loose/soft material is consecutively estimated using the shear-wave velocity (V_s) and resonance frequency (f_o) in the following equation (Murphy et al., 1971; Ibs-von Seht and Wohlenberg, 1999):

$$\begin{array}{}\text{(1)}& h=V_{\mathrm{s}}/\left(\mathrm{4}\cdot f_{\mathrm{o}}\right).\end{array}$$

In view of the similar litho-tectonic conditions and spatial proximity, the shear-wave velocity (V_s) values in the present study are based on Mreyen et al. (2021). For the loose overburden (soil) and rock mass, the V_s values are taken as ∼400 and ∼900 m s⁻¹, respectively.

The thickness of the loose material (inferred from the HVSR and V_s) at different measurement locations was later imported into Leapfrog Geo software (v5.1) along with the surface morphology (Fig. 5b). The surface morphology with a spatial resolution of ∼12 m is based on the TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurement) digital elevation model. The surface morphology and depth information of loose material were integrated using Leapfrog Geo (v5) to construct a continuous soil thickness layer and hence a 3D model of the landslide (Fig. 5c and d). This model was later used to extract the 2D slope sections (CS-1, CS-2, CS-3, and CS-4) for the slope stability evaluation (Sect. 3.2) and runout simulation (Sect. 3.3).

Figure 6 Model configuration for the slope stability analysis. (a) Landslide model. The location of the different cross sections used in the UDEC models are marked by red lines. (b–e) Configuration of the sections: CS-1 to CS-4.

3.2 Slope stability evaluation

The 2D slope sections (CS-1, CS-2, CS-3, and CS-4), shown in Fig. 6a, were used to determine the hillslope response under static (gravity) and dynamic (seismic) conditions by performing slope stability analysis in UDEC (2014) software. Each slope section comprises loose overburden (soil) over rock mass and an interface joint separating these blocks (Fig. 6b–e).

Under static conditions, the factor of safety of the slope and potential material displacement are determined, whereas under dynamic conditions, the potential material displacement, peak ground acceleration (PGA), and spectral ratio are evaluated. The factor of safety is determined using a shear strength reduction approach (Matsui and San, 1992; Griffiths and Lane, 1999). The potential material displacement under static conditions refers to displacements after the model has reached static equilibrium under gravity load. The spectral ratios are used to understand the response of the medium to the input signal by comparing the signals obtained in the monitoring points at the surface with the signal at the monitoring point at depth (base) (McCowan and Lacoss, 1978).

For the PGA and spectral ratio, material models are considered elastic, whereas for the factor of safety and material displacement (static/dynamic) calculations, elasto-plastic models are considered. The elastic material model involved modulus (elastic/shear/bulk) values of the rock mass and soil. In elasto-plastic conditions, modified Hoek–Brown (MHB) plasticity criteria (Hoek et al., 2002) and Mohr–Coulomb (M–C) plasticity criteria (Coulomb, 1776; Mohr, 1914) are used for the rock mass and soil, respectively. The joint plane is assigned Coulomb slip criteria (Coulomb, 1776) in both elastic and plastic conditions. For dynamic analysis, two different signals, i.e., the Ricker wavelet (Ricker, 1943) and a signal record of the 1976 Friuli earthquake, are used (Fig. 8).

The Ricker wavelet, a theoretical waveform, provides the advantage of being a relatively short signal marked by energy distributed over a range of frequencies. Therefore, the PGA and spectral ratios are evaluated using the Ricker wavelet to understand the ground motion amplification on the landslide surface. Notably, in many studies such ground motion amplification is found to enhance the slope instability (Lenti and Martino, 2012; Del Gaudio et al., 2019). The Ricker wavelet has been used in several studies owing to its reliable representation of seismic waves propagating through the viscoelastic homogeneous media (Gholamy and Kreinovich, 2014). Further, the displacement is determined using both dynamic signals (Ricker wavelet and Friuli earthquake, 1976) to evaluate the difference.

Soil and rock mass blocks in the cross sections (CS-1 to CS-4) were discretized into finite-difference zones of 6 and 20 m size, respectively, according to the following relation (Kuhlemeyer and Lysmer, 1973):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }l\le \mathit{\lambda }/\mathrm{10}\mathrm{or}\le \mathit{\lambda }/\mathrm{8}.\end{array}$$

Here, Δl is zone size and λ is the wavelength associated with the dominant frequency. λ can be determined using $\mathit{\lambda }=C/f$, where C is the speed of wave propagation associated with the fundamental frequency (f). For the C (or shear-wave velocity) of soil and rock mass, we used 400 and 900 m s⁻¹, respectively (Sect. 3.1). f = 2.0–4.5 Hz was considered a central frequency range. The lateral boundaries in all four slope sections (or models) were considered free-field boundaries owing to the near-surface position of the hillslope (Fig. 6). A stress-boundary condition (Joyner and Chen, 1975; UDEC, 2014) was applied at the base in which the horizontal direction is considered viscous, whereas the vertical direction is kept free. This stress-boundary condition converts seismic input from velocity waves to stress waves. To approximate the natural attenuation in the models during the seismic loading, Rayleigh damping with a 0.02 damping ratio (i.e., 2 % fraction of critical damping) and 2.5 Hz central frequency was used with both the mass and the stiffness damping. Though most of the soil types and rock mass possess the damping in the 2 %–5 % fraction of the critical damping (Biggs, 1964), plasticity models (M–C criteria) and the presence of joints result in further energy loss (UDEC, 2014). Therefore, the damping ratio was kept at the lower level of the suggested range.

Since the area is subjected to temporally enhanced rainfall (Sect. 2.2) and some studies have noted the percolation of rainfall water in the loose material resulting in the groundwater (GW) level increase and subsequent slope instability (Van Asch et al., 1999; Liang, 2020), the effect of the GW is also explored. To simulate the GW effect, coupled hydraulic (fluid-flow)–mechanical analysis was used in which mechanical deformation and joint fluid pressure affect each other as analysis progresses. Further, the model was brought to static equilibrium before performing factor of safety (FS) calculations. Notably, FS calculations were also performed under mechanical stress only (without GW). Steady-state (water table) fluid-flow analysis was used to simulate fluid flow. The GW is included in static as well as in dynamic analysis in plasticity conditions. UDEC allows the GW simulation through the joints as per the parallel plate model (Witherspoon et al., 1980). The parameters and their values used in the static and dynamic analysis are given in Table 1.

Table 1 Input parameters used in the slope stability analysis.

¹ It was inferred from the empirical equation of Kahraman (2001) using the V_S and V_P data of Mreyen et al. (2021). ² These values were inferred from the empirical equations of McDowell (1990) using the P- and S-wave velocity of Mreyen et al. (2021). ³ Based on Hoek and Brown (1997) and field observation. ⁴ It was inferred from the empirical equations of Barton (1972) and Hoek and Diederichs (2006) using the elastic modulus of rock and approximated spacing of joint sets of ∼5–10 cm. This spacing was assumed in view of the highly sheared nature of rock mass. ⁵ Based on Bednarczyk (2018) and Peranić et al. (2020) due to similar litho-tectonic conditions. ⁶ Based on Barton and Choubey (1977). ⁷ Based on UDEC (2014).

A parametric analysis was also performed to justify the selection of values of different input parameters by evaluating the change in the output parameters in response to the change in different input parameters. Out of four slope sections, CS-2 and CS-3 were chosen to perform the parametric analysis in view of their central position in the landslide and the heterogeneity in soil thickness and topography (Fig. 6c and d). In order to understand the effect of the GW level change, two GW levels were considered in the CS-2 and CS-3 sections. Since UDEC simulates the fluid flow through a joint aperture, the GW level change is manifested by different heights (h₁, h₂) of the GW at the joint. Here, the difference in h₁ and h₂, i.e., Δh, is 10 m (Fig. 6d). Among the different input parameters listed in Table 1, the angle of internal friction of soil, joint friction angle, groundwater head, and elastic modulus were used for the parametric analysis. It is of note that the bulk and shear modulus were also changed along with elastic modulus because all three modulus parameters are interrelated (McDowell, 1990). Though each parameter might have a certain effect on the output, these four have been noted to affect the factor of safety and displacement more (Kumar et al., 2021).

3.3 Runout simulation

The hillslopes affected by the seismic shaking have also been noted to be more prone to rainfall-induced slope failures, particularly in the form of debris flows (Shieh et al., 2009; Tang et al., 2011). Such debris flows can be initiated by either increased pore pressure or runoff involving entrainment (Godt and Coe, 2007). Thus, the increased frequencies of the extreme rainfall, soil moisture, surface runoff, and recent debris flows events in the region (Sect. 2.2) escalate the possibility of debris flow in the Varlaam landslide.

To ascertain the outreach of such potential debris flow during an extreme rainfall event, Voellmy–Salm (Voellmy, 1955; Salm, 1993) fluid-flow continuum-model-based Rapid Mass Movement Simulation (RAMMS) software was used. RAMMS divides the frictional resistance into a dry-Coulomb-type friction (μ) and viscous–turbulent friction (ξ) (Christen et al., 2010). The frictional resistance S (Pa) is thus

$$\begin{array}{}\text{(3)}& S=\mathit{\mu }N+\left(\mathit{\rho }g{u}^{\mathrm{2}}\right)/\mathit{\xi },\end{array}$$

where N=ρhgcos (ψ) is the normal stress on the running surface, ρ is density, g is gravitational acceleration, ψ is the slope angle, h is flow height, and u is (u_x,u_y) consisting of the flow velocity in the x and y directions. A detailed description of the governing equations is presented in the Supplement.

Generally, the values for the μ and ξ parameters are achieved using the reconstruction of real events through simulation and subsequent comparison between dimensional characteristics of real and simulated events. However, the toe of Varlaam landslide merges with the river floor, and hence there is an uncertainty in the reconstruction of the volume of previous flow events that has been washed away by the river. Therefore, μ and ξ are taken in view of the topography of the landslide slope and runout path, landslide material, and previous studies/models (Hürlimann et al., 2008; Rickenmann and Scheidl, 2013; RAMMS v1.7.0). In this study, the maximum allowable friction (μ), i.e., μ=0.4 (or ϕ=21.8^∘), was used with the turbulence (ξ) of 250 m s⁻¹ (Table 2). Different depths were considered block release depths in view of uncertainties to ascertain the exact depth of loose material that will be eroded/entrained during the debris flow. Though the landslide surface has some relics of flow channels near the left flank (Fig. 2d and e), the data pertaining to the spatial–temporal pattern of discharge at these flow channel/gullies were not available. Therefore, the release area is chosen as block release because it has been more appropriate when the flow path (e.g., gully) and its possible discharge on the slope are uncertain (RAMMS v1.7.0). The runout stopping criterion is based on the momentum threshold, which was considered 5 % of moving mass. A sensitivity analysis is also performed to evaluate the possible influence of frictional parameters on the debris flow characteristics. Further, in order to understand the influence of river channel morphology on the debris flow characteristics, their interrelationship is also sought.

Table 2 Input parameters used in the runout simulation.

¹ Considering the fact that during slope failure, irrespective of the type of trigger, the entire loose material might not slide down, the depth is taken as a variable. ² In order to keep the results of a conservative nature, we have taken a maximum allowable friction, i.e., μ=0.4 (Hungr et al., 1984; RAMMS v1.7.0). This case is considered to understand the potential impacts of debris flow even after the maximum friction. ³ This range is used in view of the type of loose material, i.e., cohesive (RAMMS v1.7.0).

4 Results and discussion

4.1 Slope stability evaluation

4.1.1 Factor of safety (FS) and displacement

The FS of the slope varies in a range of 1.17–1.32 that decreases further to 1.09–1.29 under groundwater (GW) conditions (Fig. 8). In both cases, the CS-2 model attains the lowest FS, implying more instability. The displacement in loose material was obtained in static, static with fluid (GW), dynamic, and dynamic with fluid (GW) conditions. Under the static conditions, displacement ranges between 0.4–4.0 m, which increases to 0.68–18 m under the GW conditions with its minimum at CS-1 and maximum at CS-2 (Fig. 8). Under dynamic conditions, displacement ranges from 8–60 m and further increases to 7.5–62 m by combining dynamic with GW conditions. Similarly to the static conditions, the minimum displacement is noted at CS-1, whereas the maximum is at CS-2. Further, in all sections (CS-1 to CS-4), displacement accumulated mostly at the upper part of the debris layer (i.e., landslide crown) or at the steepest portion of the slope surface. This spatial affinity of displacement and a steep gradient is caused by the influence of topography on the material displacement (Kumar et al., 2021). Notably, this dynamic displacement pattern pertains to the Friuli earthquake signal (Fig. 7c and d). A comparison of the static and dynamic displacement (caused by the Friuli earthquake signal and Ricker wavelet) is presented in Fig. 9.

Figure 7 Seismic signals in the time and frequency domain. (a, b) Ricker wavelet (as recorded at the model base monitoring point); (c, d) 1976 Friuli earthquake (Italy). Note the different timescales.

Figure 8 Factor of safety (FS) and material displacement (x direction). Panels (a)–(d) refer to original slope sections with sub-sections (red rectangle) used to represent displacement. Panels (e)–(t) refer to displacement in static, static + GW, dynamic, and dynamic + GW conditions.

Figure 9 Comparison of material displacement under different conditions. St. and Dy. refer to static and dynamic conditions, respectively. GW refers to groundwater.

As also shown in Fig. 9, the GW conditions enhanced the displacement in static as well as in dynamic conditions (Fig. 9). Static displacement showed the least scattering as evident from the median level and least difference in max and min values. Further, except for the CS-2 section, all three sections (CS-1, CS-3, CS-4) have relatively low dynamic displacement in dry and wet (GW) conditions due to the Ricker wavelet compared to the displacement caused by the Friuli signal (Fig. 9a–d). This difference may be attributed to the response of steep topography (of the CS-2 model) to the high-energy Ricker wavelet signal (Fig. 7b).

Thus, it can be understood that the Varlaam hillslope, situated in the region having frequent extreme rainfalls and earthquakes, attains more instability under saturated dynamic conditions.

4.1.2 Parametric analysis

The factor of safety (FS) of the slope increased in response to increase in the angle of internal friction of soil, joint friction, and elastic modulus (Fig. 10). Such increase in the FS (∼7 % in the CS-2) is attained by increasing the angle of internal friction of soil. This effect is attributed to the “shear strength reduction” (SSR) approach. The GW level increase resulted in a decreasing FS because the increased GW level increased the joint flow rate and thus enhanced the fluid pressure on the overlying medium, i.e., soil. This increased fluid pressure further decreased the normal stress and hence the shear stress of the overlying soil, as per Mohr's criteria (Mohr, 1914). Such decrease in the shear stress of soil resulted in the decreased FS.

Figure 10 Parametric analysis. (a–d) Variation in the FS. (e–h) Variation in the static displacement. (i–l) Variation in the dynamic displacement. Grey bar represents the values that are used in the slope stability analysis (Sect. 4.1). Notably, vertical axes in (e)–(l) are kept the same to show relative increase in displacement under static and dynamic conditions.

Owing to their spatially variable nature, static and dynamic displacements are represented in ranges of the maximum (max) and minimum (min) in Fig. 10. Static displacement decreased on increasing the angle of internal friction of soil, joint friction, and elastic modulus. Such decrease (∼40 % in CS-2 and ∼38 % in CS-3) occurred in response to the modulus increase. This decrease in the displacement refers to the fact that the increased modulus increases the normal and shear strength of the soil, and hence displacement will decrease on increasing the modulus (Hara et al., 1974). The GW level increase resulted in the increased static displacement (∼16 % in CS-2, ∼36 % in CS-3). Such increase in the static displacement is attributed to the decreased shear strength of soil due to the increased joint fluid pressure (Witherspoon et al., 1980). Similarly to the static displacement, dynamic displacement decreased on increasing the angle of internal friction of soil, joint friction, and elastic modulus and increased on increasing the GW level. Along with the modulus, the angle of internal friction of soil is also noted to decrease (∼16 % in CS-2, ∼21 % in CS-3) the dynamic displacement more. The increase in the GW level resulted in 8 % and 33 % increases in the CS-2 and CS-3 models, respectively, in dynamic displacement.

Notably, the present study utilized approximated values of the input parameters for the slope stability analysis (Table 1). Though the approximated values cannot replace the values measured in the geotechnical analysis, parametric analysis minimizes the uncertainty caused by the selection of specific values by exploring the possible output pattern.

Thus, by utilizing the central values (highlighted as grey in Fig. 10) in the slope stability analysis (Sect. 4.1.1), the present study attempted to minimize such uncertainty in the findings. Further, though the GW was also used in the UDEC models to infer the influence of saturation on slope stability, the potential response of the slope under excessive saturation (extreme rainfall) is further explored using the debris flow runout simulation (Sect. 4.2).

4.1.3 Peak ground acceleration (PGA)

Apart from the FS and displacement, ground motion (acceleration) amplification was also evaluated to understand the potential seismic deformation at the slope surface (Fig. 11). The input seismic signal for the following acceleration pattern is presented in Fig. 7a. For all four models (CS-1 to CS-4), the PGA values at the river floor (RF) ranges between 5.78–7.47 m s⁻¹ (0.58–0.74 g), whereas at the rock mass surface above the landslide crown (CR) they vary from 6.37 to 10.19 m s⁻¹ (0.65–1.03 g) (Fig. 11). At the model base (MB), maximum acceleration remains between 3.79–3.90 m s⁻¹ (0.38–0.39 g).

Figure 11 Maximum acceleration at different monitoring points.

Thus, the PGA at the RF is amplified by ∼1.5–2.0 times the maximum acceleration at the model base, whereas at the rock mass surface above the landslide crown, it is amplified by ∼1.7–2.7 times the maximum acceleration at the model base. Such amplification of the PGA at the rock mass surface above the landslide crown can be attributed to the topographic irregularity and upward propagation of seismic waves where they meet preceding waves produced on the relatively horizontal surface of the slope (Jibson, 1987; Havenith et al., 2003; Bourdeau and Havenith, 2008; Luo et al., 2020).

The debris surface, however, attains higher PGA in all four models than the rock mass surface as noted at the following three monitoring stations: DB_Lw, DB_Md, and DB_Up (Fig. 11). At the lower part of the debris (DB_Lw), the PGA ranges from 8.3 to 12.13 m s⁻¹ (0.84–1.23 g) and further grew at the middle part of the debris (DB_Md), attaining 10.17–14.40 m s⁻¹ (1.03–1.46 g). The maximum PGA is attained by the upper part of the debris (DB_Up) with a range of 7.26–18.50 m s⁻² (0.74–1.88 g). Such relatively high PGA at the debris surface can be linked to the impedance contrast between underlying rock mass and overlying soil and/or partial loss of the shear strength during seismicity (Novak and Han, 1990; Šafak, 2001).

Detailed evaluation at different monitoring points in each model is as follows: model base (MB) and river floor (RF) monitoring points have almost similar maximum acceleration values in all four models. At the lower part of the debris, i.e., DB_Lw, higher PGA is attained by the CS-3 model (∼12.1 m s⁻²) followed by the CS-2 model (∼10.8 m s⁻¹) in comparison to DB_Low points of CS-1 and CS-4. Higher PGA is attributed to lower soil thickness below this monitoring point in the CS-3 and CS-2 models that could be the main reason for acceleration amplification as also stated by Murphy et al. (1971) and Beresnev and Wen (1996).

At the middle part of the debris, i.e., DB_Md, higher PGA is attained by the CS-2 model (∼14.4 m s⁻²). Notably, despite the higher soil thickness, this monitoring point obtained a higher PGA. It possibly occurred due to irregular topography of the CS-2 model that generally results in interference of direct and scattered waves and hence amplification of ground motions (Asimaki and Mohammadi, 2018).

At the upper part of the debris, i.e., DB_Up, higher PGA is attained by the CS-1 model (18.5 m s⁻²) followed by the CS-4 model (15.8 m/s⁻²). The effect of soil thickness below this monitoring point, as explained for the lower part of debris, could be the main reason for such amplification at this monitoring point in these models. The monitoring point at the rock mass surface above the landslide crown (CR) also has almost similar PGA values in all the models except the CS-3 model. Higher PGA (10.19 m s⁻²) at the CR monitoring point of the CS-3 model might be due to its position on a steeper surface, whereas CR points in other models are at a relatively flat surface.

4.1.4 Spectral ratio

The ground motion amplifications were also explored using the spectral ratios at two central slope sections: CS-2 and CS-3 (Fig. 12). In both models, the RF (river floor) point showed no significant amplification at any particular frequency, possibly due to the flat surface positioning. In the CS-2 model, the debris lower part (DB_Lw) point shows notable amplification at 2.0–2.5 Hz with minor amplification at 4.5–5.0 Hz, whereas in the CS-3 model, the DB_Lw point shows attenuation (or de-amplification) near ∼2 Hz and slight amplification at 4.5–6.0 Hz. The contrast of amplification and de-amplification at ∼2 Hz is attributed to the geometrical variation in topography because the DB_Lw point in CS-2 is situated at a relatively elevated surface, whereas in CS-3, it is at a relatively shallow surface. Minor geometrical variations at the slope toe have been also observed to result in de-amplification at low frequencies in other studies (Bouckovalas and Papadimitriou, 2005).

Figure 12 Spectral ratio pattern. (a) CS-2 model with the position of monitoring points. (b) Spectral ratio pattern in CS-2. (c) The CS-3 model with the position of monitoring points. (d) Spectral ratio pattern in CS-3.

Notably, along with the DB_Lw point, the debris middle part (DB_Md) and debris upper part (DB_Up) points in both the models also have minor/major amplification at 4.5–6.0 Hz. This coexistence of amplification at a certain frequency range by different monitoring points at the debris surface may be attributed to impedance contrast between debris and underlying rock mass. Further, the DB_Md point in both the models showed amplification at ∼1.0 and 2.0–2.5 Hz. The amplification at a lower frequency, i.e., ∼1.0 Hz, may be attributed to the thick (40–60 m) layer of debris that possibly decreases the resonance frequency and results in amplification of ground motion as also reported by Beresnev and Wen (1996). The amplification at 2.0–2.5 Hz may be linked to the elevated topography at these points in both the models.

The DB_Up point in both the models has different responses. In the CS-2 model, it showed amplification at 1.0–1.5 Hz, whereas in the CS-3 model, the spectral ratio is relatively stagnant except minor amplification at 4.0 and 6.0 Hz. This contrast may be understood by the fact that in CS-2, this monitoring point is situated at a thicker and elevated surface, whereas in CS-3, it is situated at relatively shallow topography and on top of relatively thin landslide thickness. Finally, the crown (CR) point also has a different spectral ratio in both the models. It shows higher amplification in the CS-3 model than the CS-2 model, which may be linked to the positioning of these points. The CR in the CS-2 is situated at a relatively flat surface unlike in the CS-3 model where it is situated at a steep surface. Thus, the monitoring points showed amplification at multiple frequency ranges that are attributed to the complex topography of landslides, soil thickness variation, and impedance contrast.

4.2 Landslide runout pattern

In view of uncertainties in ascertaining the exact depth of loose material that will be eroded/entrained during the debris flow, the runout pattern was evaluated at four different release area depths: 5, 10, 15, and 20 m of the loose overburden (Fig. 13a and b). The identification of the release area was based on field and satellite imagery observations. Four factors – gullies (Fig. 2c), flow relics (Fig. 2d), signs of failure (Fig. 2e), and overburden thickness pattern (Fig. 5c) – were considered while selecting the release area. The thickness region of 60–80 m having flow relics and signs of failure was therefore selected as the potential release area (Fig. 13b). Debris flow characteristics (flow height and flow velocity) of the debris flow that will strike the river floor during such an event are also inferred along the river channel (Fig. 13c). The debris flow height and velocity at the hillslope and along the river channel are summarized in Table 3.

Figure 13 Debris flow runout pattern. (a) Soil (or debris) thickness pattern in the landslide. (b) Different release area depths (5, 10, 15, and 20 m) used for the analysis. (c) River profile section A–B used to represent the resultant debris flow along the river. (d–f) Results at 5 m depth. (g–i) Results at 10 m depth. (j–l) Results at 15 m depth. (m–o) Results at 20 m depth.

Figure 14 Relationship of channel width and debris flow characteristics (a–c) and results of sensitivity analysis at constant “maximum” friction and variable turbulence (d–f). “Ups” and “Dws” refer to the upstream and downstream sides of the river channel, respectively.

Table 3 Results of the runout simulation at different depths of the release area.

As can be seen from Fig. 13, increasing depth of the release area increases the flow characteristics at the hillslope. However, the flow characteristics vary once the flow strikes the river floor. Debris flow height increases towards downstream part of the river channel on increasing the depth of material (Fig. 13f, i, l, and o). This behavior can be linked to the gully/channel on the hillslope near the downstream part (Fig. 2c) that possibly accelerated the flow. Debris flow velocity, however, decreases on increasing the depth of material. Higher flow velocity at lower material depth (Fig. 13f) can be understood using the “turbulence (or Chézy resistance) ξ” factor used in Eq. (3). The Chézy resistance is famous as “turbulent” friction (Voellmy, 1955) since its mathematical formulations are similar to the well-known turbulent Chézy equation (Chow, 1959). According to the Chézy equation, lower material thickness results in higher flow velocity. Further, apart from such turbulence effects, river channel morphology also affects the flow characteristics. As can be seen in Fig. 14a–c, a smaller channel width (narrow sections) generally accommodates higher flow velocity and height, of course with some nonlinearity as seen at 10 m depth (Fig. 14a). Further, increasing the material depth increases the flow height more than flow velocity (Fig. 14c).

By keeping the maximum friction (μ=0.4) constant, a sensitivity analysis was also performed by using different turbulence coefficients (Fig. 14d–f). It revealed that flow height and velocity increase on increasing the turbulence coefficient (implying increasing liquid content). It is of note that the central value, i.e., 250 m s⁻², showing a moderate response, is used in the main findings, as also mentioned in Table 2. Further, flow velocity is found to be more responsive at lower turbulence, whereas at higher turbulence, flow height dominates (Fig. 14f). Further, in order to understand the extent of runout along the river channel, runout results at the maximum release area depth were also laid over Google Earth imagery (Fig. 15a and b). A top view of the landslide with the runout is shown in inset c. The predicted runout is noted to extend across the river channel mainly at two locations, one near the left flank (Fig. 15d) and the other near the right flank (Fig. 15e). At both of these locations, the river channel attains sinuosity in a range of ∼1.30–1.32 (shown through channel length measurement). The river channel might owe this sinuosity to the paleo-landslide and/or fluvial deposit that extends the slope toe at these locations. Thus, the runout findings of the present study are noted to follow the same spatial extent as possibly followed by previous landslide events.

Figure 15 Debris flow runout pattern at 20 m depth. (a) Upstream view of the landslide from the right flank. (b) Runout pattern at 20 m depth. (c) Top view of the landslide highlighting two regions where runout reached across the river. (d) Runout pattern near the left flank extending across the river channel. (e) Runout pattern near the right flank extending across the river channel.

5 Summary

The present state of the slope reveals an instability condition through the factor of safety in a range of 1.09–1.32 and potential displacement near the landslide crown (Figs. 8 and 9). Such a displacement near the landslide crown may be related to the development of shear failure in slopes (Matsui and San, 1992; Kumar et al., 2018, 2021). The possibility of shear failure becomes more viable in the case of degradation of the shear strength of slope material and/or rupture planes. Notably, both the main landslide-triggering factors – rainfall and earthquakes – have been found to degrade the shear strength of slope material through percolation and shaking-induced particle movements, respectively (Cai and Ugai, 2004; Chang and Taboada, 2009). The GW, implying the rainfall-induced percolation effect, decreases the factor of safety and increases the material displacement. This effect is attributed to the joint hydraulic pressure against the overlying loose material that decreases the normal stress and hence the shear strength of overlying loose material (Mohr, 1914; Witherspoon et al., 1980). Similarly to the GW effect in static conditions, the combined response of the dynamic force and the GW resulted in an increase in the displacement (Fig. 8). This can be attributed to the fact that seismic shaking increases the hydraulic pressure in the joints that causes enhanced material displacement in the overlying loose material (Wang et al., 2010).

Further, the ground motion amplification also revealed the slope instability (or potential deformation). The maximum value of peak ground acceleration (PGA) is attained by the upper part of the debris surface (near the landslide crown) (Fig. 11). It is linked to the impedance contrast between underlying rock mass and overlying soil and/or partial loss of the shear strength during seismicity (Novak and Han, 1990; Šafak, 2001). Further, the spectral ratio also showed signal amplification, at multiple frequency ranges, at the debris surface (Fig. 12). Such an amplification at multiple frequency ranges is attributed to the complex topography of landslides, soil thickness variation, and impedance contrast (Sect. 4.1.4). Such high amplification at the slope surface has also been considered a main cause of slope failure in many other studies (Lenti and Martino, 2012; Del Gaudio et al., 2019).

As also stated in Sect. 3.3, hillslopes affected by the seismic shaking have also been prone to rainfall-induced failures, particularly in the form of debris flows. Further, the earthquake-induced shear strength degradation of slope material may also result in enhanced entrainment during a debris flow event (Liu et al., 2020). These debris flows might be initiated by either increased pore pressure (or GW-induced hydraulic pressure) or runoff involving entrainment (Godt and Coe, 2007). Though the GW effect is seen in the slope instability (Figs. 8 and 9), the potential response of the slope under excessive rainfall is explored through debris flow runout analysis (Figs. 13–15).

The debris flow runout predictions revealed a nonlinear increase in the debris flow height and velocity along the river channel on increasing the depth of the release area (Fig. 13). This nonlinearity is attributed to the variation in the river channel width (Fig. 14) and influx of debris flow material from the slope. Though the present study noted the influence of channel morphology on the debris flow characteristics, other studies have observed the changes in channel morphology caused by debris flows (Remaître et al., 2005; Simoni et al., 2020).

Thus, there seems to be a positive feedback process between channel morphology and debris flow that is further strengthened by the finding of the debris flow extent across the river channel (Fig. 15d and e). At both of these locations, the slope toe extends towards the ESE direction, resulting in higher channel sinuosity. These extended slope toes probably represent paleo-landslide and/or fluvial deposits. Signs of flow relics at the slope surface and failure at the slope toe at these locations (Fig. 2d and e) further support the possibility of paleo-landslide deposits. Thus, the predicted extent of potential debris flow is found to follow the trails possibly created by previous landslide flow and/or slide events. Aforementioned findings – temporally increasing rainfall, soil moisture, and surface runoff (Sect. 2.2) – and frequent debris flows/flash floods in this region (Micu et al., 2013; Grecu et al., 2017; Micu, 2019) pose an increasing risk of debris flow in the study area.

Finally, there are still some uncertainties in such predictive approaches that are as follows: (1) the inclusion of the subsurface discontinuity network, spatially varying groundwater surface, and material heterogeneity in the 3D model and (2) the inclusion of variable depth and phases in the runout modeling. Despite these possible uncertainties, which will be overcome in future work, such studies are required to minimize the risk and avert possible disasters.

6 Conclusions

By utilizing field-based data and numerical simulations of a massive (∼9.1 Mm²) Varlaam landslide in the SE Carpathians (Romania), the present study explored the potential response of this landslide in a seismic and rainfall regime.

The slope revealed the factor of safety in a range of 1.09–1.32 along with a static displacement of 0.4–4 m that increases up to 8–60 m under seismic load. The groundwater, implying the saturation, further decreased the slope stability owing to enhanced joint hydraulic pressure. Ground motion amplification, during seismic shaking, further revealed the potential instability of the slope with a peak ground acceleration (PGA) on the slope surface in a range of 0.65–1.88 g. Such amplification pertains to the complex topography of landslides, soil thickness variation, and impedance contrast.

Further, though the GW effect is seen in the slope instability, the potential response of the slope under excessive rainfall is also evaluated through debris flow runout analysis. The predicted debris flow revealed a nonlinear increase in the debris flow height (9.0–26.0 m) and velocity (2.1–3.0 m s⁻¹) along the river channel. This variation along the river channel is attributed to the river channel morphology and influx of debris flow material from the slope. Owing to the predictive nature of the present study, the concept may be applied to other terrains subjected to frequent landslides mostly triggered by extreme rainfall and earthquakes.