Mass movement phenomena

In Colombia, mass movement phenomena occur relatively frequently (^{Cuanalo-Camposa & Gallardo-Amaya, 2016}) and mainly in rainy seasons. In the last few years, it has become important to conduct research on these phenomena due to the alleged increase in the frequency of their occurrence due to increased rainfall or to the perception of increased loss of life or materials each time they occur (INGEOMINAS-CVC, 2001). Not all cases of mass movement phenomena are caused by the action of nature; anthropogenic action also contributes to their occurrence, to the extent that “human activities alter landscapes’ evolutionary processes and shape their natural forms” (^{Castro-Garzón et al., 2020}, p. 78).

A mass movement or mass wasting phenomenon, can be understood as any movement of soil mass by the action of the force of gravity (^{SGC, 2016}) with a downward character (^{Lara & Sepúlveda, 2008}). These movements involve slope-forming materials under the influence of gravity, without the primary assistance of any fluid transport agent (^{SGC, 2017}). The flow of soil dissolved in water, which is known as an avalanche, should also be considered to be a mass movement phenomenon (^{Quesada & Feoli, 2018}).

The classifications of mass movement phenomena are based on morphological aspects of the fault surface, the area of occurrence, movement velocity, volume moved, and materials involved, among others (^{Alcántara & Ayala, 2002}). The most accepted and well-known classification is the one formulated by Varnes in 1978, which uses the type of movement and the material involved as criteria (^{Fell et al., 2008}).

It is important to determine the threat points of mass movement phenomena in order to regulate land use in territories and define prevention mechanisms (^{Fell et al., 2008}; ^{Aristizábal et al., 2010}). Depending on the level of precision required for the evaluation and zoning of said phenomena, there are several methodologies that always relate internal factors of the slope and external factors that trigger the occurrence of the events (^{Guzzetti, 2000}).

According to the Colombian Geological Survey (SGC), the geological, geomorphological, land cover, and land use aspects that condition landslide hazards are inherent to each study area and establish the “predisposition of the terrain to the occurrence of mass movements” (^{SGC, 2017}, p. 19).

Considering the urgency of contributing to the solution of environmental conflicts and disaster risk management (Perico-Granados et al., 2021), the correct mapping of phenomena associated with the risk of slope instability must contemplate information related to the “spatial and temporal probability of occurrence, type, magnitude, speed, travel distance and limit of retrogressive advance of mass movements” (^{Soeters & Van Westen, 1996}). However, in practice, these aspects are not taken into account due to lack of resources for a detailed characterization of the slope.

The SGC (formerly INGEOMINAS), which is part of the National System of Science, Technology, and Innovation (SNCTI), classifies the FRM zoning methods (^{Castro et al., 2014}) as implicit, semi-implicit, empirically explicit, semi-analytically explicit, or analytically explicit according to the way in which the internal or external factors are synthesized to obtain the landslide hazard map.

The Mora and Varhson methodology

This methodology was developed in Costa Rica by Sergio ^{Mora and Wilhelm-Guenther Varhson in 1994} to classify landslide hazard in a given area or region using morpho-dynamic terrain indicators (^{Mora & Vahrson, 1994}; ^{Alcántara & Ayala, 2002}). The methodology requires five factors; a combination of three of them (slope, lithology, and soil moisture) are defined as intrinsic factors and the rainfall and seismic intensity of the study area as extrinsic factors.

The SGC classification is an explicit empirical method because internal and external factors are represented by a numerical value assigned by an expert in order to indicate their degree of contribution to ground instability. These values are totaled and divided by ranges to determine the different threat levels (^{INGEOMINAS-CVC, 2001}).

The combination of these factors is expressed as Equation (1), where H stands for hazard or degree of landslide susceptibility:

In this Equation, Sr is the relative relief, Sl the lithology, and Sh the soil moisture. The triggering or external factors are seismicity (Ds) and precipitation intensity (Dp).

Landslide susceptibility due to seismicity and precipitation intensity can be established by multiplying the product of the intrinsic factors (Sr, Sl, Sh) by each of the external factors (Ds, Dp).

The relative relief factor (Sr) involves the topographic conditions of the terrain regarding the elevations and slope angles on the surface. The methodology establishes classification ranges for the different slope values that are finally represented in numerical values from 0 to 5.

The lithology (Sl) establishes the influence of the different types of soils and rocks on the occurrence of landslides. Based on the geological maps of the area, the different lithological formations are taken into account and, according to their characteristics, they are assigned a numerical value ranging from 0 to 5.

The soil moisture (Sh) factor aims to quantify the water retained in the soil. It indicates the amount of infiltrated water that remains in the subsoil, filling the soil particles. This parameter is calculated by using monthly rainfall averages from pluviometric stations located in the study area (Table 1).

The monthly averages are assigned the values in Table 2 and totaled for the twelve months of the year, thus obtaining a value between 0 and 24 units, later assigning the value of the soil moisture parameter (Sh), as established in Table 2.

Note that the soil moisture factor does not depend on any geotechnical variable. For the seismic activity (Ds) factor, the use of the modified Mercalli intensity scale is considered to be adequate, and, for the precipitation (Dp) factor, a value between 1 and 5 is selected, which corresponds to the maximum daily precipitation value, as shown in Table 3.

Mathematical formalizations

In geotechnics, soil consists of solids, water, and air. The weight of the latter is neglected. The sum of the volumes of water and air are defined as the ‘void volume’ (^{Berry & Reid, 1993}; ^{Braja, 2015}). Figure 1 shows this composition, as well as the source that originates them. Part of the water that precipitates is incorporated into the soil through runoff, and sometimes the air enclosed between solids passes into the atmosphere.

Source: ^{Quintana-Cabeza (2014)}

Figure 1 Graphical balance of water and air in the soil mass 

In geotechnics, the moisture (ω) of a soil is calculated by the ratio of the weight of water to the weight of soil, while the saturation (S) in a soil is expressed as the ratio of the volume of water to the volume of voids. This is shown in Equations (2) and (3).

Thus, there can be a saturated soil (S=100%) with a moisture content of 60%, just as there can be a soil with a moisture content greater than 100%, where the weight of water is much greater than the weight of solids in the soil mass; in this case, however, due to the large void volume, the soil does not reach 100% saturation.

The infiltration velocity (Vs) is the rate of change for the actual length of advance of water through the void volume of the soil mass per unit time, as shown in Equation (4). The discharge velocity (v) is the amount of water flowing per unit time through a total cross-section of soil perpendicular to the direction of the flow, and it is determined by the scalar product of the permeability coefficient (k) and the hydraulic gradient (i) obtained by laboratory tests according to Equation (5). For the presented model, the hydraulic gradient is 1.

Runoff, understood as rainwater flowing freely across the surface of a terrain, is estimated to occur if the intensity of the rainfall is greater than the permeability coefficient of the soil, that is, if the speed with which the water manages to infiltrate the soil is less than the intensity of the rainfall that reaches the soil surface (in the same unit of time).

Under the assumptions of this model, soil sliding on the slope is estimated to be caused by fatigue in the shear strength of the solids comprised by the soil mass, but the liquid limit (LL) of soil consistency, proposed at the beginning of the 20th century by Swedish pedologist Albert Mauritz Atterberg, is also taken into account as a presumed trigger for sliding (Figure 2).

Source: ^{Quintana-Cabeza (2014)}

Figure 2 Theoretical consistency of the soil mass according to the moisture content 

Sliding of soil on a slope due to shear strength fatigue could result in rotational and translational landslides, while exceeding the liquid limit could result in avalanche- or reptation-type landslides (Ramos et al., 2015).

The fatigue of the material (soil solids) is calculated, second by second, by a mathematical expression that sums two relationships: a) cohesion (C) with the product of the saturated specific weight of the soil (γ_sat), the saturated depth (H), and the slope slope angle (ß); and b) the product of the ratio of the effective specific weight of the soil (γ^') and its saturated specific weight affected by the ratio between the tangents of the angle of internal friction (Ø) and the slope (inclination) (^{Braja, 2015}).

Although the model does not aim to analyze the stability of the slope, the mathematical expression used for the calculation of the shear strength of soil solids defines the stability of infinite slopes with infiltration, as shown in Figure 3. A slope is an artificial conditioning of the soil to anthropic needs, whereas a hillside is the natural inclination that a soil mass has acquired at the end of its natural geological processes (Figure 3).

Source: Adapted form ^{Braja (2015)}

Figure 3 Analysis of an infinite slope, without infiltration 

In this model, the expected value to which the mathematical expression tends is 1 because, at that instant, and at a certain saturated depth, the forces tangential to the plane at the saturated depth that resist the column of saturated soil are equal to the tangential forces that cause it to slide, i.e., the failure is incipient (initiation). In slope stability, the employed safety factor is 1,5 (^{Ministerio de Ambiente, Vivienda y Desarrollo Territorial, 2010}).

The rainfall intensity values used in the model come from the intensity-duration-frequency (IDF) curves calculated by the Colombian Institute for Hydrology, Meteorology, and Environmental Studies (IDEAM), with data from meteorological station 16015020, which is located in Pamplona, Colombia (Lat. 7.373111 N, Long. 72.645 W, Alt. 2.340 masl). These curves are summarized in Table 4.

The sampled sites correspond to known landslides within the urban perimeter of Pamplona (Colombia): 13 Atrás Cariongo, 10 Cristo Rey (low), 2 Jurado (high), 7 Jurado (low), 8 La Curva, 14 Las Margaritas, 5 Los Olivos, 16 Salida Bucaramanga, 9 Santa Marta (high), and 6 Santa Marta (low) (^{Chaparro & Pedraza, 2012}). The above sites are shown in Figures 4 and 5.

Source: ^{Chaparro and Pedraza (2012)}

Figure 4 Geology and sampling points in the municipal capital of Pamplona, Colombia 

Source: ^{Chaparro and Pedraza (2012)}

Figure 5 Urban perimeter of the municipal capital of Pamplona (Colombia) 

The samples taken from the sites affected by mass movement phenomena do not correspond to the soil removed by the phenomenon but to that adjacent to it. This was done in order to avoid obtaining data on residual cohesion, the angle of internal friction, and the permeability coefficient. The codes come from the programmed sequence of field activities for sampling, and the names obey the subjective appreciation of the authors, together with urban denominations. For the purposes of this paper, the model has been run for the dataset of sites 16 and 8. In other analyses performed by the authors, it has been found that the mass movement phenomena of these sites may be directly related to irregular anthropic interventions rather than to natural actions.

Mathematical model approach

In order to obtain an operational model of the mass removal phenomenon which could be subsequently simulated under different scenarios, the methodology proposed by ^{Aracil and Gordillo (1997)}, which is based on Forrester's system dynamics (^{Forrester, 1997}), was applied with a systemic approach. One way to systemically approach this phenomenon is to propose a causal diagram that relates the geotechnical and hydrological variables and thus observe the relationships of influence between them. With this, feedback cycles and accumulation equations of the state variables can be identified.

The vegetation cover acts as a water control mechanism within the soil, as it intercepts part of the rainwater, part of which percolates into the soil. The water that is not filtered runs off, leading to erosion, which is affected by the soil per unit runoff and is in turn influenced by soil erodibility in terms of “increased sediment entrainment and deposition” (^{Montoya-Armenta & Rodríguez-Díaz, 2020}, p. 65). The initial causal diagram is shown in Figure 6.

Source: Adapted from ^{Quintana-Cabeza (2014)}

Figure 6 Causal diagram of the landslide event 

One way to control erosion is through civil infrastructure as a means of control, which decreases the erodible soil through stabilization works, which are carried out when erosion is appreciable. In Figure 6, it corresponds to the feedback loop in the upper right corner. However, other variables such as runoff continue to affect erosion. These control measures do not have an instantaneous effect, so there may be delays in their final effects. In turn, the amount of water that infiltrates and passes into the soil depends on the water that is on the ground but is diminished by the vegetation cover, which absorbs part of that water. The effect of human settlements is noticeable, as they reduce the vegetation cover, thus allowing for a greater amount of infiltrated water.

Once the variables have been identified by means of the causal diagram, they are classified, and then a Forrester diagram is generated (^{Aracil & Gordillo, 1997}). This diagram allows the system of equations to explain the phenomena to be simulated. Table 5 lists the main variables. By way of convention, it is established that subsystem A groups the variables related to stability, B contains those associated with humidity, C is made up of those related to rainfall, D is made up of those related to soil properties, and E is made up of variables related to permeability and infiltration. Table 5 lists the variables used in the model.

The model is based on the vertical calculation of the depth of soil saturated with water (H) through the determinination of the number of times that rainwater, immersed in the soil mass, occupies the unit void volume of the soil mass (AcumSaturatedDepth). It is expressed in ‘void volume units’ because it is the volume occupied by water and air in each cubic centimeter of soil mass. Figure 7 shows the structural diagram of the model, which results in a system of differential equations for the level variables, complemented by the algebraic equations that relate the other variables, which were previously established.

Source: Authors

Figure 7 Structural diagram of the landslide phenomenon 

According to the Forrester diagram, soil moisture varies with respect to the amount of rainwater (TotalRain) which is added to the amount of water already in the soil mass (WeightWaterUnit). The volume and weight of solids per unit volume is constant (SolidsUnitWeight), and the total solids value varies (TotalSolids) as the depth increases because the volume of soil mass to be analyzed increases.

The MOISTURE ALERT variable acquires a value of 1 if the moisture of the soil mass as the simulated saturated depth reaches the value of the liquid limit (LL), since it is assumed that, at that moment, the soil acquires a liquid consistency and starts its displacement. The model does not estimate any displacement length or velocity.

In sector A, the model calculates the imbalance that exists between the shear forces at the simulated saturated depth. The STABILITY ALERT variable acquires a value of 1 if there is no longer an imbalance of shear forces because it is assumed that, at that moment, the shear force resisting is less than the action of the normal component of the shear forces, which is due to the weight of the saturated soil column on the failure plane at the saturated depth.

Water infiltration into the soil is achieved by successive saturation of the surface unit volumes, which causes the water to move vertically downwards (to depth) by gravitational attraction due to the release of interstitial pressure. It is estimated that the depth of soil saturated with water depends only on the amount of rainwater that actually manages to enter the soil after being intercepted by the vegetation cover.

This leads to infer that a 100% saturation is essential because, without saturation, there is no vertical advance of water in the soil mass; if there is no saturation, one would be faced with a soil mass with such a high permeability coefficient that it does not experience pore pressure (interstitial pressure) due to the transit of rainwater between its empty spaces (volumes) without any problem.

Therefore, in sector E, we have the variables for calculating the infiltration velocity (VelocInfiltration), which is compared to the rainfall intensity (Intensity) in order to define the presence of runoff, since this physical phenomenon reduces the amount of water entering the soil mass. In sector C, the amount of water reaching the soil surface is calculated, understanding the soil surface as the boundary between the atmosphere and the soil mass. The values of solids and water weight, as well as the values of water volume and voids, were analytically calculated with processed data and information from laboratory tests. The model was implemented for the data set of two sites: 16 and 8.