Instantaneous limit equilibrium back analyses of major rockslides triggered during the 2016–2017 central Italy seismic sequence

Abstract

Among the almost 1400 landslides triggered by the shocks of the 2016–2017 central Italy seismic sequence, only a limited number, all classifiable as rockslides, involved volumes larger than 1000 m${}^{\mathrm{3}}$. Four of these failures, including the three largest among the documented landslides, were described in terms of structural and geomechanical investigations in a previous paper. In this study, the estimated acceleration time histories at the rockslide sites were evaluated through a 2D simplified numerical model accounting for the attenuation phenomena and for the topographic effect of the rock cliffs from which the slide detached. Instantaneous stability analyses were carried out to obtain insights into the variability of the instantaneous margin of safety along the motion, over the entire spectrum of mechanisms that could be activated. Finally, some general suggestions on the pseudo-static verification method for 3D cases are proposed, which represent useful indications to hazard evaluation at local and regional scales.

1 Introduction

Most of the landslides that occurred in central Italy during the earthquakes of the last 100 years were rockfalls and rockslides (Martino et al., 2017; Esposito et al., 2000). Rockfalls were few cubic meters in volume, while rockslides involved volumes lower than 40 000 m${}^{\mathrm{3}}$. Even though large rockslides have been recorded during moderate-to-large magnitude seismic events (see e.g., Sepúlveda et al., 2016), the type and size of rock slope failures that occurred during the central Italy earthquakes are compatible with their level of released energy (moment magnitude, $M_{\mathrm{w}}\le \mathrm{6.5}$) and to the lithology of the formations outcropping in the areas close to the seismogenic structures. They involved relatively small rock volumes in comparison to the relationships reported in the literature (Malamud et al., 2004; Marc et al., 2016), although similar rock slope failures were recorded during the 1976 Friuli seismic sequence (Govi and Sorzana, 1977), which was characterized by comparable seismic input and surface geology. Similarly to many other earthquake-triggered landslides (Rodriguez et al., 1999), the landslides that occurred during the central Italy earthquakes were characterized by a marked disruption of the rock mass and originated on very steep slopes, where inertial forces easily remove well-delimited and scarcely constrained blocks from the slope through rigid sliding/toppling or tensile failures of overhanging blocks. All these instability phenomena are very “brittle”; i.e., relatively small displacements develop before the constraints change abruptly, and a fast propagation phase begins with a free fall of single or multiple non-interacting blocks (Esposito et al., 2000; Lanzo et al., 2009; Stewart et al., 2018; Franke et al., 2019).

Table 1

Main features of the investigated rockslides.

Landslide		Possible triggering earthquakes	Estimated	Lithology
		(date, GMT, moment magnitude)	volume
		(m${}^{\mathrm{3}}$ $\times$ 10${}^{\mathrm{3}})$
Nera (Sasso Pizzuto Mt.)	NR	30 October 2016, 06:40:18, $M_{\mathrm{w}}=\mathrm{6.5}$	32.0	Layered limestones (Maiolica Fm.)
Costa Cattiva (Nera River valley)	CC	26 October 2016, 17:10:36, $M_{\mathrm{w}}=\mathrm{5.4}$	0.4	Layered limestones
		26 October 2016, 19:18:06, $M_{\mathrm{w}}=\mathrm{5.9}$
Rubbiano (Infernaccio gorge)	RB	30 October 2016, 06:40:18, $M_{\mathrm{w}}=\mathrm{6.5}$	15.0	Layered limestones
Piè la Rocca (Patino Mt.)	PR	24 August 2016, 01:36:32, $M_{\mathrm{w}}=\mathrm{6.0}$	15.0	Massive limestones (Calcare Massiccio Fm.)
		24 August 2016, 02:33:29, $M_{\mathrm{w}}=\mathrm{5.3}$

Since seismic loading only acts at the early detachment of earthquake-triggered rock failures (propagation is controlled only by gravity loading and slope geometry), the study of this stage is very important for hazard evaluation. Literature has evidenced that the intrinsic key factors influencing the seismic activation of these landslides include the structural features of the rock mass (e.g., joint spacing and orientation, presence of major joints) and the hydraulic conditions (e.g., Massey et al., 2017; Tsou et al., 2018). Also the topographic modifications to ground motion and the effects of step-like slope topography on seismic motion have been studied by many authors such as Ashford et al. (1997), Bouckovalas and Papadimitriou (2005), Nguyen and Gatmiri (2007), Sepúlveda et al. (2005a, b), Pagliaroli and Lanzo (2008), Lenti and Martino (2012), Li et al. (2019), and Pignalosa et al. (2022). Li et al. (2019) described parametric analyses of steep rock slopes, providing the amplification factors for only weak seismic excitations and on the assumptions of homogeneous and elastic rock materials. Nevertheless, the general influence of very steep and vertical slopes, like those from which the rockslides detached, is scarcely investigated in detail in the literature, and it is a valuable topic to investigate. Also the amplification distribution along the slope profile, that is important to estimate the inertial effects on the surface rock blocks, is not enough evidenced, and only quite rough approximations are applied, such as a linear variation between the crest and the toe (Mavrouli et al., 2009).

The fundamental method used to study and estimate the block stability in a 3D approach is the limit equilibrium method for rock wedges (Goodman, 1976; Hoek and Bray, 1977). The method was extended to the dynamic conditions by Ling and Cheng (1997) through a pseudo-static approach that was also experimentally verified by Kumsar et al. (2000). However in these works the three-dimensionality of the seismic action is not usually taken into account, as the inertial forces are applied in the direction normal to the slope. This practice can disguise the possible activation of different mechanisms during shaking, because the real mechanism of initial failure is determined by the direction of the resultant external force on the block (Goodman, 1976).

These considerations sparked an investigation into the relationship of the rock mass structure (especially the orientation of major discontinuities) with the failure mechanism and in particular with the evolution during the seismic shaking for some of the largest rockslides among those reported during the 2016–2017 central Italy seismic sequence (CISS). Analyses were conducted with a 3D instantaneous limit equilibrium method (LEM) by applying the acceleration histories of the main earthquake shocks at the rockslide sites modified by a simple local viscoelastic seismic response for topographic conditions. Such an insight, when transferred to predictive stability analyses, can lead to a better awareness of the possible mechanisms and hence to a more effective evaluation of the hazard and of the successive fall/avalanche stage, which still represent a challenging problem (Wartman et al., 2019).

2 Examined rockslides and available data

The CISS consisted, from August 2016 to mid-January 2017, of several shocks ranging from $M_{\mathrm{w}}$ 5.0 to $M_{\mathrm{w}}$ 6.5 (the latter is the highest magnitude recorded in central Italy during the last 100 years; Rovida et al., 2019) and involving nearly 1500 km${}^{\mathrm{2}}$ of the regional normal fault system affecting an area characterized by a seismic gap between the 1997 $M_{\mathrm{w}}$ 6.1 Colfiorito–Sellano earthquakes to the north and the 2009 $M_{\mathrm{w}}$ 6.1 L'Aquila earthquake to the south). The seismic events caused more than 300 casualties, heavily damaging the physical environment, buildings and historical heritage as well (Miano et al., 2020; Saretta et al., 2021), and they triggered more than 1370 landslides, mainly rockfalls and slides, affecting limestone formations and, to a lesser extent, the flysch units. The shocks with $M_{\mathrm{w}}>\mathrm{5}$ triggered the largest rockslides, which mobilized volumes up to 35 000 m${}^{\mathrm{3}}$. Their features were described by Forte et al. (2021), Lanzo et al. (2019), Franke et al. (2019) and Romeo et al. (2017). After the seismic sequence, input data on the rockslides were collected during several investigation campaigns that included the following: aerophotogrammetric surveys with unmanned aerial vehicles (UAVs); sampling of blocks and joints; and direct in situ measurements of joint orientation, spacing, and roughness. Data and results of the reconnaissance investigations are reported in detail by Forte et al. (2021), who described the local geology, the rock mass structure, the major joints delimiting the failed mass and the failed volumes.

Main features of the selected four rockslides and their possible triggering earthquakes are shown in Table 1, while their location is reported in Fig. 1 together with the epicenters of the main shocks of the CISS. Table 2 summarizes the main features of the shocks that triggered the rockslides, obtained from signals recorded at the neighboring accelerometric stations on stiff ground (Engineering Strong-Motion Database, ESM; Luzi et al., 2016).

Figure 1

Epicenters of the 2016–2017 central Italy seismic sequence (CISS) and location of the studied landslides on a simplified geological map (modified after Forte et al., 2019). Keys: carbonate bedrock (CB), marly carbonate bedrock (MCB), arenaceous flysch bedrock (AFB), clayey flysch bedrock (CFB), conglomerate bedrock (CgB), travertine (tv), debris (db), terraced conglomerates (tcg), gravels and sands (gs), Gorzano–Laga fault (GLF), and Vettore–Bove fault (VBF).

[Figure omitted. See PDF]

These rockslides were chosen because they represent four of the largest failures among those detected during the reconnaissance field surveys conducted immediately after the seismic shocks (Costa Cattiva (CC) and Nera (NR) rockslides) or the most accessible among those observed on aerial images taken soon after the end of the seismic sequence (Piè la Rocca (PR) and Rubbiano (RB) rockslides). In this way, UAV surveys, which allowed for a detailed morphological and geostructural setting, could be conducted in a relatively short time after the seismic sequence. Other large rockslides detected on aerial images, with much higher logistic issues, were successively investigated and are currently being analyzed in the framework of national research projects.

They originated from the limestone formations known as Calcare Massiccio Fm. and Maiolica Fm. and were studied by merging classical field methods with newer remote sensing approaches by UAVs. A comparison of geostructural analyses at the scale of the slope with the regional tectonic setting indicated that all the four rockslides are locally characterized by major discontinuities of the older anti-Apennine sets (NE–SW) despite a much higher frequency of the quaternary Apennine tectonic sets (NW–SE). The failures often occurred following the breakage of rock bridges during the seismic shaking, as pointed out by stability analyses and evidenced on UAV images. Figure 2a through 2d present post-collapse frontal views of the rockslides. The bedrock is made of limestone for all the four cases, which is either layered (Costa Cattiva, Nera and Rubbiano rockslides) or relatively massive (Piè la Rocca rockslide). The four rock slopes are all very steep, and three of them (Nera, Piè la Rocca and Rubbiano) are located within tectonically disturbed zones (reverse fault and associated fold hinge, a fault zone, and a thrust front, respectively). The wedges were all delimited by near-planar single major joints (labeled in Fig. 2), excepting for the Rubbiano rockslide, which was delimited at its back by a surface resulting from the combination of several discontinuities of limited extent. Figure 2 also includes stereo plots with great circles of the planes delimiting each wedge at the very beginning of the detachment, as estimated from the 3D models and point clouds obtained from UAV aerial surveys (Franke et al., 2019; Tommasi et al., 2019). Great circles refer to single major joints or planes interpolating combinations of minor joints. Low-dip joints (i.e., along which shear occurred) showed negligible intact rock bridges excepting for that delimiting at the base of the Nera rockslide; its contribution to shear strength was therefore considered. Portions of intact rock were found along the subvertical surfaces delimiting the back of two of the failed wedges, where they provided some tensile resistance (Piè la Rocca and the Rubbiano rockslides). The latter was large enough to deserve consideration in the stability analyses.

Figure 2

Frontal view of the scars of the studied rockslides with approximated limits of the mean delimiting planes: (a) Nera (NR), (b) Costa Cattiva (CC), (c) Piè la Rocca (PR) and (d) Rubbiano (RB). In the inner boxes: stereographic projections (lower hemisphere) of the discontinuity planes delimiting the failed masses (1, 2, 3) and of the local slope face (f). Details are in Forte et al. (2021).

[Figure omitted. See PDF]

Observations of Google Earth satellite images taken at different dates over 2016 indicate that the Nera and Rubbiano rockslides occurred during the strongest shock of the seismic sequence (Norcia, 30 October, $M_{\mathrm{w}}=\mathrm{6.5}$). The Piè la Rocca slide occurred during one of the two 24 August shocks: Accumoli $M_{\mathrm{w}}=\mathrm{6.0}$ or Norcia $M_{\mathrm{w}}=\mathrm{5.3}$, the latter having an epicenter very close to the site (4.2 km). The Costa Cattiva rockslide occurred during one of the two 26 October shocks: either Visso $M_{\mathrm{w}}=\mathrm{5.4}$ or Castelsantangelo $M_{\mathrm{w}}=\mathrm{5.9}$. Locations of the epicenters are reported in Fig. 1.

Since the failed rock slopes were not accessible for geophysical investigation, the shear wave velocity, $V_{\mathrm{s}}$, used in seismic response assessment was estimated based on results of borehole geophysics conducted on the same geological formations at neighboring sites having similar fracturing and loosening of the rock mass. Down holes seismic tests in the Maiolica Formation conducted in the framework of seismic microzonation of the struck area (Banca dati microzonazione sismica, https://www.webms.it/, last access: 9 March 2023) indicated that pervasively fractured rock (i.e., with RQD values close to 0) exhibits a $V_{\mathrm{s}}$ of about 600 m s${}^{-\mathrm{1}}$, which increases to 2000 m s${}^{-\mathrm{1}}$ in a fairly jointed rock mass.

The strength parameters adopted for the stability analysis of the rockslides, which are the same used by Forte et al. (2021), were derived from direct in situ investigations at the rockslide sites and on laboratory tests on samples collected both at the rockslide sites and at neighboring sites in the same geological formations. The friction angle along the sliding planes varies between 40 and 47${}^{\circ }$ depending on the local roughness and waviness. Intact rock bridges along the joints forming the slide scar, which broke during failure, were observed on close UAV images. Their contribution to joint cohesion was evaluated as the cohesive component of the shear strength of the rock mass ($c_{\mathrm{rb}}=\mathrm{570}$ kPa) multiplied by the area $A_{\mathrm{rb}}$ of the rock bridges. The parameter $c_{\mathrm{rb}}$ was estimated by linearizing the Hoek–Brown strength envelope of the rock mass, obtained from the strength envelope of the rock material scaled through the Geological Strength Index, GSI, determined on the rock outcrops at the sites (Hoek et al., 2002).

High-resolution imagery captured from UAVs and during helicopter surveys over the Nera slide also revealed that the tip of the sliding surface appeared to be irregular and paler than the surrounding rock mass. This evidence induced Forte et al. (2021) to hypothesize that the lower part of the failure surface developed through the rock mass rather than along an existing joint. Therefore, an additional contribution was considered by multiplying the cohesion $c_{\mathrm{rb}}$ by the area (800 m${}^{\mathrm{2}})$ of the failure surface at the wedge tip.

For the Rubbiano (RB) rockslide, a tensile strength equal to 10 % of the rock mass cohesion $c_{\mathrm{rb}}$ was considered as an additional strength contribution that contrasted the detachment from plane 1 (Fig. 2d). Where plane 3 is present (as at the top of the Piè la Rocca rockslide, Fig. 2c) the wedge detaches along it from the rock mass behind, thus providing no strength contribution. Geometry and strength parameters adopted in the static LEM back analyses described by Forte et al. (2021) are shown in Table 3.

Table 2

Peak ground accelerations from the available records of the shocks at stations installed on rigid outcrop within 50 km from the epicenters of the seismic events.

			Seismic	Epicentral		Horizontal	Horizontal	Vertical
Event			station	distance	$D_{\mathrm{JB}}$	E–W PGA	N–S PGA	PGA
Epicenter	Date	$M_{\mathrm{w}}$		km	km	m s${}^{-\mathrm{2}}$	m s${}^{-\mathrm{2}}$	m s${}^{-\mathrm{2}}$
Accumoli	24 Aug 2016	6.0	MTR	19.40	11.40	0.791	0.754	0.418
			LSS	26.70	22.22	0.230	0.190	0.151
Norcia	24 Aug 2016	5.3	MTR	30.80	28.13	0.295	0.305	0.137
			ANT	42.00	39.68	0.118	0.112	0.048
Visso	26 Oct 2016	5.9	T1212	18.8	12.1	0.667	0.866	0.445
			LSS	41.1	33.9	0.128	0.110	0.118
			MTR	43.8	36.2	0.174	0.168	0.091
Castelsantangelo	26 Oct 2016	5.4	T1212	15.2	12.3	1.767	1.917	0.588
sul Nera.			RQT	17.4	3.7	2.177	1.296	0.880
			LSS	37.6	23.8	0.148	0.139	0.079
			MTR	40.9	22.6	0.179	0.266	0.085
Norcia	30 Oct 2016	6.5	T1212	10.50	8.77	2.744	2.731	1.636
			LSS	32.60	25.10	0.464	0.523	0.399
			ANT	46.10	33.27	0.436	0.546	0.242

Table 3

Input parameters for the static LEM stability analyses after Forte et al. (2021).

		Plane 1			Plane 2			Plane 3	LEM (static condition)
Rockslide	Volume (m${}^{\mathrm{3}}$)	dip/dd	${\mathit{\phi }}_{\mathrm{1}}^{\prime }$	$A_{\mathrm{rb}\mathrm{1}}$	dip/dd	${\mathit{\phi }}_{\mathrm{2}}^{\prime }$	$A_{\mathrm{rb}\mathrm{2}}$	dip/dd	Mechanism	$F_{\mathrm{S}}$
		${}^{\circ }$/${}^{\circ }$	${}^{\circ }$	m${}^{\mathrm{2}}$	${}^{\circ }$/${}^{\circ }$	${}^{\circ }$	m${}^{\mathrm{2}}$	${}^{\circ }$/${}^{\circ }$		–
NR	30 940	77/337	47	570$+$800	60/270	40	0	48/95	Line $I_{\mathrm{12}}$	1.68${}^{\mathrm{a}}$
CC	400	75/330	47	0	35/090	40	0	–	Line $I_{\mathrm{12}}$	2.16
PR	14 000	75/330	47	0	40/255	42	0	72/106	On plane 2	1.07
RB	15 000	65/084	47	0	85/130	40	2880	–	On plane 1	1.11${}^{\mathrm{b}}$

${}^{\mathrm{a}}$: with the cohesive contribution of the spur at the lower wedge tip (800 m${}^{\mathrm{2}})$. ${}^{\mathrm{b}}$: with the tensile contribution of the rear wedge surface (composite surface labeled as plane 2). $A_{\mathrm{rb}}$: areas of intact rock along the sliding planes providing cohesive contribution.

3 Method of analysis

For each rockslide, after a kinematic analysis in static conditions, the seismic motion responsible of failure was estimated in two steps: a ground motion prediction equation (GMPE) was applied first to the available ground motion records to account for attenuation at the rockslide sites, then modifications to ground motion induced by the slope morphology were evaluated through a general simplified 2D numerical model that reproduces the main resonance and attenuation phenomena affecting very steep rock slopes during seismic shocks.

Since the investigated slopes can be roughly assimilated to steep flanks of deep valleys (200–500 m) separated by large and relatively flat mountain ridges, the modifications to the ground motion that the general slope morphology produces at each site were estimated through a finite difference model that simulates the viscoelastic dynamic behavior of a simplified slope: a step-like slope with a vertical cliff of height $H$ and upstream and downstream horizontal areas (Fig. 3).

Figure 3

Numerical model used for the seismic response of a vertical cliff.

[Figure omitted. See PDF]

In the frequency domain, the modification of a harmonic motion of wavelength $\mathit{\lambda }$ (frequency $f)$ propagating in a medium with shear wave velocity $V_{\mathrm{s}}$ and Poisson ratio $v$ can be expressed by the amplification ratio $A$ between the amplitude at a point at height $h$ on the slope face and the amplitude on the horizontal rigid outcrop. In a dimensional analysis approach, $A$ can be expressed as a function of the dimensionless variables $\mathit{\zeta }$, $\mathit{\eta }$ and $\mathit{\nu }$, being $\mathit{\zeta }=h/H$ and $\mathit{\eta }=H/\mathit{\lambda }=Hf/V_{\mathrm{s}}$.

If planar vertically propagating waves that are polarized are considered, two cases (P- and S-waves) are sufficient to estimate the variability of $A$ along the entire vertical wall ($\mathit{\zeta }$ in the range 0.0–1.0). To cover a sufficiently broad frequency range ($\mathit{\eta }=\mathrm{0.03}$–2.0), two analyses were performed for each wave type using Ricker waves with different frequency content as input motion.

The model, built in the 2D finite difference code FLAC (Itasca, 2011), consists of a 100 m high slope. The model bottom is an absorbing viscous boundary, while free-field boundary conditions are applied to the lateral boundaries, which are located at least 5 times $H$ from the slope face (Fig. 3). A Rayleigh formulation was assumed with a uniform critical damping ratio of $D=\mathrm{0.5}$ %. Elastic properties $V_{\mathrm{s}}=\mathrm{100}$ m s${}^{-\mathrm{1}}$ and $\mathit{\nu }=\mathrm{0.3}$ were considered, but the normalized results can be extended to a cliff of different height and stiffness thanks to the principle of linear superposition.

Finally, the behavior of the rock wedges during the shocks was analyzed using the instantaneous limit equilibrium procedure in the instrumental hypothesis that they moved rigidly with the surrounding rock mass. In fact, the volumes of the rockslides are in fact small enough to consider the primary mechanism of failure as a rigid wedge (Hungr et al., 2014) and to explore the role of the inertial forces in inducing failure. At this early stage, very small displacements occurred mainly as sliding, and the constraint configuration was dictated by the original orientation of preexisting joints. The possible subsequent kinematic evolution (e.g., toppling), the disarrangement of the wedge and the start of the propagation phase are out of the scope of this paper. The estimated seismic input was applied to calculate the time histories of the safety factors during the shocks for the rockslides under the same assumptions of static LEM back analyses conducted by Forte et al. (2021): the landslide body is subjected to gravity only, water pressure is absent, sliding surfaces are planar and a Mohr–Coulomb strength criterion is assumed. Only translational sliding mechanisms were considered, because their predominant role in driving the wedges to the collapse clearly emerged (Forte et al., 2021).

Variable inertial forces $\mathit{I}\left(t\right)=-m\mathit{A}\left(t\right)$ applied uniformly to the rigid blocks with mass $m$ were added to equilibrium equations. This procedure is mechanically consistent only as long as the block does not displace with respect to the rock mass; therefore it only provides a realistic assessment immediately before sliding begins. In fact, the relative motion alters the inertial forces with respect to those calculated with the base acceleration and furthermore reduces the strength due to progressive smoothing of the joint surface and failure of the rock bridges. For these reasons the safety factor ($F_{\mathrm{S}})$ during seismic excitation calculated through this analysis is intended to identify the most probable instants of failure initiation and the critical mechanisms. This type of calculation also represents an instrument to weight the relative importance of sliding mechanisms during shaking and thus to better handle the pseudo-static analysis method. Safety factors are calculated for the whole shaking duration, and in turn it could also assume values lower than 1.0 during some time intervals.

In resolving equilibrium and calculating $F_{\mathrm{S}}$, the activation of a different translational mechanism with respect to that occurring in static conditions was also considered. In fact, the instantaneous sliding mechanism is controlled by the current direction of the resultant external force, which in a dry slope coincides with the sum of the block weight and the inertial force. The number of passages between different mechanisms during the seismic event is related to the oscillating amplitude of the resultant force and to the distance of its pole from the kinematical region boundaries.

The instantaneous $F_{\mathrm{S}}$ for the Mohr–Coulomb strength criterion can be calculated through Eqs. () and () in the case of sliding along a single plane or both planes, respectively: 1$$F_{\mathrm{S}}={\displaystyle \frac{cA+Ntg\mathit{\phi }}{T}},$$ 2$$F_{\mathrm{S}}={\displaystyle \frac{c_{\mathrm{1}}{A}_{\mathrm{1}}+N_{\mathrm{1}}tg{\mathit{\phi }}_{\mathrm{1}}+c_{\mathrm{2}}{A}_{\mathrm{2}}+N_{\mathrm{2}}tg{\mathit{\phi }}_{\mathrm{2}}}{T_{\mathrm{12}}}},$$ where $N$ and $T$ in Eq.  () are the normal and tangential components of the resultant acting on a single sliding plane. $T_{\mathrm{12}}$, $N_{\mathrm{1}}$ and $N_{\mathrm{2}}$ in Eq. () are the components of the resultant force parallel to the intersection line $I_{\mathrm{12}}$ of planes 1 and 2 and normal to the planes, respectively. $c$, $\mathit{\phi }$ and $A$ are the cohesion, friction angle and the contact areas, respectively; subscripts refer to the plane. The passage from one mechanism to another entails an instantaneous change in the $F_{\mathrm{S}}$ value.

4 Seismic input at the rockslide sites

The estimate of peak ground acceleration, PGA, calibrated with an appropriate GMPE at each of the four rockslide sites, is presented in Fig. 4. Diamonds represent the measured PGAs of the horizontal components (i.e., geometric mean of east and north components) versus the Joyner–Boore distance ($D_{\mathrm{JB}})$ of the station. These measurements were interpolated with and calibrated against the GMPE of Bindi et al. (2011) (solid blue lines in Fig. 4). The site class A (rigid ground, according to Eurocode 8, EN 1998-5:2004) and the normal fault class were used in the GMPE, while the moment magnitude was used as the regression parameter. The PGAs at the sites were finally estimated using the $D_{\mathrm{JB}}$ of each site on the interpolated GMPE (empty blue circle symbols in Fig. 4). The input accelerograms were obtained by linearly scaling the recordings at the closest station on rock outcrop to the estimated PGAs; the same scaling factors, $S$, were used for all the components. The parameters used in both the GMPE calibration and the scaling procedure are shown in Table 4.

Figure 4

PGA estimate through a calibrated GMPE for each shock. For the PR rockslide and the CC rockslide, two shocks are considered on 24 August and on 26 October, respectively.

[Figure omitted. See PDF]

Table 4

Parameters of the triggering events utilized to calculate the motion at the rockslide sites from the available recorded accelerograms.

Rockslide
						Epicentral
Site	Lat.	Long.	Seismic event			distance	$D_{\mathrm{JB}}$	$S$	PGA ${}^{\ast }$
	${}^{\circ }$	${}^{\circ }$	Epicenter	Date	$M_{\mathrm{w}}$	km	km	–	g
NR	42.93	13.07	Norcia	30 Oct 2016	6.5	10.2	2.1	1.103	0.309
RB	42.93	13.28	Norcia	30 Oct 2016	6.5	16.8	6.7	0.926	0.259
CC	42.92	13.12	Visso	26 Oct 2016	5.9	2.3	0.0	1.862	0.146
			Castelsant.	26 Oct 2016	5.4	4.4	0.0	1.186	0.223
PR	42.82	13.13	Accumoli	24 Aug 2016	6.0	4.6	1.0	0.668	0.054
			Norcia	24 Aug 2016	5.3	16.3	14.8	1.890	0.154

${}^{\ast }$ Peak ground acceleration estimated at the site on rigid horizontal outcrop.

For the estimate of the amplification effects at the four sites the plots derived from the finite difference model described in section 2 were utilized (Fig. 5). The plots report amplification ratios of the normal ($A_{\mathrm{n}}^{\mathrm{P}},A_{\mathrm{n}}^{\mathrm{S}})$ and vertical ($A_{\mathrm{z}}^{\mathrm{P}}$, $A_{\mathrm{z}}^{\mathrm{S}})$ component, for the incoming P- and S-waves, and are functions of $\mathit{\zeta }$ and $\mathit{\eta }$. Along a horizontal line (i.e., at constant ${\mathit{\zeta }}_{\mathrm{0}})$, the diagrams give the amplitude of the transfer functions from the outcrop motion (horizontal and vertical component, respectively, for the incoming S- and P-waves) to the motion of a point at height $h={\mathit{\zeta }}_{\mathrm{0}}H$ on the cliff.

Figure 5

Amplification ratios with respect to the outcrop motion along the vertical wall of a step-like slope. Incident S-wave ($A_{\mathrm{n}}^{\mathrm{S}}$, $A_{\mathrm{z}}^{\mathrm{S}})$ and incident P-wave ($A_{\mathrm{n}}^{\mathrm{P}}$, $A_{\mathrm{z}}^{\mathrm{P}})$.

[Figure omitted. See PDF]

According to literature results (Ashford et al., 1997; Assimaki et al., 2005), for the incident S-waves the most amplified wavelength corresponds to the first normalized modal frequency of 0.2, with a peak of $A_{\mathrm{n}}^{\mathrm{S}}$ greater than 1.4. For the Nera rockslide ($H=\mathrm{400}$ m), the main resonance frequency, which is about 1.0 Hz, is critical along almost the entire cliff, although amplification decreases as elevation decreases. In addition, at medium and lower elevations, the higher frequencies are reduced overall. The vertical component produced by incident S-waves has significant amplitude ratios $A_{\mathrm{z}}^{\mathrm{S}}$ only at the crest and for normalized frequencies in the wide range 0.4–1.4 (e.g., 2.0–7.0 Hz for the Nera case with about $V_{\mathrm{S}}=\mathrm{2}$ km s${}^{-\mathrm{1}}$).

The amplification ratios for incident P-waves, $A_{\mathrm{n}}^{\mathrm{P}}$ and $A_{\mathrm{z}}^{\mathrm{P}}$, reveal a main amplification of the vertical component at the crest and at almost the whole vertical wall for a normalized frequency of about 0.1 (0.5 Hz for the Nera case). Conversely, the horizontal component is flattened all along the cliff wall for all frequencies ($A_{\mathrm{n}}^{\mathrm{P}}<\mathrm{0.8}$).

The linear process used to assess the motion at the elevation of the center of gravity of the rockslide proceeds as follows. Since the normal (horizontal) and vertical components $a_{\mathrm{n}}$ and $a_{\mathrm{z}}$ of the outcrop acceleration can be considered equivalent, respectively, to an S- and a P-wave in a vertical plane normal to the slope face, their Fourier transform $a_{\mathrm{n}}\left(\mathit{\eta }\right)$ and $a_{\mathrm{z}}\left(\mathit{\eta }\right)$ are multiplied to the transfer (amplification) functions to obtain the output components on the cliff (i.e., after morphological modifications). These are successively combined: $$\begin{array}{}\text{3a}& {\displaystyle }{a}_{\mathrm{out},\mathrm{n}}\left(\mathit{\eta }\right)=A_{\mathrm{n}}^{\mathrm{S}}(\mathit{\eta },{\mathit{\zeta }}_{\mathrm{0}})a_{\mathrm{n}}\left(\mathit{\eta }\right)+A_{\mathrm{n}}^{\mathrm{P}}(\mathit{\eta },{\mathit{\zeta }}_{\mathrm{0}})a_{\mathrm{z}}\left(\mathit{\eta }\right),\text{3b}& {\displaystyle }{a}_{\mathrm{out},z}\left(\mathit{\eta }\right)=A_{\mathrm{z}}^{\mathrm{S}}(\mathit{\eta },{\mathit{\zeta }}_{\mathrm{0}})a_{\mathrm{n}}\left(\mathit{\eta }\right)+A_{\mathrm{z}}^{\mathrm{P}}(\mathit{\eta },{\mathit{\zeta }}_{\mathrm{0}})a_{\mathrm{z}}\left(\mathit{\eta }\right),\end{array}$$ where ${\mathit{\zeta }}_{\mathrm{0}}=h_{\mathrm{0}}/H$ is the normalized height of the rockslide gravity center. Finally the acceleration vector $\mathit{A}({\mathit{\zeta }}_{\mathrm{0}}$, $t$) is obtained by applying the inverse Fourier transform to the (3) and assuming a shear wave velocity of 2200 m s${}^{-\mathrm{1}}$ at all sites. The hypothesis that the seismic response develops in plane strain conditions can be assumed for a very long cliff or valley, and therefore the component $a_{\mathrm{p}}\left(t\right)$ can be considered unmodified. The geometrical and morphological features that control the response calculations are shown in Table 5 for each rockslide.

Table 5

Parameters utilized to calculate the topographic modification of the seismic motion at the rockslide sites.

			Height of	Period
	Dip direction	Cliff	the rockslide	of first
	of the slope	height,	center of	mode,
Rockslide	face, $\mathit{\alpha }$ (${}^{\circ }$)	$H$ (m)	mass, $h$ (m)	$T_{\mathrm{0}}$ (s)
Nera (NR)	330	400	250	0.91
Costa Cattiva (CC)	330	300	90	0.68
Rubbiano (RB)	115	530	170	1.20
Piè la Rocca (PR)	330	300	180	0.68

In Fig. 6 the amplification/attenuation effects are described by comparing the acceleration response spectra (damping 5 %) of the two motions for all the considered shocks and for both the horizontal (normal) and vertical components. The modification of the motion is usually significant for periods lower than 1–2 s, while it is negligible for periods higher than the fundamental period $T_{\mathrm{0}}$ of the first vibration mode of the cliff ($\mathit{\lambda }=\mathrm{5}H)$ indicated by vertical dotted lines.

These analyses aim to estimate the general modifications to the seismic motion caused by large-scale morphological features. Nonetheless, local irregularity of the slope surface like sharp ridges, spurs and pinnacles, which can induce significant further local amplifications/attenuations especially for small volumes, is not considered in the present research.

Figure 6

Response spectra of the three components of the acceleration at the rigid horizontal outcrop and on the vertical rock cliff at the elevation of the rockslide center of mass as estimated through the seismic response of the numerical model for the four case studies (for each possible triggering earthquake).

[Figure omitted. See PDF]

5 Seismic back analysis

The visual conception of the possible mechanism switches calculated during the instantaneous LEM analysis is represented by the kinematical regions reported in Fig. 7 for each rockslide. The regions are spherical triangles identified by the directions of the normal vectors to the planes and the directions of the plane intersection lines (Londe et al., 1969). Therefore the calculated (instantaneous) direction of the resultant external force on the block defines a different sliding mechanism depending on which triangle it belongs to. The directions of the initial static resultants (block weight) are indicated through red circles, but during the seismic analyses they move around and can cross over the triangle boundaries.

Figure 7

Stereographic conform projections (lower hemisphere projected from upper focal point) of the trihedral/dihedral regions (solid lines) that define different sliding mechanisms depending on the direction of the resultant force. $I_{ij}$, sliding along the intersection line between the planes $i$ and $j$ and pi, sliding along the plane $i$. Light dashed lines are the projections of the average local slope face. Red stars indicate the resultant orientations corresponding to the minimum $F_{\mathrm{s}}$ during the seismic shocks (see Fig. 8). Full triangles and squares indicate the orientation of the intersection line ($I_{\mathrm{12}})$ between the two planes and the dip direction of the slope face, respectively.

[Figure omitted. See PDF]

For each analysis, the instantaneous activated mechanisms and the time histories of $F_{\mathrm{S}}$ are reported in Fig. 8. For instance, the mechanism of the Nera rockslide (Fig. 8a) changes from sliding along the $i_{\mathrm{12}}$ to sliding along the single plane 2. The latter mechanism has quite a lower level of safety, and the $F_{\mathrm{S}}$ repeatedly crosses the critical threshold $F_{\mathrm{S}}=\mathrm{1}$ during the shock. This means that the available strength was reached since the very first oscillations and irreversible displacements grew up towards collapse. In cases similar to the Nera slide, a high $F_{\mathrm{S}}$ calculated in static conditions is not meaningful in evaluating the “distance” from failure in seismic conditions, as a very small deviation of the resultant force from the vertical direction can be sufficient to activate a less safe sliding mechanism also for moderate shaking.

Figure 8

Time histories of the safety factor $F_{\mathrm{s}}$ during the triggering shocks and of the instantaneous active mechanism: $I_{\mathrm{12}}$, sliding along the intersection line between the planes and $p_{\mathrm{1}}$ and $p_{\mathrm{2}}$, sliding along plane 1 and 2, respectively. Stars indicate minimum $F_{\mathrm{s}}$. Empty circles highlight instants with peak values of specific acceleration components: $x,y$, geographic components (E, N); $z$, vertical component; $h$, horizontal component; $n$, component normal to the slope face; i12, intersection line of the sliding planes; p1 and p2, dip directions of the sliding planes; and $m$, instant of maximum acceleration magnitude.

[Figure omitted. See PDF]

Due to the geostructural setting and the amplitude of the examined seismic shocks, the analysis of the Costa Cattiva rockslide yields a sliding mechanism along the line $i_{\mathrm{12}}$ both in static and dynamic conditions (Fig. 8b). The position of the blocks produced by the rock avalanche that followed the wedge failure confirms this mechanism. For both the 26 October shocks, the analyses do not justify the Costa Cattiva failure in seismic conditions (i.e., the computed $F_{\mathrm{S}}$ is $>\mathrm{1}$). Since an overestimation of the strength is improbable due to the simple structural conditions and the low joint roughness, this result is likely explained by having neglected the small-scale amplification. The wedge was in fact located on top of a narrow sharp ridge protruding from the slope.

The analysis of the Piè la Rocca slide (Fig. 8c) also helped to assess that the wedge likely failed during the Norcia event ($M_{\mathrm{w}}=\mathrm{5.3}$). The earlier event (Accumoli, $M_{\mathrm{w}}=\mathrm{6.0}$) $F_{\mathrm{S}}$ trespasses the stability threshold only once and for a very short time span, which could cause only very small displacements without reaching full collapse. Although the displacements experienced by the rockslide and the maximum available displacement before the collapse were not estimated, the geometric layout of the rockslide scar suggests that the wedge should have experienced displacements as large as to break, at least partially, a constraining rock spur at its highest part, whose failure surface is however small (3 %) compared to the area of plane 1. At the Piè la Rocca rockslide, the frequent switches between the two sliding modes determine a change in $F_{\mathrm{S}}$ that is not as important as for the Nera slide, because the two mechanisms have quite similar safety margins against failure.

The Rubbiano rockslide (Fig. 8d) maintains a unique mechanism during the application of the acceleration history, and despite the significant epicentral distance (7 km), the available strength was exceeded several times during the strong $M_{\mathrm{w}}$-6.5 shock. The structural layout of the slide scar and the slenderness of the detached wedge (small thickness normal to the cliff in comparison to the large extent parallel to the cliff) indicate that very small displacements were sufficient to reach the collapse, likely favored by a rocking effect.

Along the time histories shown in Fig. 8, the values of $F_{\mathrm{S}}$ are highlighted for particular instants: when the components along the geographical directions (E–W, N–S, up–down) and the horizontal component of the acceleration reach their maximum absolute values (respectively x, y, z, h points in Fig. 8) and when the acceleration vector magnitude reaches its maximum value ($m$ points in Fig. 8). It is apparent that, despite the fact that at these instants the inertial force is quite high, $F_{\mathrm{S}}$ is not always near its minimum. For some cases (e.g., $x,h,m$ conditions of the Costa Cattiva slide and the Piè la Rocca slide during the 26 October shock) these instants even correspond to the maximum values of $F_{\mathrm{S}}$.

Different instantaneous $F_{\mathrm{S}}$ values are also highlighted when the maximum values of particular acceleration components (all towards the slope and related to the geometry and orientation of the slope) are reached. These are the dip direction of the rock face ($n$ points in Fig. 8), the intersection line between the two main sliding planes (i12 points in Fig. 8) and the dip directions of the two planes (p1 and p2 points in Fig. 8). The results show that both $n$ and i12 conditions give the minimum $F_{\mathrm{S}}$ or a value near to the minimum of the entire shock. The only exception is represented by the i12 condition for the Nera slide. In this case, the resultant force falls in the region corresponding to the safer of the two possible mechanisms, and the i12 condition gives $F_{\mathrm{S}}=\mathrm{1.5}$, i.e., much higher than the minimum value ($F_{\mathrm{S}}=\mathrm{0.83}$).

These observations provide some clues for a rational choice of the direction of the inertial force to be applied on a 3D rock wedge in pseudo-static stability analyses. Due to the significant anisotropy of the mechanical problem, the inertia calculated at instants when the magnitude of the acceleration vector or that of some pre-defined components are maximum can have a negligible or favorable influence on stability. Conversely, the resistance to sliding can be overcome when the component along more adverse directions, either that along p1/p2 or i12, is significant. The orientation of the resultant forces determining the minimum $F_{\mathrm{s}}$ (stars in Fig. 7) always falls within sectors delimited by these two directions (thick dashed lines in Fig. 7). Application of the pseudo-static inertial force along these two directions yields the most conservative result only on the condition that they involve all the possible sliding mechanisms. Otherwise, if both these directions correspond to the same mechanism, other orientations of the inertial force should also be tested to verify the activation of different mechanisms.

6 Discussion

Stability analyses in static and seismic conditions were performed on four rockslides that occurred during the main shocks of the 2016–2017 central Italy seismic sequence. The failed masses can be realistically sketched as wedges delimited by two intersecting planar discontinuities and possibly by a detachment surface at their back. The activated primary mechanisms were sliding along either one plane or the intersection line between two planes. These mechanisms developed until the wedges lost their constraints and rockfalls/avalanches started. The volume of the rockslides (not exceeding 32 000 m${}^{\mathrm{3}})$ is small enough to assume an initial rigid motion of the wedges.

The available ground motion measurements were interpolated with an attenuation law with fixed source mechanism and stiffness class. Then a simple viscoelastic model was implemented in a parametric finite difference stress–strain analysis to calculate motion modifications due to the morphologic conditions, i.e., a step-like rock slope. Both the normalized results and the applications to the actual rockslide sites show that significant horizontal amplification is expected almost only at the crest, while at intermediate heights the main effect is a reduction of the horizontal component and an amplification of the vertical one. Cyclic strength degradation is another important issue that seems to have played an important role in most of the major rockslides described in the previous sections. The high number of loading cycles applied during the main earthquakes seems to have especially affected the rock bridges along persistent joints of the limestone formations, both under shear and in tension. In this respect, static limit equilibrium back analyses of the Nera rockslide indicate that rock bridges were necessary to ensure stability even in static conditions and also provided sufficient strength to maintain the wedge stability during the 26 October $M_{\mathrm{w}}\mathrm{6.0}$ event. Shear strength was most likely overcome during the successive 30 October $M_{\mathrm{w}}$ 6.5 shock.

The examined rockslides, which represent four of the largest landslides that occurred during the 2016–2017 sequence, are all characterized by a highly asymmetric wedge shape. This entails a low factor of safety due to the reduced (even null) strength contribution along one of the wedge planes. Even the static LEM stability analyses showed that the potential mechanisms are often not univocally established. In fact, either minor modifications of the geometrical layout or a small deviation of the resultant of external forces can activate a different mechanism with respect to that initially hypothesized. A clear representation of this problem can be obtained in the stereographic projection by subdividing the direction space into regions associated with different mechanisms. The variability of the mechanism is particularly significant in seismic conditions when the inertial force is added to the weight of the blocks thus making the resultant force fluctuate around its initial orientation.

The instantaneous LEM back analyses, carried out under the hypothesis that blocks are rigidly connected to the underlying bedrock, showed that also the safety margin can deeply fluctuate during the shock as a function of the mechanisms that are potentially activated. The minimum safety factor during the shock does not necessarily coincide with the typical directions of the pseudo-static force in a classic pseudo-static analysis (normal to the slope face or along the line of intersection between the sliding planes). Therefore, direction is to be varied through a rational and complete examination of all the possible mechanisms.

7 Conclusions

From the analysis of the local seismic response and of the stability in static and seismic conditions of the four investigated rockslides, the following considerations can be highlighted. All the rockslides experienced significant horizontal amplification only at the crest, while at intermediate heights the main effect is a reduction of the horizontal component and an amplification of the vertical one.

Instantaneous LEM stability analyses showed that the potential mechanisms are often not univocally established, as either minor modifications of the geometrical layout or a small deviation of the resultant of external forces can activate a different mechanism, and the safety margin can deeply fluctuate during the shock.

An increased awareness on the actual failure mechanism from instantaneous stability analyses can better orient more complex analyses both for the triggering stage and the propagation phase.

Cyclic strength degradation, due to the high number of loading cycles, likely played an important role in weakening both shear and tensile strength of the rock bridges along persistent joints of the limestone formations, at least at two sites.

The study indicates that the actual failure mechanism largely depends on the knowledge of the specific structural features of the slope. This issue is to be accounted for in risk analysis not only at local scales but also for long stretches of valley flanks overlooking transportation infrastructures in mountainous regions. In this respect, the extensive application of UAV surveys gave the possibility to obtain quantitative data even of slopes inaccessible to remote terrestrial surveys and demonstrated to be applicable at affordable times and costs along infrastructure stretches of considerable length. Quantitative data consist not only of the determination of geometry and structural setting of the slope but also of geomechanical parameters as medium- to large-scale roughness and extent of the rock bridges along major joints.

Code availability

The software used for this paper (FLAC) is of commercial use and can be acquired from the producing company ITASCA.

Data availability

Data can be accessed upon request to the authors.

Author contributions

LV: conceptualization, formal analysis and software, writing – original draft, writing – review, visualization; GF and MDF: geological investigation, writing – original draft, writing – review, visualization; PT: conceptualization, geotechnical investigation, funding acquisition, writing – review; GL: supervision, writing – review, funding acquisition; KF: UAV investigation, writing – review, funding; AS: investigation, supervision, funding acquisition.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Earthquake-induced hazards: ground motion amplification and ground failures”. It is not associated with a conference.

Acknowledgements

Valentina Tuccio carried out preliminary seismic response analyses for the Nera River valley rockslide. The authority for the Sibillini Mts. National Park (Alessandro Rossetti) and the Dipartimento della Protezione Civile (Paola Pagliara and Paola Bertuccioli) are acknowledged for allowing surveys in the park area and for access to landslide areas during the seismic emergency, respectively. BYU undergraduate students Bryce Berrett, Nicole Hastings, Jeffery Derricott, Bridgette Ostrum and Doug Graff, as well as graduate student Michael Freeman, contributed significantly to the UAV fieldwork and development of the 3D reconstructions of the rockfalls. The authors thank the two anonymous reviewers for their helpful suggestions.

Financial support

Melania De Falco, Giovanni Forte, and Antonio Santo were partially funded by the Italian Civil Protection Department's RELUIS project (2018), PR8-UR18-WP2 UNINA (Principal Investigator Antonio Santo); Giuseppe Lanzo, Paolo Tommasi, and Luca Verrucci by the projects Progetto di Ateneo Sapienza 2017 “Site investigations, monitoring and modelling of earthquake induced rockslides triggered by the 2016–2017 Central Italy seismic sequence” and Progetto di Ateneo Sapienza 2020 “Failures occurred during the 2016–2017 Central Italy seismic sequence based on site investigation and numerical modelling: hints for hazard mitigation in high seismicity areas” (Principal Investigator Giuseppe Lanzo); Paolo Tommasi by Progetto Fra.Si. funded by the Italian Ministry for Ecological Transition (Principal Investigator  P. Reichenbach); and Kevin W. Franke by a Mentored Research Grant (MRG) from Brigham Young University's Ira A. Fulton College of Engineering and the Center for Unmanned Aircraft Systems (C-UAS) and the National Science Foundation Industry/University Cooperative Research Center (I/UCRC) under NSF award no. IIP-1650547, along with significant contributions from C-UAS industry members.

Review statement

This paper was edited by Céline Bourdeau and reviewed by two anonymous referees.