Under natural conditions, slopes can generally be found in different dynamic states: (group 1) stable slopes with no prior or current movement, (group 2) creeping slopes, which exhibit permanent small displacements, with periodic intermittent sudden increases in the movement of the whole slope or some parts of the slope, and (group 3) the so-called conditionally stable slopes, whose stability strongly depends on their initial conditions. From the point of view of nonlinear dynamics, the conditions that leads to the instability of stable slopes (group 1) were first studied by Davis (1992), who initially proposed that the dynamics of the two-block system sliding on an inclined slope can be described by a system of three ordinary differential equations. In particular, Davis (1992) proposed such a model for certain classes of debris flows to explain the episodic surge motions. However, he was unable to estimate the onset of surging because he could not accurately determine the interaction between the steep feeder section of the slope and the lower and flatter accumulation section of the slope. In this paper, Davis (1992) assumed a velocity-dependent frictional force and suggested that there is a time delay between the movement of the feeder and the accumulation section of the slope. Our previous research, Kostić et al. (2023b), Kostić et al. (2023a), and Kostić and Stojković (2023) started from the model by Davis (1992) and introduced the explicit effect of time delay and random background noise Kostić et al. (2023b), with variable friction law Kostić et al. (2023a) and under the influence of river level oscillations, which were modeled as an external effect of colored noise Kostić and Stojković (2023). In the paper Kostić et al. (2023b) we first justified the introduction of random noise into a landslide model based on the results of recorded ambient noise within an existing landslide in Serbia. The results of further analysis showed that the background noise has a positive effect on the dynamics of the landslide, i.e., the increase in noise intensity requires a lower deformability of the slope and a higher displacement delay for bifurcation to occur. This confirms the positive stabilizing effect of increasing the noise intensity on the dynamics of the analyzed landslide model. In Kostić et al. (2023a), we first excluded the presence of significant noise in the observed landslide displacement. In a further analysis, we examined the effects of delayed failure and different friction laws on the change of dynamic regimes. It appears that under homogeneous geological conditions, sliding stabilizes with the increase of friction parameters, while rapid irregular sliding along the slope occurs with the increase of certain friction parameters. In the paper Kostić and Stojković (2023) we first show that the noise part of the deterministic model for the oscillation of the Kolubara River and the Ibar River (Serbia) has the properties of colored noise. In further analysis, it is found that landslide dynamics is sensitive to the change in noise intensity and that the increase in noise intensity leads to the onset of unstable landslide dynamics, while landslide dynamics is rather robust to the change in correlation time ϵ.

Concerning the dynamics of creeping slopes (group 2), some initial studies on creeping slopes were carried out by Chau (1995), who proposed a system of first-order ordinary differential equations, but where friction depends on both on the rate and state (Dieterich-Ruina friction law). This model is proposed for creeping slopes, and Chau (1995) investigated the specific bifurcation that can occur in the analyzed dynamical system depending on the state parameters of the sliding surface. He also showed that a small perturbation can lead to instability of the previously stable slope. Our earlier work Kostić et al. (2014) built on the work of Chau (1995) and further investigated the role of time delay in the emergence of deterministic chaos in the model of landslide dynamics. In particular, we investigated the dynamics of a single-block model on an inclined slope with Dieterich–Ruina friction law under the variation of two time delays and the initial shear stress, assuming that a periodic perturbation of the initial shear stress mimics the external triggering effect of long-distance earthquakes or some non-natural vibration source. The results obtained show conditions under which complex landslide dynamics occur, with a complete Ruelle–Takens–Newhouse route to deterministic chaos.

Chau (1999) also investigated the creeping slope model, in which friction depends on two state variables. Chau applied Reyn’s stability classification for the three-dimensional space and showed that the change of the nonlinear parameters of the slip surface, which can be triggered by rainfall or human activities, causes the occurrence of bifurcations in the studied dynamical system. One conclusion of his research is that, assuming a two-state variable friction law, the weakening of the velocity observed in the laboratory does not necessarily imply an unstable creeping slope with the same slip surface in the field (in the case of an infinite slope), while an amplification of the velocity does not necessarily imply a stable slope with the same slip surface. On the other hand, Qin et al. (2002) has proposed a nonlinear dynamical model for landslide evolution in the form of a system of three first-order nonlinear differential equations, where stress, displacement and precipitation change over time. The proposed model was applied to the case study of Xintan landslide, successfully reproducing the observed data.

The third dynamical state of the natural slope, i.e., conditionally stable slopes, has not been investigated so far. From a purely theoretical point of view, dynamic systems are in some cases sensitive to the change of initial conditions, i.e., this change can lead to the occurrence of global bifurcations Fuentes et al. (2011). We have already investigated such sensitivity from the point of view of theoretical neurology Franović et al. (2016). In the present work, we investigate the effect of changing the initial conditions on the onset of a complex dynamical behavior, which is supposed to correspond to an active landslide movement.

From an engineering point of view, landslides are indeed sensitive to changes in the original soil conditions. For example, Schilirò et al. (2019) conducted various flume tests in order to investigate the triggering process of rainfall-induced shallow landslides, focusing on the role of initial hydraulic conditions by changing the slope. They found that at a slope of 35°, the initial water content increases before the triggering occurs. This result and the lack of failure at gradient of 27° indicate the remarkable sensitivity of the tested material to even small variations in gradient. The data derived from the pore water pressure and soil moisture sensors point to two potential trigger mechanisms for variations in initial water content, namely, the advance of the wetting front at relatively high initial soil moisture content and the rise of a temporary perched water table, respectively, in relation to the failure mode. In general, the triggering of landslides was suggested to be strongly dependent on the initial water content and the initial slope angle. Also, initial porosity of the soil was already pointed out by Iverson et al. (2000) as an important initial condition for the occurrence of very different landslide rates. In their research, they raise a simple question: Can small differences in initial conditions cause some landslides to accelerate catastrophically and others to creep downslope only intermittently? Their analysis shows that landslides with a rate of 1 m/s can be triggered in moist sandy soil with a porosity of about 0.5, while the landslide rate drops to 0.002 m/s for the same soil with a porosity of about 0.4. Similar conclusions were drawn by Iverson et al. (2015), whose results from the Oso 2014 landslide indicate the primary dependence of landslide mobility on initial water content, initial porosity and initial sediment conditions.

Numerous previous studies have examined the stability of rock and soil slopes under various conditions, employing diverse models and utilizing a range of analytical, numerical and experimental methodologies (Xu et al., 2017; Zhang and Zhou, 2018; Zhou and Chen, 2019; Cheng et al., 2021; Shou et al., 2022; Li et al., 2023). In the present work, we start from the phenomenological spring-block model of Davis (1992), introduce a time delay and assume a velocity-dependent friction law proposed by Morales et al. (2017). We then simulate the model numerically and investigate how it behaves under the change of initial conditions. The paper is organized as follows. Section 2 provides a description of the investigated model. The methodology, results and related discussion can be found in Section 3, while Section 4 is dedicated to the main conclusions and directions for further research.

Let us consider a spring-block model on an inclined plane with angle θ, consisting of block-units with mass δm connected by springs with stiffness κδm and natural length δl which are acted upon by friction F(v)δm F v = a v 3 − b v 2 + c v . (1) A steady-state motion with sliding velocity v is achieved when the tangential component of gravity is balanced by friction, as the spring effects within the block cancel each other out: δ m g ⁡ sin ⁡ θ = F v δ m ⇔ F v = g ⁡ sin ⁡ θ . (2) This means that in the simplest case of a general motion of the block-unit i with the position p _i and the velocity v _i , Newton’s second law d v i d t = g ⁡ sin ⁡ θ + κ p i + 1 − p i − δ l − κ p i − p i − 1 − δ l − F v i , (3) reduces to d v i d t = F v − F v i + ∑ j = i − 1 j ≠ i i + 1 κ p j − p i . (4)

We will now proceed as in Vasović et al. (2016). If we go beyond the nearest-neighbor interaction, introduce a spring-effect delay τ and a Wiener noise d-term, and switch to a unit’s displacement x _i and velocity y _i = v _i − v relative to the steady-state values, we can model our system of block-units by the following stochastic delay differential equation: d y i d t = F v − F v + y i + ∑ j = 1 j ≠ i N κ x τ , j − x i + 2 d d w i d t . (5) Here subscript τ indicates that the quantity is evaluated at t − τ.

After applying the mean-field (MF) approximation to Eq. 5 (see Supplementary Appendix), we obtain the following system that describes our spring-block model. d m x d t = m y , (6a) d m y d t = − ( 3 a v 2 − 2 b v + c ) m y − ( 3 a v − b ) m y 2 − a m y 3 − ( 3 a v − b ) s y − 3 a m y s y + k ( m x , τ − m x ) , (6b) d s x d t = 2 u ; (6c) d s y d t = − 2 s y [ ( 3 a v 2 − 2 b v + c ) + 2 ( 3 a v − b ) m y + 3 a ( m y 2 + s y ) ] − 2 k u + 2 d , (6d) d u d t = − u [ ( 3 a v 2 − 2 b v + c ) + 2 ( 3 a v − b ) m y + 3 a ( m y 2 + s y ) ] − k s x + s y . (6e)

The system follows the time evolution of the MF variables.

• m _x –mean displacement of the block elements,

The friction in our model is given by F v = a v 3 − b v 2 + c v > 0 , (7) where the default values of the parameters are v = 0.1 , a = 3.2 , b = 7.2 , c = 4.8 , (8) and is shown in Figure 1. Since we would like to track the behavior of the system when the friction parameters are varied individually, positivity of the friction for the default sliding velocity F (v = 0.1) > 0 puts the following restrictions on the intervals of the friction parameters: a > b / v − c / v 2 = − 408 , for  fixed  default   b   and   c , (9a) b < a v + c / v = 48.32 , for  fixed  default   a   and   c , (9b) c > b v − a v 2 = 0.688 , for  fixed  default   a   and   b . (9c) If we wish the friction to keep the same shape (as a function of v), the friction parameters must also fulfill the following conditions:

• a > 0 so that F(v) generally increases with v;

Let us now find fixed points of system (6a)–(6e) describing the time evolution of our model. Equating all time derivatives to 0 immediately leads to m _y = 0 and u = 0 ((6a) and (6c)). Furthermore, dm _x = 0 implies m _x,τ = m _x , so that the delay term in (6b) disappears. This leaves us with the following system for the fixed point: 0 = 3 a v − b s y , (14a) 0 = 3 a s y 2 + 3 a v 2 − 2 b v + c s y − d , (14b) 0 = s y − k s x . (14c) For k = 0 the system is solved by s _y = 0 provided that d = 0. For k > 0 we have s _x = s _y /k and. 0 = 3 a v − b s y , (15a) 0 = 3 a s y 2 + 3 a v 2 − 2 b v + c s y − d . (15b) If 3av − b ≠ 0, as satisfied by the default values of the friction parameters and sliding velocity, the system is again solved by s _y = 0 and d = 0. However, if 3av − b = 0, a second non-trivial fixed point with s y = − 3 a v 2 − 2 b v + c ± 3 a v 2 − 2 b v + c 2 + 12 a d 6 a (16) appears, which can be simplified into s y = 3 v 2 − c / a ± 3 v 2 − c / a 2 + 12 d / a 6 . (17) Since s _x , s _y ≠ 0, this fixed point presumably corresponds to oscillations within the system that do not affect the motion of the center of mass. For the default friction parameters, the second fixed point appears for v = b/(3a) = 0.75. As we can see from (10), this value lies in the middle of the interval between v ₋ and v ₊, in which the friction shows a weak dependence on the sliding velocity. It is interesting to note that our friction law has up to three sliding solutions for the same friction value in a certain range of inclination angles. The behavior of the system in this regime could be an interesting topic for further research.

We will examine the behavior of the system around the fixed point: m y = 0 , s x = 0 , s y = 0 , u = 0 , d = 0 . (18) For the sake of simplicity, we introduce an offset parameter δ and look at the initial conditions m x 0 = δ , m y 0 = δ , s x 0 = δ , s y 0 = δ , u 0 = δ , (19) and d ≠ 0. A system with d ≪ 1 and δ ≪ 1 is considered to be near the fixed point.

System (6a)–(6e) was solved numerically in MATLAB using the forward Euler method with time step dt = 0.001 for a dense selection of system parameters and a wide range of initial conditions δ. The late-time series of solutions were then automatically classified as describing (asymptotically) sliding, oscillatory or complex oscillatory regime. The classification was done by looking at the time series in a window t ∈ [99000, 100000]. Examples of solutions in these regimes are shown in Figure 2. The solution is classified as unstable if the oscillation amplitude increases exponentially and reaches numerical infinity within a very short time, which terminates the numerical simulation.

We see from Figure 2 that m _x generally increases with time, indicating a higher sliding velocity than the assumed steady-state value v or perhaps landslide straining. It is important to note that this change δv is typically orders of magnitude smaller than the sliding velocity v. In the sliding regime, the effect disappears when we remove the noise d. However, this is not the case in the oscillating regimes. A possible explanation could also be the non-isotropic friction law. Indeed, if the friction parameter b is switched off, the effect is drastically reduced and often even completely eliminated. The interaction with other parameters of the system, such as the delay time, deserves further investigation.

The regime of complex oscillations is considered as a phenomenological model of landslide dynamics. It is particularly interesting because it can manifest itself in various different forms. Period-doubling is the simplest example of this regime (Figure 2C). The doubling is clearly visible when comparing the phase portraits in Figure 2F and Figure 2G. More complicated examples of complex oscillations such as stick-slip and pulse-like oscillations can be seen in Figure 3.

After recognizing the possible solution types, we independently varied all parameters of the system and investigated their influence on the extent of the different regimes in the parameter space of the system. We performed this analysis both near and away from the fixed point of the system in order to work out the separate influence of the initial conditions on the dynamics of the landslide.

When analyzing in detail the transition between the different dynamical regimes of the landslide model under study (Figure 4), it becomes clear that a Ruelle-Takens-Newhouse route to deterministic chaos occurs. We start with stable sliding at low values of k and τ. As the control parameters are increased, an oscillatory regime first appears through a direct Hopf bifurcation, followed by a period-doubling bifurcation and finally the emergence of a complex oscillatory regime.

The role of k and τ could be further investigated for very high values of τ. We have found that for extremely high values of τ, the change of k leads to a regular periodic stick-slip behavior. Nevertheless, it should be noted that such behavior could be considered transient, at least in some cases, since the widening of the investigated time window transforms the regular periodic stick-slip behavior into regular pulse-like oscillations.

As for the change in initial conditions, it has the following effects on the dynamics of the studied landslide model under investigation. From the point of view of spring stiffness, i.e., the deformability of the slope, and the time delay, distant initial conditions “stabilize” the landslide dynamics (Figure 5). In particular, when the initial conditions are away from the fixed point, a complex oscillatory regime is pushed towards higher values of k and τ.

Let us now present some analytical arguments for the qualitative behavior of the transition between the sliding and the oscillatory regime in the (k, τ)-diagram for small τ. Since m _y = dm _x /dt, the oscillations in m _x should also be reflected in oscillations in m _y , therefore we can consider Eq. 6b. Firstly, assuming that the amplitude of the oscillations is small, we can disregard higher powers of m _y . Secondly, the amplitude of s _y oscillations is typically orders of magnitude smaller than that of m _y , so we can also disregard s _y -terms which yields d m y d t = − 3 a v 2 − 2 b v + c m y + k m x , τ − m x . (20) If we consider the small τ regime and remember that m _y = dm _x /dt, we can expand m _x,τ = m _x (t − τ) as follows m x , τ = m x − τ m y + τ 2 2 d m y d t − τ 3 6 d 2 m y d t 2 + O τ 4 , (21) which transforms (20), up to O (τ ⁴), into a second-order linear differential equation k τ 3 6 d 2 m y d t 2 + 1 − k τ 2 2 d m y d t + 3 a v 2 − 2 b v + c + k τ m y = 0 . (22)

If the discriminant of its characteristic polynomial Δ = 1 − k τ 2 2 2 − 2 k τ 3 3 3 a v 2 − 2 b v + c + k τ (23) is negative, we can expect an oscillating solution, and if it is positive, a sliding solution. The discriminant can be further rewritten as a quadratic polynomial in kτ ²: Δ = − 5 12 k τ 2 2 − 1 + 2 τ F ′ v 3 k τ 2 + 1 , (24) which has one positive zero k τ 2 = 6 5 1 + 2 τ F ′ v 3 2 + 5 3 − 1 + 2 τ F ′ v 3 , (25) that is k = 4 6 − 6 5 τ 2 − 4 − 6 F ′ v 5 τ + F ′ v 2 4 6 + O τ = 0.76 τ 2 − 1.07 τ + 1.22 + O τ , (26) which defines the boundary between the sliding and the oscillating regime. For the data used in Figure 6A, this boundary is well fitted by k = 4.56 2 τ 2 + 8.5 2 τ + 1.6 4 , for 0.1 < τ < 0.5 , 5.4 1 τ 2 + 5.2 2 τ + 5.19 8 , for 0.5 < τ < 5.0 . (27) We see that (26) only qualitatively explains the sharp drop of the transition line in Figure 5 between these two regimes for small τ. It also explains the nearly vertical small-delay transition line in Figure 6 and Figure 9. Namely, since the leading term in (26) does not depend on the friction parameters, a constant k leads to a nearly constant τ, regardless of the changes in the friction parameters. This applies in particular to a, as it has the least influence on F′(v).

As far as the effect of the friction parameters is concerned, the influence of the variation of the higher values of the parameter a is similar to the effect of the slope deformability (Figure 6). However, at very low values of the parameter a, a new dynamical regime appears, namely, the unstable regime (Figure 7). The investigated system exhibits such dynamics regardless of the values of the time delay for the values of the parameter a in the range of 0.03–0.04. The influence of the initial conditions can further destabilize the dynamics of the landslide, i.e., the sliding regime with small time delay (Figure 7B) “vanishes” when the initial conditions are away from the fixed point. It should be noted that low values of a are far outside the reasonable parameter intervals (13a)–(13c), so it is possible that the unstable regime is caused by a velocity feedback loop where the system is accelerated by the increasing “unnatural” negative friction. Other parameters have also been considered at low values, but a has received more attention because its phase transition line is at low values.

Turning to the next friction parameter, it turns out that shifting the initial conditions away from the fixed point interacting with b destabilizes the dynamics of the landslide. Indeed, the oscillatory regime occurs at much lower values of b and τ when δ is closer to 0. Moreover, at higher values of b, the system enters a mixed sliding/oscillatory regime, the extent of which is largely independent of the initial conditions (Figure 8). Here, even a small change in b or τ can cause the “random” switching between the oscillatory and the sliding regime. If we look closer at Figure 8C, we see that in the case of the mixed regime λ _max remains positive even for very long but finite time series and only converges asymptotically to 0. This could perhaps explain the “chaotic” nature of the mixed regime. This is also in contrast to other regimes where λ _max ends up oscillating tightly around 0 within a finite time (Figure 3J–L; Figures 8D, E).

As for the joint effect of the friction parameter c and the change in initial conditions, the dynamics of the landslide seems to stabilize for initial conditions away from the fixed point. In particular, it can be seen in Figure 9 that in this case the complex oscillatory regime occurs for lower values of c than when the initial conditions are close to the fixed point.

The interaction of noise and change in initial conditions is shown in Figure 10. It appears that moving away from the fixed point stabilizes the dynamics of the landslide, as complex oscillations only occur for a very narrow range of parameters when δ is away from 0 (Figure 10B).

In this paper we investigate the effects of initial conditions on the dynamics of the phenomenological spring-block model. The initial model is described by a system of stochastic delay differential equations in which both the delay of the displacement and the random background noise are taken into account. By applying the mean-field approximation method, we have formulated a corresponding mean-field variant of the initial model. Such a deterministic mean-field model has qualitatively the same dynamics as the starting stochastic one. The analysis of the changes in the dynamics of the landslide model caused by the varying initial conditions was carried out numerically.

When the initial conditions start near the fixed point, the following dynamical regimes occur: (1) sliding regime, (2) oscillatory regime, (3) complex oscillatory regime, (4) unstable regime, (5) mixed sliding/oscillatory regime. The sliding regime represents the stable equilibrium regime of the creeping landslide. Regular periodic oscillations are considered as the transient unstable regime, while the complex oscillatory regime as well as unstable and mixed sliding regime correspond to the unstable landslide dynamics.

For the initial conditions away from the fixed point, the following can be observed: (1) the change in the initial conditions does not lead to a new dynamical regime, (2) the change in the initial conditions acts as a stabilizing or destabilizing factor, depending on which system parameter it interacts with. The stabilizing effect occurs when interacting with the noise intensity d, the friction parameter c, higher values of the friction parameter a and the spring stiffness (slope deformability) k. This stabilizing effect manifests itself as follows: (1) almost complete suppression of the occurrence of an unstable landslide (interacting with d), (2) partial limitation of the effect of the friction parameters a and c in the sense that an unstable landslide occurs at lower (parameter c) or higher values (parameter a) of the control parameters. A similar effect is captured by the interaction with the spring stiffness (slope deformability) k. The destabilizing effect of changing the initial conditions occurs when interacting with the parameter b and lower values of the parameter a. The destabilizing effect is shown by the fact that a transient oscillatory regime occurs at lower values of the parameter b and the stable sliding regime is completely suppressed at lower values of the parameter a.

As far as the interaction between the displacement delay and the change in the initial conditions is concerned, the change in the initial conditions only has an influence on the effect of the delay if it interacts with the friction parameter b, lower values of the friction parameter a or the noise intensity. In the case of the noise intensity, a change in the initial conditions leads to a complex oscillatory regime that only occurs at higher values of the displacement delay. For the friction parameter b, a change in the initial conditions leads to the appearance of the oscillatory regime for much lower values of the delay, while for the lower values of the parameter a, a change in the initial conditions suppresses the appearance of the stable sliding regime for low values of the time delay.

From an engineering point of view, the results of our research underline the importance of initial conditions for landslide dynamics. The nature of this effect (stabilizing or destabilizing) strongly depends on the friction parameters, the noise intensity, the displacement delay and the deformability of the slope.

In our research, we found that at extremely high values of displacement delay (which are hardly to be expected under real conditions), long-lasting stick-slip-like motion occurs (e.g., τ = 60, k = 6, δ = 10^–6, t up to 10⁵), which can be considered as the most representative phenomenological model of real landslide dynamics. Further research should investigate the possible occurrence of these stick-slip dynamics for a much smaller delay, which should bring the results closer to those expected in practice. However, we have observed short-term transient stick-slip behavior for a delay as small as τ = 5 (e.g., τ = 5, k = 10, δ = 10^–6, t ≲ 50 and τ = 15, k = 5, δ = 10^–6, t ≲ 1000).

The experimental verification of the proposed model and its conclusions could in principle be performed with a similar setup as in Pajalić et al. (2021) or Cheng et al. (2021). To compute the variables of the MF model and observe their evolution over time, a precise and fast mapping of the large number of points on the slope would be required. The initial steady-state velocity of the landslide could be controlled by adjusting the speed of the artificial rain as well as the steepness of the slope. On the other hand, the variation of the friction parameters could be realized by using different soil samples. This could be an interesting topic for future research.