1 Introduction

Landslides are important natural phenomena that reshape mountain and hill areas in subaerial and subaqueous zones. These processes can occur in several forms depending on numerous aspects such as the environmental characteristics, the trigger mechanism, the properties of the sliding material. Several kinds of classifications have been proposed. The most widely used is the one due to Varnes (1978), that, basically, categorizes the slides according to the material type (rock, debris, earth) and to the motion characteristics. In this work, we will focus our attention on rocks, without any constraints about the kind of motion. Commonly, rocks motion is simulated by means of block-based approaches where the rock slides are portioned in a set of blocks that can interact with each other according to various laws (e.g. Hungr 1995; Grilli and Watts 1999; Tinti et al. 1999, 2006; Lucas et al. 2011; Efremidis et al. 2015; Stamatopoulos and Di 2015).

In this paper, we introduce a new approach for modeling the dynamics of sliding bodies. We describe the motion of a mechanical system formed by \(N > 2\) (with \(N \in {\mathbb{N}}\)) mutually interacting point masses forming a 2D grid over analytical and real surfaces. Material point masses have been used in the description of a wide number of physical mechanisms, including modeling of explosions (Hu and Chen 2006) and landslides (Andersen and Andersen 2010; Liano-Serra et al. 2016). These applications make use of extensions of the so-called Material Point Method (Sulsky et al. 1994), employing a mutual Lagrangian and Eulerian description in which the landslide is divided into several discrete material points at which specific physical properties are evaluated. Particle conceptualization has been used also in the so-called Particle Finite Element Method where a combination of the standard finite element analysis and the particle approach is adopted to solve continuous problems like in fluid mechanics (Idelsohn et al. 2004) and, more recently, even in the analysis of the landslides processes (see Wang et al. 2019 for a synthetic description of the method).

The adoption we propose here to describe the motion of rocks represents a novelty, since it merges the traditional block approach with the material point masses technique, resulting in an easy tool for the general treatment of landslide processes. We make use of the geometry of the points grid to evaluate the internal forces acting on the system and we solve the equations of motion consequentially. In this way, we can estimate the inter-particles forces depending on the structure of the grid. In the above-mentioned particle techniques, this is not possible. The similarity with the block approach relies on the connection between the grid geometry and the general dynamics. This is often used to evaluate the shape or volume of the moving blocks (Tinti et al. 1997).

The main goals of this work are: (1) to introduce a new approach in the description of sliding bodies; (2) to prove that the basic principles of the particle-based techniques and the block-based approach can be synergistically applied resulting in a semi-rigid body description; (3) to validate the model through the simulation of the 1783 Scilla landslide, taken as a benchmark case, since it has been already deeply investigated in the literature.

The scheme of this paper is as follows. In the next section, we introduce the basic assumptions and governing equations. Section 3 serves to illustrate the method through simple analytical examples. Section 4 treats the landslide that occurred close to Scilla on 6th February 1783 during the 1783–84 Calabria seismic crisis. The landslide was triggered by an earthquake and caused a catastrophic tsunami, killing more than 1500 people. The discussion of the results in Sect. 5 will close the paper.

2 Formulation of the Problem

We assume that a landslide can be represented as a set of \(N\) point masses, whose total mass equals the mass of the landslide body and that slide on a surface. The point masses are the nodes of a grid and interactions take place only between the pairs of masses connected by grid sides, that are immaterial and are introduced only to indicate the pairs of interacting masses. Typically, each mass interacts only with neighbours and this local form of interaction can be accounted for by grids formed by triangles or quadrilaterals. This, however, is not a limitation for the model that can deal with any kind of grids. The forces acting on the point masses are gravity, which is the driving force, the reaction force exerted by the sliding surface, the basal friction and the drag force due to the action of the ambient fluid, the latter being practically relevant only when the masses move underwater. In addition, each mass interacts with the neighbour masses through forces that can be seen as internal forces for the \(N\)-mass system. As mentioned before, these forces are assumed to act along the sides of the grid. In a sense, each grid side is associated with a pair of masses and to a corresponding pair of mutual opposite forces, that is the forces exerted by one mass of the pair on the other. We remark further that in our model, the point masses move remaining strictly adherent to the surface, i.e. they are allowed to slide down the surface, but can neither leap nor roll.

In this paper, we explore the dynamics of the \(N\) point-masses mechanical system under the assumption that the internal forces can be of two kinds, we denote as rigid-body (RB) forces and elasto-plastic (EP) forces, though the formulation is more general and can cover also other interaction laws. As RB forces we mean forces that act to keep the distance between two interacting masses constant in time. The reason for such denomination is that a rigid body does not change shape while moving, which implies that any two points of the body will remain at the same distance at any time. However, we point out that the systems dealt with by the model UBO-Inter (where “UBO” stands for University of Bologna and “Inter” stands for Interaction) do deform while sliding, since they adapt themselves to the geometry of the sliding surface, and that this can occur even if the masses interact only through forces of the RB type. To better understand this issue, let us consider a simple system formed by 4 point masses that lie on a sliding surface and that are located at the vertices of a quadrilateral. Let us further assume that interaction takes place through RB forces only between pairs of masses that lie on contiguous vertices, but not between pairs at the ends of diagonals. As the effect of the RB forces, the quadrilateral will conserve the length of its four sides but could change shape since one of its diagonals might get shorter while the other might become longer. In all the applications presented in this paper, we will consider systems of \(N\) point masses that occupy the vertices of a grid formed of quadrilaterals and we will assume that RB forces are active along the sides of the quadrilaterals whereas EP forces are active along the diagonals. As we will see, the effect of the along-diagonal EP forces is that of allowing the quadrilaterals to deform, but of avoiding that they may degenerate to a single line. We can refer to such systems as quasi-rigid bodies and we believe that they provide a good representation of rockfalls dynamics.

Let us consider \(N > 2\) point masses that are set over a surface \(\varSigma\) given by \(z = f\left( {x,y} \right)\). The mass positions are given by the coordinates \(\vec{r}_{i}\), in a \(\left( {x,y,z} \right)\) Cartesian reference system. We assume that the point masses are virtually connected through a grid, forming triangles or quadrilateral elements, and interact with each other through forces acting along the sides of the grid. Therefore, the masses occupy the nodes of the grid, and the grid is immaterial. Nonetheless, for convenience of illustration, we will refer to such a system as an \(N\)-element grid-based body, or simply as a body. The reason for this is that the \(N\) point masses of our system can be viewed as the centers of mass (CoMs) of \(N\) elements that altogether form a partition of a continuous body. More precisely, since all points lie on the surface \(\varSigma\), then they can be seen as the projection of the element CoMs on \(\varSigma\). Consequently, we assume that studying the dynamics of an \(N\) point-mass system on \(\varSigma\) allows one to study the evolution of the corresponding continuous body that moves on that same surface. As mentioned in the previous introductory section, we simplify the inter-element interaction process by introducing RB and EP forces. A sketch with the subdivision of a sliding body in interacting elements, point masses, and grid sides is shown in Fig. 1 for a generic body shape. Though our approach is general, in this paper we consider virtual immaterial grids formed by triangles or quadrilaterals.

Fig. 1
figure 1

Planar view of a generic body. The thick solid line is the boundary of the body formed by several irregular elements. The point masses represent the CoMs of the elements, and the interaction forces act along the dashed lines joining the CoMs. The resulting immaterial grid is composed of triangular elements

Let us start the formulation of the problem from the equations of motion of the point mass i (\(1 \le i \le N\)), that can be expressed in the compact form:

$$m_{i} \ddot{\vec{r}}_{i} = \vec{F}_{i} + \mathop \sum \limits_{j} \vec{h}_{ij} + \mathop \sum \limits_{j} \vec{p}_{ij}$$
(1)

where \(m_{i}\) is the mass, \(\ddot{\vec{r}}_{i}\) is the acceleration and \(\vec{F}_{i}\) is the resultant of the external forces. The last two terms represent respectively the contributions of the RB forces that conserve the distance between pairs of point masses i and j (\(\vec{h}_{ij}\)) and of the EP forces that do not conserve it (\(\vec{p}_{ij}\)). The summations extend over the subset of all point masses interacting with the point mass i. Both internal types of forces, by assumption, respect the reciprocity law, that is \(\vec{h}_{ij} = - \vec{h}_{ji}\) and \(\vec{p}_{ij} = - \vec{p}_{ji}\).

The external forces include gravity, which is the driving force of the whole sliding process, the reaction of the surface \(\varSigma\), the surface friction, and, when appropriate, the buoyancy force. They can be given the form:

$$\vec{F}_{i} = m_{i} \vec{g} + m_{i} \vec{N}_{i} \left( {\ddot{\vec{r}}_{i} - \vec{g}} \right) \circ \hat{n}_{i}$$
(2)

where \(\hat{n}_{i}\) represents the unit vector normal to \(\varSigma\), pointing upwards, and \(\vec{N}_{i} = \hat{n}_{i} - \mu \frac{{\dot{\vec{r}}_{i} }}{{\left| {\dot{\vec{r}}_{i} } \right|}}\) is the normal vector corrected for the effect of the basal friction (\(\mu\) being the basal friction coefficient) acting against the velocity \((\dot{\vec{r}}_{i} )\). Notice that, when buoyancy is active, the expression (2) still holds, provided that the mass \(m_{i}\) is replaced by the reduced mass, as will be clear in the application of the method later on. Further, the dot represents a scalar product. The drag force due to the ambient acting on the system will be introduced later. Eq. (2) can be rewritten in the more convenient form:

$$\vec{F}_{i} = - m_{i} g\hat{k} + m_{i} \left[ {\left( {g\hat{k} + \left( {\dot{\vec{r}}_{i} \circ \dot{\vec{r}}} \right)/R_{i} } \right) \circ \vec{N}_{i} } \right]$$
(3)

where \(\hat{k}\) is the vertical unit vector pointing upward, \(g\) is gravity, and \(R_{i}\) is the local radius of curvature of \(\varSigma\). The latter expression can be deduced employing Eq. (2) by splitting the acceleration \(\ddot{\vec{r}}_{i}\) in its components. Indeed, the normal acceleration and the velocity of the point mass are linked together through the surface derivatives:

$$\ddot{\vec{r}}_{i} \circ \vec{n}_{i} = \frac{{\dot{\vec{r}}_{i} \circ \dot{\vec{r}}_{i} }}{{R_{i} }} = \dot{\vec{r}}_{i} \circ \left( {\dot{\vec{r}}_{i} \circ \nabla } \right)\hat{n}_{i} = \frac{{v_{ix}^{2} f_{ixx} + 2v_{ix} v_{iy} f_{ixy} + v_{iy}^{2} f_{iyy} }}{{\sqrt {1 + f_{ix}^{2} + f_{iy}^{2} } }}$$
(4)

where \(\nabla\) is the gradient operator. Here the velocity components in the \(x\) and \(y\) directions are denoted by \(v_{ix}\) and \(v_{iy}\), respectively. \(f_{ix}\) and \(f_{iy}\) represent the surface first derivatives in the position of the point mass i. Analogously, \(f_{ixx}\), \(f_{ixy}\) and \(f_{iyy}\) denote the second derivatives.

In order to simplify the illustration of the problem, let us consider that example shown in Fig. 2, where the \(N = 12\) masses of the system are set in the nodes of a grid formed of \(N_{Q} = 6\) quadrilaterals. We further assume that in each quadrilateral the side edges are non-deformable (which means that the pair of masses at the end of these edges do not change their initial distance), while the diagonal edges can change length according to an elasto-plastic law (which means that the corresponding masses are allowed to come closer or to move away from each other).

Fig. 2
figure 2

Sample of \(N = 12\) point masses placed at the nodes of a grid with \(N_{x}\) points in the \(x\) direction and \(N_{y}\) in the \({\text{y}}\)-direction. The RB forces \(h_{j}\) (\(j = 1:N_{E} )\) are assumed to act on the side edges (solid lines), while the EP forces \(p_{k}\)(\(k = 1:N_{D} )\) are assumed to act on the diagonals (dashed lines) of the \(N_{Q}\) quadrilaterals. Notice that during the motion the structure changes shape to adapt itself to the sliding surface but remains topologically invariant. This simple semi-rigid structure is useful to understand the forces acting on the system. Nonetheless, it is merely one of the possible configurations

In general, if a grid like this is formed by \(N_{x}\) and \(N_{y}\) point masses in the respective \(x\) and \(y\) directions, it results that the total number of point masses is \(N = N_{x} N_{y}\), the number of quadrilaterals is \(N_{Q} = \left( {N_{x} - 1} \right)\left( {N_{y} - 1} \right)\), the number of side edges is \(N_{E} = 3N - 2N_{x} - 2N_{y} + 1\), while the number of diagonal edges is \(N_{D} = 2N_{Q}\), if one considers both diagonals.

On the masses two kinds of interaction forces are active. The RB forces denoted as \(\vec{h}\) and acting along the \(N_{E}\) side edges are evaluated by imposing that the 3D Cartesian distance between the point masses connected by side edges are constants of the motion. It can be proven that these forces do not perform any work on the system (Gallotti and Tinti 2019; Tinti and Gallotti 2019). These forces change their magnitude during the motion, and the algorithm to evaluate them is strictly related to the solving process, as we will show later on. The second type of forces (EP), denoted as \(\vec{p}\), act along the \(N_{D}\) diagonal edges and their magnitude can be expressed by:

$$p_{ij} = \frac{{m_{i} m_{j} }}{{m_{i} + m_{j} }}g k\left( w \right)$$
(5)

where \(k\left( w \right)\) is a function depending on the relative lengthening (shortening) of the edge length joining the point masses i and j. More precisely, \(w = \left( {l - l_{0} } \right)/l_{0}\) where \(l\) is the instantaneous diagonal distance and \(l_{0}\) is its initial value. The function \(k\left( w \right)\) expresses the rheological behavior of the system material. In our present approach we assume it is given in the form:

$$k = \hbox{max} \left( {\frac{\alpha w}{{w - w_{1} }};k_{max} } \right)\quad if\quad w < w_{1}$$
$$k = \frac{\alpha w}{{w - w_{1} }}\quad if\quad w_{1} \le w < 0$$
$$k = - \frac{{k_{min} }}{{w_{2} }}w\quad if\quad 0 < w < w_{2}$$
$$k = - \frac{{k_{min} }}{{w_{3} - w_{2} }}\left( {w - w_{3} } \right)\quad if\quad w_{2} \le w < w_{3}$$
$$k = 0\quad if\quad w \ge w_{3}$$

which corresponds to a simplified form of elasto-plastic deformation. Here the values of \(k_{max}\) and \(k_{min}\) are chosen properly to reduce or enhance the degree of deformability. In Fig. 3 we show a plot of the function \(k\left( w \right)\) for the sample values: \(w_{1} = - 0.5, w_{2} = 0.5, w_{3} = 1, k_{max} = 10, k_{min} = 1.\). Observe that imposing a fixed-distance bound on all the edges (sides and diagonals of the quadrilaterals) generally would be impossible because it would lead to an overdetermined problem, giving the system no sufficient degrees of freedom. On the other hand, if we consider the other extreme, that is if we impose EP forces on all the edges of the grid, this would lead to body conditions too far from rigidity.

Fig. 3
figure 3

Function \(k\left( w \right)\) for the sample values: \(w_{1} = - 0.5, w_{2} = 0.5, w_{3} = 1, k_{\hbox{max} } = 10, k_{\hbox{min} } = 1.\) The function represents a simplified rheological behavior of the system material

To solve the system of Eq. (1), the RB interaction forces \(h_{ij}\) must be determined in advance. This is one of the key aspects of the model UBO-Inter. To this purpose, we use the bound on the length of the side edge linking the point masses i and j, that can be written as:

$$l_{ij}^{2} = \left( {\vec{r}_{i} - \vec{r}_{j} } \right) \circ \left( {\vec{r}_{i} - \vec{r}_{j} } \right) = const.$$
(6)

After time-deriving the expression (6) twice, we get:

$$\left( {\ddot{\vec{r}}_{i} - \ddot{\vec{r}}_{j} } \right) \circ \vec{r}_{ij} = - \dot{\vec{r}}_{ij} \circ \dot{\vec{r}}_{ij}$$
(7)

The difference of accelerations in Eq. (7) can be expressed in terms of the equations of motion (1) related to the point masses i and j. After substitution we obtain:

$$\left( {\frac{1}{{m_{i} }}\vec{F}_{i} - \frac{1}{{m_{j} }}\vec{F}_{j} } \right) \circ \vec{r}_{ij} + \frac{1}{{m_{i} }}\left( {\mathop \sum \limits_{n} \vec{h}_{in} \mathop \sum \limits_{n} \vec{p}_{in} } \right) \circ \vec{r}_{ij} - + \frac{1}{{m_{j} }}\left( {\mathop \sum \limits_{q} \vec{h}_{jq} + \mathop \sum \limits_{q} \vec{p}_{jq} } \right) \circ \vec{r}_{ij} = - \dot{\vec{r}}_{ij} \circ \dot{\vec{r}}_{ij}$$
(8)

Analogous equations hold for all pairs of point masses interacting via RB forces and all together they form a linear system of equations where the unknowns are the magnitudes of such forces, i.e. \(h_{ij}\). Since these forces are as many as \(N_{E}\) (that is one for each side edge of the grid), one can redefine the unknowns as \(h_{k}\) with k (\(1 \le k \le N_{E}\)), and, after some algebraic manipulations, provide the solution of the system as:

$$H = A^{ - 1} C$$
(9)

where \(H = \left[ {h_{1} ,h_{2} , \ldots ,h_{{N_{E} }} } \right]\) is the vector of the unknown magnitudes of such forces, \(A\) is an \(N_{E} \times N_{E}\) matrix whose elements can be derived from Eq. (8), and \(C\) is the vector containing all the terms in Eq. (8) that do not depend on \(H\). Indeed, every row of the matrix \(A\) corresponds to a couple of masses in the system that interact through forces of type RB and are thus connected by a constant-length side of the grid. For a more detailed explanation one can refer to Gallotti and Tinti (2019), and Tinti and Gallotti (2019), where the basic process leading to the linear system with unknown \(H\) and to the above inversion is explained in detail for couples of interacting masses.

Once the internal forces are known, one can solve the equations of motion. In the model UBO-Inter, these are solved by means of an explicit fourth-order Runge–Kutta scheme, implemented through a MATLAB code. These equations imply first- and second-order derivatives of the sliding surface (see Eq. (5)), that are known analytically only in theoretical cases. In real applications where digital topographic databases are available, the space derivatives can be computed only by employing discrete differential operators. Though this is not explicitly treated in this paper, an optimised algorithm to compute space derivatives is incorporated within the code UBO-Inter, with the double aim of reducing the computational time and of reducing swift derivative variations along the point-mass trajectories, since these might lead to occasional instabilities.

3 An Analytical Example

In the aim of applying the model to landslide-like processes, it is essential to understand the nature of the internal forces acting on the system. Therefore, we will address our attention to the interaction forces \(h_{n}\). To this goal, we consider a simple-geometry example formed by \(N = 3\) masses (a, b, c in Fig. 4) set at the vertices of an equilateral triangle and that interact through RB forces. The masses are located in a sphere cup of radius \(R_{s}\), whose bottom point is denoted by the letter O, as shown in Fig. 4. All masses are set at the same height \(z = R_{s} \left( {1 - cos\theta } \right)\), where \(\theta\) is the colatitude (i.e. the angle between the \(z\)-axis and the line joining the sphere center with any one of the masses). The distance between the masses coincides with the sides of the triangle and is equal to \(d = \sqrt {3 } R\), with \(R = R_{s} \sin \theta\). By considering an initial horizontal velocity \(\vec{v}_{0}\) tangent to the sphere and assuming a frictionless surface, the motion is forced to occur on the circumference of radius \(R\). In this peculiar configuration, all masses are subject to a reaction force with the same magnitude, and this equality holds also for the gravity force and the interaction forces. In the direction of \(R\), the centripetal force has to balance the component of the reaction force as well as the components of the internal forces.

Fig. 4
figure 4

Front (upper sketch) and side (bottom sketch) view of the 3-mass system with masses posed at the vertices of an equilateral triangle in a sphere cup of radius \(R_{s}\). The masses are denoted by a, b, c. The gravity force and the centripetal force \(F_{r}\) are represented by red arrows and blue arrows respectively. We denote the radius of the circle where the motion occurs with \(R\) and with \(\theta\) the colatitude. For the sake of clarity, in the bottom figure, we just show the forces acting on the two front masses a and b

Considering that the interaction forces form an angle of \({{\pi }}/6\) with the radius R, we can write the following equation:

$$m\frac{{v_{0}^{2} }}{R} = F_{r} \sin \theta - \sqrt {3 } H$$
(10)

where the centripetal acceleration term is on the left side member, \(H\) is the interaction force, \(m\) is the mass and \(F_{r}\) is the reaction force exerted by the sliding surface. While running on the horizontal circumference, the mass runs also on the great circle of radius \(R_{s}\) passing through the instantaneous position of the mass. Even on this circumference the centripetal force \(F_{r}\), directed towards the center of the sphere, has to balance the sum of all active forces, including gravity. It is thus possible to write the following equation:

$$m\frac{{v_{0}^{2} }}{{R_{s} }} = F_{r} - mg\cos \theta - \sqrt {3 } H\sin \theta$$
(11)

where \(g\) is the gravity acceleration.

The two equations together form a linear system in the unknowns \(F_{r}\) and \(H\), that has the solution:

$$F_{r} = \frac{mg}{\cos \vartheta }\quad 0 \le \vartheta \le \frac{\pi }{2}$$
(12.1)
$$H = \frac{mg}{\sqrt 3 }\left( {\tan \vartheta - \frac{{v_{0}^{2} }}{rg\sin \vartheta }} \right)$$
(12.2)

This basic formula can be easily extended to a generic number of equal masses set at an equal mutual distance and spinning on the same horizontal circular plane in a sphere cup. The expressions for the RB forces differ only by a numerical coefficient. Indeed, it is easy to show that for a system of N equal masses, equally spaced, the above formula (12.2) generalizes to:

$$H = \frac{Mg/N}{{2sin\left( {\pi /N} \right)}}\left( {\tan \vartheta - \frac{{v_{0}^{2} }}{rg\sin \vartheta }} \right)$$
(13)

where \(M = Nm\) is the total mass of the system. It’s interesting to see that when the number of points tends to infinity, the system can be seen as a continuous ring of mass M and negligible cross-section, and the coefficient before the parenthesis tends to the constant \(\sigma g\), where \(\sigma = \frac{M}{2\pi }\) can be considered the linear density per radian of the ring. This example is very useful to understand the repulsive/attractive nature of the RB interaction force. Once the angle \(\theta\) is set, the force \(H\) will be positive, and then repulsive, for low values of the initial velocities. In this case, in absence of interaction, the masses would normally tend to fall towards the center O of the sphere cup, that is towards the lowest admissible point and therefore would tend to reduce their mutual distances. On the contrary, for higher \(v_{0}\) the force \(H\) is attractive (negative) because the masses would normally escape the sphere cup and would tend to increment their mutual distances. Formally \(H \ge 0\), if \(v_{0}^{2} \le rg\sin \theta \tan \theta\). In the special case where \(v_{0}^{2} = rg\sin \theta \tan \theta\), which can be designated as the escape velocity, the internal forces are all zero, and the masses move as they were independent.

The analytical expression (12.2) can be used as a benchmark for the interaction force computed by the numerical code UBO-Inter. Pointedly, we computed a \(300 s\) simulation on a sphere of radius \(R_{s}\), with the point masses set at \(\theta \cong 64^\circ\). An initial purely horizontal velocity \(v_{0}\) was imposed. The difference between the computed interaction force and the expression (12.2), properly normalized to the analytical value, resulted to be negligible, in the order of \(10^{ - 13}\). The distance between the point masses are similarly conserved during the motion showing relative errors as small as \(10^{ - 13}\). The total numerical energy is likewise conserved.

The same experiment was repeated with systems with more point masses (up to \(N = 6\)) and the numerical results were equally good.

4 The 1783 Scilla Landslide

Through the model UBO-Inter, we have simulated the 1783 Scilla landslide taken as a benchmark case. The slide detached on 6 February 1783 from Mt. Pacì some 20 min after a violent shock that was one of the long sequence of earthquakes that shocked the Calabria region, Italy, in the period 1783–1784 (Boschi et al. 2000). The landslide triggered tsunami waves that killed more than 1500 people on the close beach of Marina Grande, where they had gathered far from houses to escape the earthquake destructive shaking. There are several detailed historical accounts of this event (see e.g. the coeval sources such as Sarconi 1784; Minasi 1785; Vivenzio 1788, and the modern reconstructions, such as Tinti et al. 2004; Graziani et al. 2006; Gerardi et al. 2008). Recent geological and geomarine investigations have identified the initial detachment niche of the landslide as well as the final deposit offshore (Bosman et al. 2006; Mazzanti and Bozzano 2011; Bozzano et al. 2011).

Owing to the relevance of the case in terms of life toll and to the availability of historical observations and of recent data, this case is a reasonable choice to validate the model UBO-Inter. For the sake of completeness, following a strategy often used to validate simulation models (see for example Wang et al. 2019) we will compare the UBO-Inter simulation not only with the experimental data but also with the results of simulations carried out through a reliable landslide model. In this work, we selected the model UBO-Block 1D that was applied to the Scilla landslide case by Zaniboni et al. to investigate the tsunami propagation in the source area near Scilla (2016) and a larger area covering the Messina Straits (2019). It is worth noting that the model UBO-Block 1D (together with the more sophisticated version 2D) was validated through laboratory data and several applications to real landslide events (e.g. Tinti et al. 2006; Zaniboni and Tinti 2014).

The zone of the landslide occurrence is represented in Fig. 5 through a topo-bathymetric grid that covers an area of \(4.0 \times 4.5\;{\text{km}}^{2}\), including the Mt. Pacì slope, its surrounding, and the final deposit. It was obtained from the SRTM database for land topography and from the GEBCO database. The grids were complemented by local nautical charts published by the IIM (Istituto Idrografico Militare, the Italian Navy hydrographic institute) for the shallow-water and offshore bathymetry.

Fig. 5
figure 5

Present-day topo-bathymetry of the Scilla area. The initial sliding area is delimited by the green line. The initial thickness is shown through the black-red palette. The final detected deposit is delimited by the red boundary. The sliding surface used in the simulation by Zaniboni et al. (2016) is shown by black diagonal lines. Coordinates are in UTM-32 format

In the model UBO-Inter, where masses are considered dimensionless points, the thickness of the slide body is accounted for by varying the point masses, provided that the constraint of the total slide mass is respected: larger masses are placed in points where the slide is thicker. We assume that the rock density is uniform, more precisely \(\rho_{r} = 2200 \;{\text{kg}}/{\text{m}}^{3}\), and that the total volume of the landslide is \(V = 6.4 \times 10^{6} {\text{m}}^{3}\). This volume is deduced by filling the discovered detachment niche and reconstructing the original morphology. Note that the final deposit was seen to involve a volume slightly smaller than the total reconstructed landslide volume, but the discrepancy can be easily interpreted as the effect of mass removal due to erosion processes (see Zaniboni et al. 2016) that started to act soon after the landslide emplacement.

In Fig. 5 we show the initial landslide thickness distribution and the final deposit area. These can be considered observation data. The area that was assumed by Zaniboni et al. (2016, 2019) to be potentially swept by the landslide is portrayed as hatched in Fig. 5 and crosses the deposit. The model UBO-Block 1D approximates the landslide as a chain of 10 deformable blocks that are allowed to spread laterally within the prescribed sliding surface and whose CoMs run along a common prescribed trajectory. Given these constraints, the code computes accelerations and velocities of the blocks along the CoMs trajectory from the stage when the slide starts moving until the slide comes to rest.

In the model UBO-Inter, the N point-masses are set at the nodes of a grid formed by quadrilaterals that can deform during motion. We have tested several kinds of initial mass distributions, all under the additional constraint that the CoM of the entire system be coincident with the CoM of the slide. To build viable sets of point masses we have used the criterion that the distance between the masses is larger than, or in the order of, the resolution of the topo-bathymetric dataset (\(50\;{\text{m}}\)) since this ensures an efficient interpolation process when computing the first and second derivatives of the sliding surface (see Eq. (4)). Taken this into account, the initial sliding area determines the maximum number of masses that can be accommodated within it. In our case, the initial landslide footprint is \(1.4\;{\text{km}}^{2}\), which allowed us to vary \(N\) between 9 and 42.

In the model UBO-Inter, when the point masses move on land (i.e. when they are above the sea level) we set the basal friction \(\mu = 0.23\), which represents the friction coefficient in dry conditions. When the slide enters the sea, the basal friction diminishes, and we use the value \(\mu_{w} = 0.04\). Underwater, the slide is slowed down not only by the basal friction but also, and mainly, by the drag exerted by the water. In the model, this is accounted through an equivalent basal friction coefficient that is added to the usual friction \(\mu_{w}\) through the following expression:

$$\mu_{eq} = \mu_{w} + 0.5\frac{{\rho_{r} }}{{\rho_{w} }}\frac{{C_{d} }}{{g l_{v} }}v^{2}$$

Here \(\rho_{r} = 2200\;{\text{kg/m}}^{3}\) is the rock density, \(\rho_{w} = 1030\; {\text{kg/m}}^{3}\) is the sea water density, \(C_{d} = 1\) is the drag coefficient for a block with size given by \(l_{v}\), that is the length of the edge connecting the point masses in the motion direction and \(v\) is the relative velocity of the water with respect to the point mass. This, assuming that the slide moves in still water, identifies with the velocity of the point mass. Further, underwater one has to account for buoyancy, which is simply implemented by introducing the factor \(\alpha = 1 - \rho_{w} /\rho_{r}\) and replacing the mass \(m_{i}\) with the reduced mass \(\alpha m_{i}\) as soon as the point mass enters the sea water. The simulation has been run until \(t = 100 s\) with a time step \(dt = 0.05 s\).

In Fig. 6, we compare the time histories of the speed of the slide CoM of different mass configurations computed by means of the code UBO-Inter against the corresponding curve resulting from the model UBO-Block 1D. The initial accelerating phase of the motion is quite similar for all the tested configurations. Discrepancies can be noticed in the velocity peaks. Notably, the \(N = 9\) and \(N = 16\) systems share the same behavior and attain the same peak as the reference model (\(32\; {\text{m/s}}\) at \(t = 20\, {\text{s}}\)). For \(N = 20\) and \(N = 25\) the velocity peak (\(37\; {\text{m/s}}\)) is reached at \(t = 22\;{\text{s}}\). For the configurations with more point masses, the velocity peaks are higher (about \(42\;{\text{m/s}}\)) and are reached later (at \(t = 23\;{\text{s}}\)). The decelerating phase of the motion shows also some dissimilarities. Slides represented by fewer mass points show a quasi-linear negative speed slope (corresponding to an almost constant deceleration). Configurations with more nodes exhibit a first strong deceleration stage followed by a lower decrease rate. All the configurations come to rest at \(t = 110\;{\text{s}}\), matching the UBO-Block stopping time.

Fig. 6
figure 6

Slide CoM velocity vs. time. UBO-Block 1D (10 blocks) speed (red) is compared against the results obtained through UBO-Inter for configurations of different number of point masses (\(N = 9\), black; \(N = 16\), green; \(N = 20\), blue;\(N = 25\), gray; \(N = 36\), magenta, \(N = 42\), cyan). Curves differ somewhat, but the time of arrest matches the UBO-Block result

Judging from the results of Fig. 6, the UBO-Inter speed curves can be divided into three groups. The \(N = 9\) and \(N = 16\) systems behave similarly and match the UBO-Inter result. Configurations with \(N = 20\) and \(N = 25\) show slightly higher maximum velocity some seconds later. Eventually, the systems with a larger density of point masses show similar higher velocity maxima. In the following, we select two configurations belonging to the most different classes (namely, \(N = 16\) and \(N = 42\)) to investigate the landslide motion in more detail. We start with the RB forces, focusing on the first \(50\;{\text{s}}\) of the motion when the slide experiences the strongest deformations, due to the impact with the sea water. Figure 7 illustrates the mean of all the RB forces acting on the point masses, normalized to the weight corresponding to the average mass value, i.e. \(\rho_{r} V/N\). Both curves show a narrow peak (7–8 s wide) and a following descent, but the peak of the configuration with more nodes is smaller (\(0.6 \;{\text{g}}\) against \(0.7\;{\text{g}}\)) and delayed (\(21 \;{\text{s}}\) against \(16 \;{\text{s}}\)). It is interesting to note that the RB forces possess a considerable magnitude, being a significant fraction of the mass weight and that they can pass through large variations during the motion, ranging  from 0.2 to 1.75 times the initial value, which means that the sliding mass is subjected to very strong stress changes.

Fig. 7
figure 7

Mean RB force for the \(N = 16\) and \(N = 42\) point-mass systems (red and blue lines, respectively) in the first \(50 \;{\text{s}}\) of the landslide motion

Bearing in mind the previous results, the slide trajectory is shown in Fig. 8, for the \(N = 16\) (left panel) and \(N = 42\) (right panel) configurations. The system structure is shown at \(t = \left[ {0, 30, 100} \right] {\text{s}}\). Observe that at the \(30 \;{\text{s}}\) picture, the landslide is fully underwater. Notice that the initial shape is partially lost during the motion in both cases, which makes more manifest the deformation experienced by the landslide. The \(N = 16\) slide spreads in the underwater part of the motion with the effect that the deposit covers an area wider than the initial area. For the \(N = 42\) system, this effect is less evident, while more evident is a slight overall clockwise rotation. For both systems the final position is located in the central-eastern part of the observed deposit.

Fig. 8
figure 8

Motion of the Scilla landslide, modelled with UBO-Inter. The landslide is represented by \(N = 16\) (left panel) and by  \(N = 42\) (right panel) point masses. The initial niche is colored in cyan, while the cyan dashed line delimits the observed final deposit. Positions of the point masses (red dots) are given at times \(t = 0, 30, 100 \;{\text{s}}\). The \(100 \;{\text{s}}\) configuration practically coincides with the end of the simulation and can be taken as the final deposit. Altitude values are given in meters

5 Discussion and Conclusions

In this study, we have presented a new approach (the model UBO-Inter) to simulate the motion of rockslides based on the dynamics of systems of interacting point masses that are subjected to external and internal forces. This model is suitable for mildly deformable bodies and couples the system dynamics with the structure geometry.

We consider two types of interaction forces between pairs of point masses: the RB forces maintain constant the distance of the masses of the pair like in a rigid body, while the EP forces contrast distance changes through an elasto-plastic behavior. The nature of the link (either RB or EP) for pairs of interacting point masses has to be specified a priori. This allows great flexibility in describing the behavior of complex mechanical systems like landslides since it influences the level of deformability of the body. Some caution is needed, however, in assigning the link type of interacting pairs. Indeed, if the number of the RB links is too high, it follows that the fixed-distance constraints exceed the degrees of freedom of the \(N\) point-mass mechanical system, which leads to the impossibility of any motion and the inconsistency of the mathematical formulation. On the other hand, if the EP links are too many then large deformations can occur leading to the coalescence of point masses and other irregularities that do not reflect the behavior of the sliding rock. A convenient and simple way to describe a semi-rigid body and to avoid these problems is to set the point masses at the nodes of a grid formed of quadrilaterals where the side edges are fixed (RB), and the diagonal elements are elasto-plastic (EP). This is the structure we have used in the application to the 1783 Scilla landslide case.

The nature of the interaction forces of the RB type has been clarified through a simple analytical example dealt with in Sect. 3. We have set three point masses at the vertices of an equilateral triangle inside a sphere cup on the same horizontal plane. By imposing the same initial horizontal velocity to each mass, we describe a motion where the RB forces acting along the edges of the triangle have exact analytical expressions. Furthermore, due to the specific balance between external and internal forces, this example allows the reader to understand the repulsive/attractive nature of the RB interaction. It has been found that numerical and analytical results show a perfect match.

In Sect. 4 we have simulated the Scilla 1783 event, a case widely studied in the literature. Among the various investigations, we have selected the simulation performed via the block-based model UBO-Block 1D in Zaniboni et al. (2016, 2019) as a reference study to test the performance of the model UBO-Inter. The reference model depicted the Scilla landslide as a chain of 10 blocks with prescribed lateral spreading and prescribed slide CoM trajectory.

The model UBO-Inter can compute at any time step the kinematical quantities (position, velocity, and acceleration) of all constitutive point masses as well as the internal forces of type RB and EP linking the interacting pairs. We remind that the EP forces are computed through the simple formula (5), while the computation of the RB forces involves the setting up and inversion of a time-dependent matrix (see Eq. (9), which is a more complex and time-consuming process. It follows that the trajectories of the point masses and, as a consequence, of the CoM of the system are the result of the computations and are not prescribed a priori by UBO-Inter.

In this study, we have used a number of \(N\) point-mass configurations with \(9 \le N \le 42\). The performance of the UBO-Inter simulations is judged considering two elements: the position of the final deposit (that derives from observations) and the comparison of the kinematics of the UBO-Inter configurations against the result of the UBO-Block 1D. As regards the final deposit, we have found that in all the cases we considered, the point masses end their trajectories within the area that was identified as the depositional area of the landslide. We remark that this latter was only broadly defined by geomarine surveys (see discussion in Zaniboni et al. 2016) and therefore it has to be taken more as an indicative rather than a strict constraint. Therefore, even if the point masses do not cover uniformly the depositional area, but only a part of it (see Fig. 8), the results are quite satisfactory.

As regards the slide kinematics, we selected the velocity curve of the CoM of the landslide as a term of comparison, since it permits a very easy judgement. The main result is that all curves have approximately the same time length (about \(110 \;{\text{s}}\)) as the reference model and show similar behavior with a quick velocity increase, a single velocity peak, and a lower velocity decrease phase. Going in more detail, we see that the configurations with fewer point masses (\(N = 9, N = 16\)) have velocity curves showing a better match with the UBO-Block 1D simulation, especially in catching the height and timing of the velocity peak. Instead, configurations with more nodes (\(N \ge 20\)) overestimate the slide velocity maximum and delay it by some seconds. Noting that these slight dissimilarities occur after the landslide enters the sea water, one can conclude that they are related to the way the model UBO-Inter computes the water drag, that is related to the system geometry and is incorporated in the basal friction. Though discrepancies are minor, nonetheless they deserve attention and need to be taken into account in the future implementations of the model. It is relevant to underline that differences in the slide velocities do not influence the trajectory of the landslides considerably.

Indeed, from Fig. 8 one can see that the \(N = 16\) and \(N = 42\) systems show a similar path.

As regards the interaction forces, they are not computed explicitly by the UBO-Block model, since the lateral mass spreading is predefined geometrically and the main utilization of that model is to provide a reasonable output for tsunami computation (mostly influenced by the velocity and the geometry of the landslide underwater). Therefore, no comparison with the UBO-Inter results is possible. Nevertheless, it is relevant to stress that this code can compute the RB and EF forces as a function of time. The example shown in Fig. 7 serves only to show that they can change remarkably during the sliding process and in a sense, they influence the motion of each point mass and are influenced by the motion itself.

To summarize, we conclude that: (1) the model UBO-Inter solving in a Lagrangian way the equation of motion of point masses interacting in pairs through rigid-body-like (RB) and elastoplastic (EP) forces is adequate to describe the dynamics of semi-rigid (i.e. mildly deformable) landslides; (2) a simple analytical example involving RB forces has proven that the code computes these forces exactly; (3) the model can simulate the Scilla 1783 landslide adequately even with a small number of point masses (see the good agreement obtained with \(N = 9\) and \(N = 16\) systems).

Further studies will refine the way the drag forces are treated in the model, and great efforts will be devoted to extending the systems from 2D to 3D, allowing one to handle cases of thick landslides with point masses distributed into a volume, rather than a surface.