Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru

Ronda, Gonzalo; Santi, Paul; Pope, Isaac E.; Vargas Luque, Arquímedes L.; Paria, Christ Jesus Barriga

doi:10.3390/geosciences14110315

Open AccessArticle

Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru

by

Gonzalo Ronda

^1,*,

Paul Santi

¹

,

Isaac E. Pope

¹

,

Arquímedes L. Vargas Luque

² and

Christ Jesus Barriga Paria

²

¹

Department of Geology and Geological Engineering, Colorado School of Mines, Golden, CO 80401, USA

²

Escuela Profesional de Ingeniería de Minas, Universidad Nacional de Moquegua, Moquegua 18001, Peru

^*

Author to whom correspondence should be addressed.

Geosciences 2024, 14(11), 315; https://doi.org/10.3390/geosciences14110315

Submission received: 25 September 2024 / Revised: 6 November 2024 / Accepted: 9 November 2024 / Published: 18 November 2024

(This article belongs to the Special Issue Landslide Monitoring and Mapping II)

Download

Browse Figures

Versions Notes

Abstract

:

Since the times of the Incas, farmers in the remote Andes of Peru have constructed terraces to grow crops in a landscape characterized by steep slopes, semiarid climate, and landslide geohazards. Recent investigations have concluded that terracing and irrigation techniques could enhance landslide risk due to the increase in water percolation and interception of surface flow in unstable slopes, leading to failure. In this study, we generated an inventory of 170 landslides and terraced areas to assess the spatial coherence, causative relations, and geomechanical processes linking landslide presence and Inca terraces in a 250 km² area located in the Ticsani valley, southern Peru. To assess spatial coherence, a tool was developed based on the confusion matrix approach. Performance parameters were quantified for areas close to the main rivers and communities yielding precision and recall values between 64% and 81%. On a larger scale, poor performance was obtained pointing to the existence of additional processes linked to landslide presence. To investigate the role of other natural variables in landslide prediction, a logistic regression analysis was performed. The results showed that terrace presence is a statistically relevant factor that bolsters landslide presence predictions, apart from first-order natural variables like distance to rivers, curvature, and geology. To explore potential geomechanical processes linking terraces and slope failures, FEM numerical modeling was conducted. Results suggested that both decreased permeability and increased surface irrigation, at 70% of the average annual rainfall, are capable of inducing slope failure. Overall, irrigated terraces appear to further promote slope instability due to infiltration of irrigation water in an area characterized by fluvial erosion, high relief, and poor geologic materials, exposing local communities to increased landslide risk.

Keywords:

landslides; Andes; irrigation; terraces; Peru; confusion matrix; logistic regression; GIS; geohazards

1. Introduction

Terracing techniques have been used globally to increase the area and quality of croplands in mountain communities. Since the times of the Incas, farmers in the remote Peruvian Andes have constructed and irrigated terraces to grow crops in a landscape characterized by steep slopes, semiarid climate, and landslide geohazards. Recent investigations have concluded that irrigation techniques could enhance landslide risk due to the increase in water saturation in unstable slopes, leading to failure [1,2]. Several studies concluded that intensive flooding irrigation was the most likely cause of landslides in agricultural regions with case studies from Peru [3,4,5,6,7,8,9] and China [10,11,12,13,14,15], among others [16]. These landslides commonly threaten local agricultural activities that rely on the terraces for continued economic development [6,17,18]. In contrast, Inca hydrological and terracing techniques have been described as generating sufficiently drained engineered slopes that would prevent slope failures if appropriate irrigation techniques were used [19], casting doubts on the hypothesis that terraces themselves may be the cause of landslides.

In the valley of the Ticsani, Moquegua district, Southern Peru (Figure 1), communities including Carumas and San Cristóbal rely on terrace agriculture as their main economic income in an area that has been severely affected by landslides. Every year, valuable communal croplands are lost to these events without sufficient knowledge of potential causes and triggers, hindering the application of effective mitigation measures. In this contribution, we address the relation between landslides and agricultural terraces by assessing the spatial coherence between both. The aim of this research was to assess potential causative relations between landslide occurrence and irrigated terraces in a 250 km² area located in the valley of the Ticsani (Figure 1). The scientific question leading the investigation is whether irrigated terraces are promoting landslides in the area. We hypothesize that a combination of excess water due to irrigation, poor geologic materials, and high relief may be behind landslide occurrence there.

Most of the published research has focused on the temporal relation between irrigation and landslide occurrence [2,5,20] or the geomechanical processes leading to failure [7,8,9,10,11], while the application of spatial statistics to assess the spatial correlation between both processes has remained largely unexplored, with some examples from Indonesia [21,22]. To gain further understanding of the potential links between these features, we applied a novel statistical method based on the confusion matrix to assess the spatial correlation between landslides and irrigated terraces. In addition, a spatially derived multivariate regression approach was used to investigate the contribution of other processes to the occurrence of landslides in the area. Finally, FEM modeling was conducted to explore the role of frequent landslide-generating geomechanical processes including excess pore water pressure and poor quality materials with decreased residual strength, both present in the study area.

2. Background

The area of interest (AoI) is part of the Southwestern Sierra, with a temperate semiarid climate characterized by an average annual rainfall of 260 mm/yr, and with rain deficiency from May to December [23]. Hydrologically, the area is part of the Río Tambo watershed, which drains to the Pacific Ocean and is composed of the local rivers Carumas and Putina (Figure 1). Route 102 connects the communities of Carumas, Sacuaya, and San Cristóbal with Moquegua, the main city of the region, located further southwest; approximately 12,000 people lived in the area in 2015 [24] (Figure 1). The main economic activity of the AoI consists of cooperative agriculture, with families living on an average USD 1000 per year obtained from annual harvests [25]. The main crops include potatoes, quinoa, corn, and alfalfa [25]. To sustain crops under a hydrologically deficient climate during three-quarters of the year in a high relief area, local farmers rely on both terracing techniques and artificial irrigation, with 94% of the water being used for irrigation [24].

2.1. Geomorphology

The geomorphology of the area is characterized mainly by fluvial, volcanic, and mass wasting landforms in mountainous terrain with average slopes of 25° and maximum slopes of 75°. Some periglacial landforms have been recognized in the higher altitude portions of the AoI. Within fluvial landforms, rivers and streams can be observed including the Carumas and Putina rivers (Figure 1). Fluvial terraces and scarps flank these rivers, and are frequently associated with landslides affecting croplands, as shown in Figure 2. Mass-wasting landforms include different types of landslides, further described in the following sections.

2.2. Geology, Volcanism, and Seismicity

The geology of the area of study consists of a sequence of Mesozoic and Cenozoic sedimentary rocks, along with modern deposits, as shown in Figure 3. The oldest rocks of the area are exposed as small outcrops of the Yarabamba monzodiorite, part of the Paleozoic coastal batholith found in the coastal cordillera. More ubiquitous are the Cachios and Labra formations, part of the Jurassic Yura group, a marine and coastal sequence including sandstones, bituminous shales and limestones that extensively outcrop in southwest Peru (Figure 3). These units are unconformably overlain by the volcano-sedimentary Matalaque Formation (andesites, dacites, and volcanic breccias) and the Toquepala Group (andesites, tuffs, and volcanic breccias) of the Cretaceous age (Figure 3). In the AoI, the spatially restricted Cretaceous Arcurquina Formation (conglomerates, sandstones, mudstones, and limestones) can be found overlying the Toquepala group through a slight angular unconformity (Figure 3).

Sandstones and conglomerates of the Paleocene age follow, part of the Puno Group, overlain by Eocene andesites and volcanoclastic rocks of the Pichu Formation (Figure 3). Further volcanoclastic rocks overlie, including sandstones and conglomerates interbedded with tuffs of the Moquegua Formation, Oligocene–Miocene in age. These rocks are overlain by friable rhyolitic and rhyodacitic tuffs of the Miocene Hualyllas Formation, and dacites, andesites, and volcanoclastic deposits of the Barroso Group (Figure 3). Quaternary rocks and deposits are mainly volcanic lavas, pyroclastic deposits including ignimbrites and fall deposits, and fluvial deposits with increased volcanoclastic content. To a lesser extent, glacial deposits can be found in the eastern portion of the AoI (Figure 3).

The main structures in the area include faults and folds with strikes trending NW–SE (Figure 3). The Carumas river and the domes of the Ticsani volcano are distributed along this main direction, suggesting that these structures concentrate erosion, volcanism, and seismicity in the area.

The region is part of the Central Volcanic Zone generated by the subduction of the Nazca plate below the South American plate. Two active volcanoes are located less than 20 km from the AoI: Ticsani and Huaynaputina (Figure 1). Ticsani, closer to the AoI, produced one of the largest debris avalanches in Peru when part of its edifice collapsed and flowed down the Putina River, generating a deposit between 250 m and 350 m thick (Figure 1 and Figure 3) [26]. The town of San Cristóbal and extensive irrigated croplands are founded over those deposits (Figure 1 and Figure 3).

Historical records document the abandonment of old settlements in the AoI, likely related to earthquake damage and suggesting that seismicity is a recurrent geohazard there [27]. In October 2005, an earthquake identified as part of the Calacoa-Moquegua seismic crisis occurred in the AoI. The seismic crisis caused structural damage to houses and facades in most of the constructed areas in Calacoa, San Cristóbal, Cuchumbaya and Carumas (Figure 1). During the seismic crisis, rockfalls and landslide reactivations were registered [27], suggesting that seismicity can be an important landslide trigger in the area.

2.3. Landslide Geohazards

The term landslide has been defined as a movement of a mass of rock, debris, or earth down a slope [28]. In Peru, landslides, popularly referred to as huaicos, occur frequently, and as reported in written media are usually associated with large life and economic losses [29]. The landslide with the highest recorded losses in Peru was the landslide-driven Yungay–Ancash debris flow, triggered by the 1970 Ancash earthquake, which buried two towns, killing more than 18,000 people [30]. Identified landslide triggers in Peru include seismicity [27], increased rainfall [31,32], and more recently irrigation [1,4,5,8,9].

Figure 3. Geology of the AoI, modified from [33,34]. Jm-ca = Cachios Fm.; Jm-p = Puente Fm.; Js-g = Gramadal Fm.; Js-l = Labra Fm.; Ki-mat = Matalaque Fm. Ki-hu = Hualhuani Fm. Kis-a = Arcurquina Fm. KsP-bc/y-mdcz = Batolito de la Costa Yarabamba monzodiorite. KP-tn,di = Tonalite, diorite. P-Pu = Puno Gp. P-Pi = Puno Fm. Nm-huay = Huaylillas Fm. NQ-b-and = Barroso Gp. Andesite. Q-gl = glacial deposits. Q-pl = coluvial deposits. Qh-al = alluvial deposits. Qh-vl-ce = tuff deposits. Qp-b-agand = Barroso Gp. Andesitic agglomerate. Qp-b-and = Barroso Gp. andesite. Qp-b-andp = Barroso Gp. Porphyritic andesite. Qp-b-da = Barroso Gp. Dacite. Qp-vl-pi = pyroclastic deposits.

The effect of irrigation on landslide occurrence has been studied in the region by several research groups. Ref. [5] studied historical images and determined movement for a set of landslides in southwestern Peru. These authors concluded that the development of large irrigation programs enhanced landslide movement in the region.

Ref. [4] studied the Siguas landslide, a 34 million cubic meter rotational landslide affecting unconsolidated fluvial sediments caused by water saturation linked to extensive irrigation in the Siguas irrigation project in Arequipa (Figure 1). Saturation was considered by the authors as an underlying cause, while seismicity could have triggered the landslide.

Ref. [1] provided a review of irrigation-induced landslides studies around the world, including cases from southern Peru. Their results showed that 80% of the landslides documented in the study correspond to dry areas with less than 600 mm/yr of rainfall, and that 95% of them happened in areas with summer seasonal rains, such as the AoI. The authors suggest that higher infiltration at warmer temperatures, causing higher effective permeability and higher saturation, coupled with irrigation, leads to increased landslides during summer rainfalls. Most of the studies reviewed included weak materials like loess, sedimentary, alluvial, and volcanic—all materials present in the AoI (Figure 3).

Ref. [8] used geophysical methods to investigate subsurface groundwater to improve landslide modeling around the Majes irrigation program in southern Peru. These authors identified an increase in groundwater levels due to irrigation that may have led to the instability of slopes producing landslides.

Ref. [9] studied the 2020 Achoma landslide, located in the Colca valley in a similar setting to the area of this study. The slide affected 40 hectares of valuable croplands and was estimated to be 14.5 million cubic meters in volume. These authors conducted limit-equilibrium and numerical modeling of the slide, with results pointing to increased groundwater levels potentially caused by deficient and leaking irrigation infrastructure.

Regarding global climatic patterns, ref. [31] have suggested that there is a correlation between landslide occurrence and increase in precipitation during El Niño years in central Peru.

2.4. Terraces and Irrigation Techniques

Terracing techniques have been used from ancient to modern times in Asia, Africa, and Europe [19]. In the Americas, terraces have been found from the southern US to Argentina [35]. The main purpose behind terracing is to increase cropland area in regions of high slopes, reduce surface flow and soil erosion, and increase water retention [19]. In the AoI, non-technified irrigation techniques are used, such as flood irrigation.

Terraces are defined as a horizontal surface sustained by two walls aligned, one in the direction of maximum slope, and one wall cross-cutting it [19]. A typical profile of a terrace found in southern Peru is shown in Figure 4. Typical terracing techniques used in the region according to [19] include:

Excavation of trenches to bedrock or to a depth of 0.25 the height of the wall.
Piling up rocks to build the walls, with decreasing size upwards. The base rock should be flat to increase friction with foundation.
Rocks not used to build the wall due to their smaller size or roundness should be placed as fill on the internal side of the wall in an upwards decreasing-size fashion. Any bedrock outcropping inside of the area to be terraced should be removed or excavated.
Filling the sub-rhizosphere with 15 cm layers of fertile soil and compacting each before adding a new layer, then depositing the uppermost layer with fertile soil, without compacting it.
Leveling the surface of the terrace in such a way that there is 30 cm of thickness of arable fertile soil.

3. Methods

3.1. Landslide Terraces Spatial Coherence Assessment

The proposed study relied on a landslide inventory of the AoI (Figure 1). To generate this inventory, Google Earth satellite imagery was used to manually map landslides in an area of 250 km². Landslide criteria followed [36,37,38]; at least 3 evidence indicators were used to define landslides, including, in order of importance:

Loss of vegetation cover;
Oversteepened slopes;
Circular shapes possibly linked to scarp or flanks;
Cracks and irregular terrain constrained within circular shapes;
Changes in average slope related to a zone of depletion and a zone of accumulation;
River channel offset due to accumulated material.

Main landslide parameters like elevation, surface area, and classification (based on [36]) were recorded in the inventory. The landslide inventory also included landslides reported in Peru’s national geohazards repository, INGEMMET [34]. Quality assurance was conducted by performing a field visit to the AoI, to observe scarps, loss of vegetation, and changes in slopes as indicators of landslides.

In addition to the landslide inventory, irrigated features in the AoI were mapped. The main identification criteria were based on [19] and included the recognition of walls and regular patterns mimicking the topography, presence of crops, and shallower slope angles. To assess slope effects, terraces and flat irrigated areas were discriminated based on their slope, with terraces being classified as areas outlined by slopes exceeding 25°.

Once the landslide and terrace inventories were completed, a study was performed aimed at quantifying the spatial correlation between both datasets. To do so, a human-inspected semi-quantitative tool was tested based on the confusion matrix method to assess model performance in raster classification [39]. The original method is based on dividing ground-truth and modeled pixels from two raster datasets into a series of classes, shown in Figure 5, which are then used to quantify model performance according to a series of statistical parameters.

The traditional confusion matrix approach is based on comparing overlapping pixel values to assess classification performance. In addition, when trying to relate terraces and landslides, it is necessary to assess downstream proximity and not only overlap between both classes. A novel approach was used that relies on visual classification of polygon features. The following criteria were used for classification, as shown in Figure 5:

True positives (TP) included landslides that had terraces or irrigated features upstream, to a distance of 500 m, and included landslides that were partially overlain by terraces or irrigated areas whenever such features persisted upstream of the landslide.
True negatives (TN) included areas that did not have landslides or irrigated features. Lacking both features of interest, these areas were considered irrelevant for the analysis.
False positives (FP) included terraces and irrigated areas that did not have any landslides associated. These were considered a type I error or a false alarm.
False negatives (FN) included landslides that did not have any associated terraces or irrigated areas upstream and landslides that had terraces or irrigated areas built over the landslide deposits, whenever no other upstream irrigation features were present. These were considered a type II error or a missed alarm.

Every landslide polygon was classified into TP or FN. To quantify the FP, i.e., areas where terraces were not linked to any landslides, terrace polygons were used. Since polygons are highly variable in size due to the changing extent of landslides and terraces, an area-normalization calculation was performed according to Equations (1)–(3). This procedure avoided counting a larger polygon the same as a smaller polygon in the statistics calculation. Square-root of area calculations were used to avoid skewed-distribution effects on the statistics calculation.

N o r m a l i z e d T P = \frac{\sqrt{A r e a T P L a n d s l i d e s}}{\sqrt{T o t a l A r e a L a n d s l i d e s}}

(1)

N o r m a l i z e d F N = \frac{\sqrt{A r e a F N L a n d s l i d e s}}{\sqrt{T o t a l A r e a L a n d s l i d e s}}

(2)

N o r m a l i z e d F P = \frac{\sqrt{A r e a F P T e r r a c e s}}{\sqrt{T o t a l A r e a T e r r a c e s}}

(3)

With all the polygons classified and normalized by the square-root area, statistics were calculated to quantify the spatial association of landslides with terraces. Different parameters have been used to quantify performance from confusion matrices: a review can be found in [40]. Recall, precision, and F-1 factor are commonly used metrics that are presented together [40] and were calculated according to Equations (4)–(6). Precision is the number of correctly predicted landslides divided by the total number of identified landslides, whether correct or not. Recall is defined as the number of correctly predicted landslides divided by the total number of true landslides. F-1 is a combination of both [40].

P r e c i s i o n (P R) = \frac{N o r m a l i z e d T P}{N o r m a l i z e d T P + N o r m a l i z e d F P}

(4)

R e c a l l / S e n s i t i v i t y (R E) = \frac{N o r m a l i z e d T P}{N o r m a l i z e d T P + N o r m a l i z e d F N}

(5)

F_{1} = \frac{2 \times N o r m a l i z e d T P}{2 * N o r m a l i z e d T P + N o r m a l i z e d F P + N o r m a l i z e d F N}

(6)

Preliminary observations of the AoI showed that large areas within the AoI were not irrigated but still appeared to have landslides, and thus a moderate correlation was expected. Most landslides are concentrated close to the main rivers, communities, and irrigated areas. If subareas with higher landslide density were to be analyzed with this method, then higher spatial coherence would be expected.

3.2. Logistic Regression Analysis of the Landslide Inventory

It is possible that spatial coherence analysis will not completely explain landslide presence, and there may be additional factors apart from the presence of terraces that explain the presence of landslides. Logistic regression is a method frequently used to estimate the probability of occurrence of an event or dependent variable based on a series of predictors or independent variables set by the user. It constitutes a binary regression method applied in spatially referenced datasets that include a binary dependent variable such as the presence (Y′ = 1) or absence (Y′ = 0) of landslides. The main output of the analysis would be similar to Equation (7), where the dependent variable (Y′) is explained by a constant y-intercept (β₀) and a set of independent variables (V_i) multiplied by a scale coefficient (β_i) representative of the relative importance of the variable in the probability of occurrence of Y′.

Y^{'} = β_{0} + β_{1} V_{1} + β_{2} V_{2} + \cdot \cdot \cdot β_{n} V_{n}

(7)

A second output of the regression is the cumulative distribution function shown in Equation (8):

P (1) = \frac{e^{Y^{'}}}{1 + e^{Y^{'}}}

(8)

P(1) depicts the cumulative probability of occurrence of Y′, as a function of the independent variables included in the regression, and is thus considered a measure of model performance.

Logistic regression has been extensively used to produce landslide susceptibility maps and is often based on different but overlapping predictive variables that commonly include elevation, slope, and curvature [41,42,43,44,45,46]; summary presented in Table 1). Depending on the geological variability, geology has also been used as a predictor for landslides [41,43,44,46]. Other authors have argued that topographic wetness index (TWI) is representative of wetness concentration and thus landslide occurrence [42,43,46].

The natural variables included in our analysis are listed in Table 1. They were processed as raster datasets in ArcPRO, and were chosen based on the literature review and on the availability of spatially distributed datasets in remote areas such as southern Peru. A scheme representing the collection of parameters values is shown in Figure 6. Variables like rainfall distribution, crop type, and root depth were considered but not included in the analysis due to the absence of spatial distribution data within the AoI.

In this original approach, logistic regression was applied to the entire AoI to assess the change in cumulative probability, P(1), due to the inclusion of terraces as an independent categorical variable in comparison to simpler models including natural variables. P(0) for the dependent variable “landslide presence” was represented by including 60 control areas without landslides in them, as depicted in Figure 6. Control areas were manually generated as circular areas, with variable sizes, and were spread across the entire AoI (Figure 6).

The Zonal Statistics tool (Figure 6) was used in ArcPRO to obtain mean pixel values of geology, elevation, slope, curvature, and TWI for all the landslide and control area polygons in the inventory. Geology was converted into a 4 categories raster where 1 correlates to intrusive rocks (most stable geologic conditions); 2 to lavas and hard sedimentary rocks; 3 to soft sedimentary rocks; and 4 correlates to unconsolidated fluvial and pyroclastic deposits (least stable geologic conditions). Once categorized, the geology variable was treated as a continuous variable for mean calculation. Distance to rivers was obtained by calculating the minimum distance from the boundary of each landslide and control area polygons to the closest river or stream using the Generate Near Table tool in ArcPRO (Figure 6). Distance to rivers was intended to be used as a proxy for river erosion, a known process that promotes landslides. However, river erosion is also likely affected by other variables such as geology and groundwater conditions.

Once the mean datasets had been compiled, Minitab Statistical Software was used to conduct the logistic regression. Main outputs of the regression included a regression equation and performance statistics such as p-values for each parameter. p-values determined if a given parameter was significant for the regression, and parameters yielding p-values > 0.05 were considered not significant for subsequent regressions. For each parameter, a coefficient β (Equation (7)) was also calculated, considered an indicator of the importance of the independent variable in the regression.

Overall model performance of the regression was assessed by calculating the area under the curve (AUC) of the receiver operating curve (ROC), with highest values indicating better model performance. In addition, Akaike’s information criterion (AIC) was used to assess regression performance, with smaller values indicating better performance.

The use of this methodology is not intended to generate a prediction model to be applied in other areas of the region but rather to assess the sensitivity of the prediction performance to the inclusion of human-made variables like terraces to a set of natural predictors. Hence, the method was applied to the entire AoI, whereas in typical applications a subset is used to train the model and the rest to test its performance.

3.3. Finite Element Modeling of the San Cristobal Landslide

The methods proposed so far constitute spatial statistical approaches to assess the correlation between landslides and variables such as the presence of terraces, geology, and slope, among others. To explore the mechanisms behind landslide occurrence in the AoI, finite element numerical modeling was conducted. The San Cristobal landslide was selected for modeling, which occurred on a slope bounded by the Putina River, located below irrigated croplands and close to one of the main populations of the area (Figure 7), putting at risk valuable croplands, infrastructure, and lives.

The geology surrounding the San Cristobal landslide consists of Cretaceous volcano-sedimentary rocks and debris avalanche deposits topped by pyroclastic deposits (Figure 7). The debris avalanche deposits that crop out on the scarps of the Putina River (Figure 7, Q-Pl) were the main constitutive material in the model. Other materials included pyroclastic deposits (Figure 7, Q-vl-pi) and volcanoclastic deposits of the Matalaque Fm. (Figure 7, Ki-mat).

The selected materials were modeled using a Mohr–Coulomb constitutive model, appropriate for predicting soil failure. The main inputs included Young’s modulus, Poisson ratio, unit weight, friction angle, and cohesion. A continuum modeling approach was used with a finite-element grid, an approach frequently used to model soils without discrete discontinuities. Two-dimensional finite-element modeling (FEM) of the San Cristobal landslide was conducted in RS2 [47]. FEM has been successfully applied in the region to model landslides and slopes [7,8] and was selected over limit-equilibrium approaches that do not account for internal deformation of materials. The geometry of the model was constrained from Google Earth data, as shown in Figure 7. A pre-landslide topography was assumed to constrain model boundaries. Strength reduction factor (SRF) analyses were used to estimate the factor of safety (FoS) for each model.

Material properties were obtained from the literature based on field descriptions of the outcrop units (Figure 7) [27], rather than from subsurface sampling, which was impractical and dangerous in parts of the AoI. In the absence of material testing, sensitivity analyses were performed to assess the influence of material properties including cohesion, friction angle, Young’s modulus, and Poisson ratio. Few studies exist on the geotechnical characterization of debris avalanche (DA) deposits associated with volcanic edifice collapse (Figure 7, Q-pl), with most of them related to the modeling of the actual flow and not the properties of the final deposits. In a similar setting to that of the Ticsani volcano, with similar climatic and weathering conditions, the Cotopaxi volcano in the Andes of Ecuador produced a debris avalanche [48] with similar characteristics to the one modeled. Although material properties on DA deposits can be highly variable, the parameters obtained for the Cotopaxi DA by [48] were used as best-estimate inputs, based on the similarity of the events, to model the Ticsani DA deposit (Table 2). To model the pyroclastic deposits that outcrop towards the top of slope (Figure 7, Q-vl-pi), material properties measured by [48] in similar deposits from Ecuador were used.

Regarding the Matalaque Fm., ref. [49] described it as a volcano-sedimentary sequence mainly composed of consolidated coarse sandstones, limestones, and conglomerates with intercalated lava flows. This unit is not clearly discriminated in satellite images and its extent may have been overestimated in the 1:50,000 geologic map. For modeling purposes, it was modeled as an intact cemented sandstone, considered an average material for the sequence. The parameters used are summarized in Table 2.

Groundwater levels in the AoI are unknown. Since the slope analyzed consists predominantly of an assumed homogeneous debris avalanche deposit, only one free groundwater surface was used for modeling (Figure 7). Different groundwater levels were tested to estimate a threshold level governing slope stability, following the approach by [7]. To test for sensitivity, groundwater depth was changed between 2703 m and 3055 m, modifying the start value of 2890 m by −50%, −25%, −12.5%, +25%, +45%, covering all possible levels in the slope analyzed.

The main limitations of the slope stability modeling procedure used include restricting the analysis to a 2D transect of the investigated landslide, the usage of reference values for input materials, and the assumption that the modeled materials are homogeneous, with uncertain extents. Groundwater is also an important source of uncertainty due to the absence of data constraining its depth and distribution into one or more aquifer levels. To deal with the uncertainties and limitations mentioned, sensitivity analyses were conducted, changing groundwater levels and material parameters including friction angle, cohesion, Young’s modulus, and Poisson ratio. A complete description of the sensitivity analysis conducted can be found in [50].

3.4. Seepage Analysis

To test the influence that increased irrigation-driven infiltration would have on the stability of the San Cristobal landslide slope (below the irrigated terraces in Figure 7), steady-state seepage analyses were conducted using a FEM approach in RS2. Models were run to estimate groundwater levels after equilibrium is established, as well for slope stability to estimate FoS under changing conditions of infiltration and permeability. The main limitation of this approach was the utilization of infinite time to reach seepage equilibrium, assuming constant infiltration at the rates set.

Seepage models were set by defining vertical surface infiltration and a base groundwater level corresponding to the Putina River under flood conditions (Figure 7). Infiltration was set initially at 1 × 10⁻⁵ m/s, which would be equivalent to 31.5 mm/yr, and increased up to 6 × 10⁻⁵ m/s, equivalent to 189 mm/yr, and decreased down to 1 × 10⁻⁶ m/s, equivalent to 3.15 mm/yr. In those ranges, infiltration was changed by factors of 25, 50, 100, 200, 300, 500, −25, −50, −62.5, −75, and −90. Not counting other infiltration inputs like irrigation, the region receives an average input of 260 mm/yr from annual rainfall alone, making the infiltration rates tested conservative.

Permeability values for the different materials in the models were estimated from the literature and are summarized in Table 2. Permeability varied by a factor of 0.01 and 100 for all the materials modeled to cover typical variations for common geomechanical materials, which are assumed to vary by as much as two orders of magnitude from a presumed median value. A complete description of the seepage sensitivity analysis conducted can be found in [50].

4. Results

4.1. Spatial Coherence Analysis Results

A semiquantitative confusion-matrix based analysis was conducted for the entire AoI. In the area of study, irrigated terraces comprise 18% of the area, whereas landslides account for 3% of the area. To test the influence of slope in irrigation and landslide presence, terraces and flat irrigated areas were separated for the analysis.

Statistics obtained for the analysis included precision, recall, and F1 factor, summarized in Table 3. Precision ranged between 0.4 and 0.5, recall ranged between 0.4 and 0.6, and F1 ranged between 0.4 and 0.6 (Table 3). Maximum model performance to predict the ground truth class (landslides) from predicting classes (terraces and/or flat irrigated areas) was achieved when integrating terraces and flat irrigated areas (Table 3). Minimum performance statistics were obtained when using only flat irrigated areas as the predicting class (Table 3), suggesting that terraces resulted in slightly better predictors for landslides than flat irrigated areas. The full dataset including features classification and calculations can be found in [50].

The performance statistics obtained did not surpass 0.63 (Table 3), pointing to poor model performance and spatial coherence when predicting landslide presence from the presence of terraces and flat irrigated areas. These results agree with prior expectations that considered the existence of landslides that were not coupled to irrigation features in the entire AoI.

A higher apparent coherence between landslides and irrigated features was observed surrounding the main rivers and communities in the AoI. To assess the spatial coherence there, a second confusion matrix-based analysis was conducted limited to a buffer zone of 500 m surrounding the Carumas and Putina rivers, as shown in Figure 8. In this case, terraces and flat irrigated areas were integrated into a single class given the low variability in performance statistics obtained when discriminating them (Table 3). Model performance statistics yielded a precision of 0.6, a recall of 0.8, and an F1 of 0.7 (Table 3) that point to increased model performance when predicting landslide presence coupled to irrigation features close to main rivers and communities. False negatives, representing landslides not linked to irrigation features, showed a decrease from 0.4 to 0.6 when analyzing the whole AoI and to 0.2 when restricting the analysis to main rivers. These results agree with expected increased spatial coherence between landslides and irrigation features in areas characterized by increased irrigation, slopes, and fluvial erosion.

The results obtained for both analyses point to a decreased spatial coherence when analyzing the entire AoI, and to increased spatial coherence when restricting the analysis to the main rivers within the AoI. The main communities in the AoI, and thus irrigated fields, are concentrated along the Carumas and Putina rivers, where most of the landslides have been observed (Figure 8). This agrees with the observed increased spatial coherence there and supports a potential causative relation between irrigation and landslide occurrence. Still, the decreased spatial coherence observed in the entire AoI suggests that other processes may be linked to landslide presence.

4.2. Logistic Regression Results

To assess the relative contribution of natural and anthropogenic processes to landslide occurrence, a logistic regression analysis was performed. Eight models were analyzed with various combinations of the seven continuous variables and one categorical variable representing terrace presence (Table 1). Complete regression equations, coefficients and ROC curves are included in [50]. The models yielded an AUC of the receiver operator function between 64% and 88%, and AIC values between 284 and 348 (Table 4), representing poor to acceptable predicting performances. In every model, the inclusion of terrace presence as a variable slightly increased regression performance according to both AUC and AIC indicators. The most relevant variables were, in order of importance according to the coefficients obtained: distance to rivers, curvature, geology, and the presence of terraces. Slope and TWI, frequently important predictors of landslides [41,42,43,44,45,46] yielded high p-values in all cases, suggesting limited predictive power, most likely due to the averaging of highly variable values within polygons analyzed, which points to a limitation of the approach used. Due to this limitation, and being the main goal of this analysis testing the influence of terrace presence in predictive power, slope and TWI were not included in further models.

Model 1 was run using all the available variables except for terrace presence and yielded an AUC of 0.81 and an AIC of 289, and Model 2 was run adding terrace presence, and showed a slight improvement in AUC and AIC (Table 4). Variables representing slope, elevation, and TWI yielded p-values > 0.05 and were considered limited predictors of landslide presence, most likely due to the averaging of highly variable values within the polygons used.

Model 3 was conducted without the variables slope, elevation, and TWI, considered limited in prior models. All the variables used yielded p-values < 0.05 and were considered relevant terms of the regression. The regression equation for Model 3 can be seen in Equation (9). Model 4, was run including the variables tested in Model 3 plus terrace presence, producing the regression Equation (10).

P (1) = e x p (Y^{'}) / (1 + e x p e x p (Y^{'})) Y^{'} = - 1.04 + 0.512 G e o l o g y + 4.07 L o n g . C u r v a t u r e - 1.66 P e r p . C u r v a t u r e - 150 D i s t a n c e t o r i v e r s

(9)

\begin{matrix} Including Terraces variable P (1) = e x p (Y^{'}) / (1 + e x p e x p (Y^{'})) \\ {0 Y}^{'} = - 1.22 + 0.428 G e o l o g y + 4.07 L o n g . C u r v a t u r e - 1.81 P e r p . C u r v a t u r e - 138 D i s t a n c e t o r i v e r s \\ {1 Y}^{'} = - 0.49 + 0.428 G e o l o g y + 4.07 L o n g . C u r v a t u r e - 1.81 P e r p . C u r v a t u r e - 138 D i s t a n c e t o r i v e r s \end{matrix}

(10)

Despite being considered statistically significant according to p-values, further models were run without geology and curvature to assess changes in regression performance. Model 5 and 6, without geology, and Models 7 and 8, not including geology and curvature, showed a slight fall in performance (Table 4).

The results of Model 4 were further explored by calculating the probabilities of landslide presence associated with the variables measured for each landslide or control area. The multivariable regression fits were then plotted against the variables value used in Model 4 (Table 4) and separated by terrace presence. The plots are shown in Figure 9. Overall, trends towards higher landslide probabilities were found for observations linked to terrace presence for every variable explored in Model 4. These recognized trends sustain the idea that terraces are statistically more frequently associated with landslides than non-terraced areas, independent of natural variables like geology, curvature, and distance to fluvial erosion.

4.3. Finite Element Modeling Results

FEM modeling of the San Cristobal landslide was conducted for a three-layered model (Figure 7), as shown in Figure 10. Strength reduction factor analyses yielded FoS between 1.3 and 1, whereas several models did not converge. Groundwater depth was changed between 2703 m and 3055 m, yielding a FoS between 1 and 1.3 for the model (Figure 10 and Figure 11). The lowest groundwater depth that yielded a FoS of 1, pointing to transient instability of the slope, was 2885 m (Figure 10 and Figure 11). This value of groundwater depth is considered a threshold for slope stability based on the input materials used for modeling.

Most models showed a concentration of maximum shear strain close to the toe of the slope that propagated upwards when incorporating lower strength reduction factors. Strain seemed to concentrate at the boundary between the debris avalanche deposits and the harder Cretaceous sedimentary rocks (Figure 10). The obtained maximum shear strain zone is in good agreement with interpreted failure surfaces derived from satellite image observations (Figure 10).

The models yielded increased sensitivity to changes in friction angle, related to material strength, and to changes in groundwater depth, potentially related to irrigation, as shown in Figure 11. Cohesion, Young’s modulus, and Poisson ratio had limited sensitivity in affecting slope stability. Only in extreme cases did models show no convergence pointing to failure (Figure 11).

4.4. Seepage Analysis Results

Seepage analyses were conducted under changing permeability and infiltration conditions (results are shown in Figure 10 and Figure 12). In most cases, groundwater flow resulted in two piezometric surfaces, one perched within the fine top pyroclastic deposits and a main one within the debris avalanche deposits, as shown in Figure 10. Most scenarios tested yielded a FoS over 1.2, with models displaying a FoS of 1 under increased infiltration rates of 6 × 10⁻⁵ m/s, equivalent to 189 mm/yr, and below 260 mm/yr, the average annual rainfall of the region (Figure 12). Regarding changes in permeability, most models yielded a FoS > 1.2. Only when permeabilities were decreased by a minimum of a factor of ten did the models not converge, pointing to failure (Figure 12). In this case, failure is considered a consequence of decreased permeability and increased water saturation and decreased effective stresses promoting failure.

5. Discussion

This research was aimed at assessing the influence that irrigated terraces have on landslide occurrence in an area where communities are losing valuable croplands to landslide geohazards. The lack of sufficient knowledge of the processes promoting landslides prevents identification of effective mitigation measures to reduce hazards. The main hypothesis of this project was that terraces were promoting landslide occurrence in the area. Three methods were used to address the scientific questions formulated, two of them focusing on spatial statistical approaches, and a third one focusing on potential mechanisms behind slope failures in the area.

Overall, the presence of terraces appeared to be spatially correlated to landslide presence, with a correlation of 80% for areas surrounding the main rivers and communities in the AoI where irrigation and fluvial erosion were concentrated (Table 3, Figure 8). This result points to an increased exposure to landslide geohazards of valuable croplands and communities, concentrated along the Carumas and Putina rivers (Figure 8). When analyzing the entire AoI, including areas distant from the river, a poorer correlation of 63% was calculated, suggesting the existence of other natural processes accounting for landslides in areas where human activity is minimal or inexistent (Table 3).

To investigate the role of natural and anthropogenic variables, a logistic regression was performed. The results suggest that terrace presence also appears to improve landslide prediction, not as a governing factor but as a contributing factor in a system apparently dominated, in order of importance, by distance to rivers, high slope curvature, and weak geology (Table 4, Figure 9). In every pair of the models tested, the inclusion of terraces as an additional variable increased the performance indicators AUC and AIC (Table 4). This was explored in detail for Model 4 by plotting the data against natural variables but discriminating them by terrace presence or absence, and a clear trend was identified supporting increased landslide probability linked to terraces (Figure 9). We interpret that terraces decrease slope stability in a setting already preconditioned to slope instability through fluvial erosion, high relief, and poor geologic materials.

Both the spatial coherence and logistic regression analyses addressed the role of terraces in landslide occurrence through statistical methods. FEM numerical modeling of a failed slope was conducted to explore the geomechanical processes potentially linking terraces and landslides. Initially, the effect of terrace loading on the slope was considered an important factor to test. However, the scale of the terraces, up to 2 m in height (see Section 2, Figure 2, Figure 4, and Figure 7), was negligible in comparison to the scale of the slopes analyzed, which were hundreds of meters high and wide (Figure 7 and Figure 10). In addition, the fact that terraces were mainly composed of reworked slope material (see Section 2, Figure 4) meant that the additional loading imposed to slopes by terraces construction could be considered negligible. Thus, the main factor linking terraces and slope instability would be irrigation rather than loading.

The main modeling results point to both groundwater and friction angle as governing factors in slope stability (Figure 10, Figure 11 and Figure 12). The successful back-modeling of the failed slope supports the selected material properties, especially the most sensitive friction angle (Figure 11). Seepage models showed a plausible groundwater increase due to surface recharge at infiltration values of 189 mm/yr, which is below the average annual surface irrigation input of 260 mm/yr [23], with this seepage being one of the most sensitive factors leading to failure (Figure 10, Figure 11 and Figure 12).

In summary, terraces are considered spatially correlated to landslides and associated with increased probabilities of landslide presence, whereas irrigation excess water infiltrating from terraces appears to be a driving factor of slope instability. An in-depth discussion of the results is included in the next sections, along with some informed recommendations for communities in the AoI.

5.1. Spatial Coherence Analysis

The spatial coherence analysis conducted for the entire AoI showed a poor performance, yielding performance statistics values between 36% and 63%, when trying to predict landslide presence from irrigation features (Table 3). Separating irrigation features showed that landslides correlated with 10% more precision to in-slope terraces than to flat irrigated areas (Table 3). Still, the highest predictive performance was obtained when lumping all types of irrigation features together, achieving a precision of 63%, which can still be considered unconvincing in terms of predictive power and spatial coherence (Table 3).

Since most of the landslides and irrigation features recognized in the area appeared to be concentrated along the main rivers in the AoI, a focused spatial coherence analysis was conducted (Table 3, Figure 8). The results showed that correlation and model performance improved to performance statistics values between 64% and 81% (Table 3). Based on these results, it is concluded that close to the main rivers and communities, landslides are correlated to irrigation features in up to 81% of the cases (Table 3). This higher correlation surrounding the main rivers agrees with an increased exposure to fluvial erosion but also to increased irrigation close to the main communities where most of the crops are grown (Figure 8). Since correlation does not imply causation, further analyses were needed to assess the relative contribution of natural and anthropogenic factors to landslide presence in the area.

The main limitations associated with this methodology include the use of human mapped features potentially subject to observer biases. In addition, the selection of the buffer zone restricting the analysis can have an impact on the correlation.

5.2. Logistic Regression

According to the regression equations obtained for all the models (Equations (9) and (10); Table 4), the most relevant predictors of landslide presence are, in order of importance, distance to rivers, longitudinal curvature, perpendicular curvature, and geology. We interpret that decreasing the distance to rivers would correlate to higher landslide presence due to an increase in fluvial erosion, while longitudinal and perpendicular curvature indicate incised topography, a typical indicator of landslides. Finally, poorer geological categories indicate lower-strength materials that logically would be related to increased landslide occurrence. Results of the logistic regression analysis showed that not all the variables selected for the regression were relevant. Slope, elevation, and TWI did not pass the statistical tests, yielding p-values well above 0.05, although it is possible that the high variability of their raster values within the polygons used had masked such influence. For instance, slope and TWI have been considered relevant predictors of landslide presence in previous studies [41,42,43,44,45,46]. We argue that the averaging of highly variable predictors within landslide polygons masks the statistical relevance of such predictors, constituting a significant limitation of the approach used. Despite such limitations, it was still possible to test the change in predictive power by including terraces.

One of the main predictors of landslide presence was distance to rivers (Equations (9) and (10), Table 4). This variable was intended to be used as a proxy for fluvial erosion, a known promoter of landslide occurrence in fluvially dominated landscapes like the AoI. The high coefficients obtained for this variable (Equations (9) and (10)) suggest that other processes, apart from fluvial erosion, may also have been captured, such as poorer geologic conditions for unconsolidated sediments, lower depth to groundwater, and steeper slopes, among others. The potential contribution of different, and not fully accounted for, variables within the distance to rivers proxy constitutes an important limitation of this approach. Nevertheless, we conclude that, despite the influence of overlapping variables, distance to rivers captures the effects of fluvial erosion as a main predictor of landslide presence.

The best regression performance, considering the number of variables and AUC/AIC performance indicators, seems to be achieved by Models 3 to 6, with increased performance shown by Models 4 and 6 that included terraces (Table 4). The remaining models (Table 4) are considered either too complex (Models 1 and 2), or too simple (Models 7 and 8) and yielded decreased performance.

In Models 2, 4, 6, and 8, terrace presence was added to the regression as a categorical variable (Table 4), yielding a p-value below 0.05 in all cases and thus considered a significant variable in predicting landslide presence. The addition of terraces to the analyses slightly improved the regression in every model by approximately 1–2%, as shown by both AUCs and AICs indicators (Table 4). Such change points to a secondary influence in landslide prediction in comparison to natural factors like geology, curvature, and fluvial erosion. Yet, when discriminating calculated probabilities based on actual observation data for Model 4 (Figure 9), observations linked to terraces systematically yielded higher probabilities of landslide presence. It is concluded that, despite not being a primary driver or predictor of landslide presence, terraces increase the probability of landslide presence in the area.

5.3. Finite Element Numerical Modeling

The selection of material parameters based on the literature review, rather than based on subsurface sample testing, limited model accuracy. Nevertheless, the results of stability analyses appear to be robust, with most resulting in FoS values close to 1 (Figure 10). The Mohr–Coulomb constitutive model selected also seemed to be appropriate, yielding convergent results. The sensitivity analysis showed that variations in input parameters produced FoS results that were either far too stable or lacked convergence for most of the models tested.

In terms of limitations and error sources, some models exhibited limited variability towards the boundaries, suggesting that future analyses should model larger areas to avoid boundary effects.

Regarding the seepage analysis, the steady-state approach used is limited by an infinite time assumption leading to equilibrium. Still, by testing different infiltration rates, the influence that irrigation could have on groundwater levels and instability was tested with acceptable outputs (Figure 10 and Figure 11). Infiltration was shown to influence instability, promoting failure at values of 189 mm/yr (Figure 11), a conservative value that represents 70% of the 260 mm/yr annual rainwater [23] the area receives without taking into consideration artificial irrigation. Further data on irrigation magnitudes used in the area would help better constrain the models. In addition, decreases in permeability also led to sensitive increases in instability (Figure 11). Both the increase in infiltration and the decrease in permeability allowed for water build-up (Figure 11), increased saturation, and decrease in effective stresses promoting failure, all of which indicate the influence that groundwater has on slope instability.

Based on the conducted analysis, some recommendations are offered. Both the sensitivity and the seepage analysis support the need for further data on groundwater levels, material geomechanical properties, and irrigation measurements (Figure 10). According to Figure 11, stability appears to be more sensitive to changes in permeability than to changes in infiltration, suggesting that mitigation measures that increase percolation and water movement could have a larger effect on stability than decreasing the irrigation, which is crucial to sustain crops. It is thus recommended to consider the possibility of draining heavily irrigated slopes by installing pattern drains. These drains would probably have to extend at least 100 m into the slope to reduce saturation in the most likely slip surfaces for the San Cristobal landslide, as exemplified in Figure 10.

Finally, the geomechanical processes leading to slope failure have been explored. Previous studies have concluded that flood irrigation surpluses are the root cause of irrigation-induced landslides due to the saturation and rise in groundwater and decrease in effective stresses [7,13,15]. Slope stability models showed increased sensitivity to changes in groundwater levels (Figure 11), with rising levels linked to increased pore water pressures and decreased effective stresses. Seepage and stability models indicate that flood irrigation surplus may be percolating into susceptible slopes beyond the shallow terraces, thereby promoting landslide occurrence (Figure 10 and Figure 12). The adoption of water saving techniques to avoid slope saturation by using drip irrigation instead of flood irrigation has been proposed to reduce irrigation-induced landslide risk globally [1,13,15,51].

6. Conclusions

In this study, three methodologies were used to investigate the relation between irrigated agricultural terraces and landslide occurrence in the Ticsani valley, Moquegua region, Peru. Local communities there are exposed to increased landslide geohazards in an area characterized by steep slopes, fluvial erosion, and weak geologic materials, where intensive flood-irrigation has been used to support crop growth. Based on the results of this study and its discussion, we arrived at the following conclusions:

There is a spatial correlation between irrigated terraces and landslides surrounding the main rivers and communities of the Ticsani valley.
Irrigated terraces enhance landslide presence predictions, pointing to a contributing role in landslide occurrence.
Terraces do not likely change large-scale slope loading conditions, but they likely change groundwater conditions due to irrigation surplus.
Groundwater elevation and friction angle are considered the most sensitive factors influencing slope stability, according to performed modeling.
Vertical infiltration at conservative rates can raise groundwater levels to values causing slope failure, according to the models performed.
Rising groundwater due to irrigation surplus is considered the main anthropogenic cause behind landslide occurrence in the area.
Drainage of susceptible slopes and adoption of water saving techniques like drip-irrigation are recommended in increased risk areas surrounding communities.
Mitigation efforts should be concentrated in slopes surrounding main communities, underlain by poor geologic materials such as recent pyroclastic deposits, that receive increased infiltration water due to the presence of irrigated terraces upslope.

Author Contributions

Conceptualization, G.R. and P.S.; methodology, G.R.; validation, G.R., and P.S.; formal analysis, G.R.; resources, P.S.; data curation, G.R. and I.E.P.; field work, G.R., P.S., A.L.V.L. and C.J.B.P., writing—original draft preparation, G.R.; writing—review and editing, P.S. and I.E.P.; project administration, P.S.; funding acquisition, P.S. All authors have read and agreed to the published version of the manuscript.

Funding

Essential funding to carry out this research was provided by a Fulbright/Ministry of Education of Argentina scholarship granted to G.R. Colorado School of Mines provided tuition funds to G.R and provided research fellowship funds to I.E.P. through the First-Year Innovative & Research Scholar Training (FIRST) fellowship. Field work was funded by the Institute for Initiatives in Latin America, Colorado School of Mines.

Data Availability Statement

Supporting data and additional results can be found in [50].

Conflicts of Interest

The authors declare no competing interests.

References

Garcia-Chevesich, P.; Wei, X.; Ticona, J.; Martínez, G.; Zea, J.; García, V.; Alejo, F.; Zhang, Y.; Flamme, H.; Graber, A.; et al. The impact of agricultural irrigation on landslide triggering: A review from Chinese, English, and Spanish literature. Water 2020, 13, 10. [Google Scholar] [CrossRef]
Liu, Z.; Qiu, H.; Zhu, Y.; Huangfu, W.; Ye, B.; Wei, Y.; Tang, B.; Kamp, U. Increasing irrigation-triggered landslide activity caused by intensive farming in deserts on three continents. Int. J. Appl. Earth Obs. Geoinf. 2024, 134, 104242. [Google Scholar] [CrossRef]
Ballón Marroquín, A. Deslizamientos de tierras en el distrito de Aczo (Provincia de Huari, departamento de Ancash). In Comisión Carta Geológica Nacional; Compilación de Estudios Geológicos. Boletín N_13; Instituto Minero, Geológico y Metalúrgico: Lima, Peru, 1996; pp. 177–190. [Google Scholar]
Araujo Huamán, G.E.; Taipe Maquerhua, E.L.; Miranda Cruz, R.; Valderrama Murillo, P.A. Dinámica y Monitoreo del Deslizamiento de Siguas. Región Arequipa, Provincia Caylloma y Arequipa, Distrito Majes y San Juan de Siguas. 2017. Available online: https://hdl.handle.net/20.500.12544/791 (accessed on 11 November 2024).
Lacroix, P.; Dehecq, A.; Taipe, E. Irrigation-triggered landslides in a Peruvian desert caused by modern intensive farming. Nat. Geosci. 2020, 13, 56–60. [Google Scholar] [CrossRef]
Lacroix, P.; Handwerger, A.; Bievre, G. Life and death of slow-moving landslides. Nat. Rev. Earth Environ. 2020, 1, 404–419. [Google Scholar] [CrossRef]
Graber, A.; Santi, P.; Meza Arestegui, P. Constraining the critical groundwater conditions for initiation of large, irrigation-induced landslides, Siguas River Valley, Peru. Landslides 2021, 18, 3753–3767. [Google Scholar] [CrossRef]
Flamme, H.E.; Krahenbuhl, R.A.; Li, Y.; Dugan, B.; Shragge, J.; Graber, A.; Sirota, D.; Wilson, G.; Gonzales, E.; Ticona, J.; et al. Integrated geophysical investigation for understanding agriculturally induced landslides in southern Peru. Environ. Earth Sci. 2022, 81, 309. [Google Scholar] [CrossRef]
Lemus, O.; Santi, P.M.; Colque, P.; Meza, P.; Salas, G. Failure conditions and triggers of the Achoma landslide, central Andes region, Arequipa Peru. In 2023 Graduate Research and Discovery Symposium (GRADS) Posters and Presentations; Colorado School of Mines, Arthur Lakes Library: Golden, CO, USA, 2023. [Google Scholar]
Zhang, Z.; Zeng, R.; Zhao, S.; Meng, X.; Ma, J.; Yin, H.; Long, Z. Effects of irrigation projects on the classification of yellow river terrace landslides and their failure modes: A case study of heitai terrace. Remote Sens. 2023, 15, 5012. [Google Scholar] [CrossRef]
Peng, D.; Xu, Q.; Zhang, X.; Xing, H.; Zhang, S.; Kang, K.; Qi, X.; Ju, Y.; Zhao, K. Hydrological response of loess slopes with reference to widespread landslide events in the Heifangtai terrace, NW China. J. Asian Earth Sci. 2019, 171, 259–276. [Google Scholar] [CrossRef]
Wen, Y.; Gao, P.; Mu, X.; Li, M.; Su, Y.; Wang, H. Experimental Study on Landslides in Terraced Fields in the Chinese Loessial Region under Extreme Rainfall. Water 2021, 13, 270. [Google Scholar] [CrossRef]
Xu, Q.; Zhao, K.; Liu, F.; Peng, D.; Chen, W. Effects of land use on groundwater recharge of a loess terrace under long-term irrigation. Sci. Total Environ. 2021, 751, 142340. [Google Scholar] [CrossRef]
Zhang, F.; Wang, G.; Peng, J. Initiation and mobility of recurring loess flowslides on the Heifangtai irrigated terrace in China: Insights from hydrogeological conditions and liquefaction criteria. Eng. Geol. 2022, 302, 106619. [Google Scholar] [CrossRef]
Wang, J.; Liu, L.; Zhao, K.; Wen, Q. Farmers’ adoption intentions of water-saving agriculture under the risks of frequent irrigation-induced landslides. Clim. Risk Manag. 2023, 39, 100484. [Google Scholar] [CrossRef]
Deng, C.; Zhang, G.; Liu, Y.; Nie, X.; Li, Z.; Lui, J.; Zhu, D. Advantages and disadvantages of terracing: A comprehensive review. Int. Soil Water Conserv. Res. 2021, 9, 344–359. [Google Scholar] [CrossRef]
Gao, X.; Roder, G.; Ding, Y.; Liu, Z.; Tarolli, P. Farmers’ landslide risk perceptions and willingness for restoration and conservation of world heritage site of Honghe Hani Rice Terraces, China. Landslides 2020, 17, 1915–1924. [Google Scholar] [CrossRef]
Rozaki, Z.; Wijaya, O.; Rahmawati, N.; Rahayu, L. Farmers’ Disaster Mitigation Strategies in Indonesia. Rev. Agric. Sci. 2021, 9, 178–194. [Google Scholar] [CrossRef]
Salas, D. Andenes, Agrosistema Frágil; Llerena, C.A., Inbar, M., Benavides, M.A., Eds.; Conservación y Abandono de Andenes; Universidad Nacional Agraria La Molina: Lima, Peru, 2004; pp. 23–44. [Google Scholar]
Xu, Y.; Lu, Z.; Leshchinsky, B. Kinematics of irrigation-induced landslides in a Washington Desert: Impacts of basal geometry. J. Geophys. Res. Earth Surf. 2022, 127, e2021JF006355. [Google Scholar] [CrossRef]
Bradley, K.; Mallick, R.; Andikagumi, H.; Hubbard, J.; Meilianda, E.; Switzer, A.; Hill, E.M. Earthquake-triggered 2018 Palu Valley landslides enabled by wet rice cultivation. Nat. Geosci. 2019, 12, 935–939. [Google Scholar] [CrossRef]
Watkinson, I.M.; Hall, R. Impact of communal irrigation on the 2018 Palu earthquake-triggered landslides. Nat. Geosci. 2019, 12, 940–945. [Google Scholar] [CrossRef]
Huanca, E.; Marisol, S. Ciclos Horarios de Precipitación en el Perú Utilizando Información Satelital. 2016. Available online: https://hdl.handle.net//20.500.12542/112 (accessed on 12 November 2024).
Marroquín Liu, A.M. Balance Hídrico Superficial de la Subcuenca del Río Paltiture. Bachelor’s Thesis, Universidad de Piura, Piura, Peru, 2016. [Google Scholar]
de Moquegua, G.R. Plan de Competitividad Región Moquegua 2012–2021. 2012. Available online: https://consultas.regionmoquegua.gob.pe/wp-content/uploads/transparencia/InformacionAdicional/Plan%20Competitividad%20Region%20Moquegua%202012-2021%20GRM-UPSM.pdf (accessed on 12 November 2024).
Mariño Salazar, J.; Thouret, J.C. Geología, Historia Eruptiva y Evaluación de Peligros del Volcán TICSANI (Sur del Perú). 2003. Available online: https://repositorio.igp.gob.pe/items/3db71979-77c9-4bdc-9fef-26ba41b068ca (accessed on 12 November 2024).
Mariño Salazar, J.; Rivera Porras, M.A.; Cruz Pauccara, V.; Cacya Dueñas, L. Efectos Geológicos y en Las Construcciones Producidos Durante la Crisis Sísmica de Calacoa-Moquegua (Octubre 2005); Sociedad Geológica del Perú—SGP: Lima, Perú, 2006. [Google Scholar]
Cruden, D. A simple definition of a landslide. Bull. Eng. Geol. Environ. 1991, 43, 27–29. [Google Scholar] [CrossRef]
Comercio, E. Search of Periodistical Articles Tagged with the Word “Deslizamientos”. 2023. Available online: https://elcomercio.pe/noticias/deslizamientos/ (accessed on 12 November 2024).
Vílchez Mata, M.S. Casos Históricos de Movimientos en Masa Que Causaron Grandes Daños en Perú; Instituto Geológico, Minero y Metalúrgico—INGEMMET: Lima, Peru, 2018. [Google Scholar]
Vilímek, V.; Hanzlík, J.; Sládek, I.; Šandov, M.; Santillán, N. The share of landslides in the occurrence of natural hazards and the significance of El Niño in the Cordillera Blanca and Cordillera Negra Mountains, Peru. In Landslides: Global Risk Preparedness; Springer: Berlin/Heidelberg, Germany, 2013; pp. 133–148. [Google Scholar]
Margirier, A.; Audin, L.; Carcaillet, J.; Schwartz, S.; Benavente, C. Tectonic and climatic controls on the Chuquibamba landslide (western Andes, southern Peru). Earth Surf. Dyn. 2015, 3, 281–289. [Google Scholar] [CrossRef]
Quispevisana, L.; Zapata, A. Memoria Descriptiva de la Geología del Cuadrángulo de Omate (34-u); INGEMMET: Lima, Per, 2000. [Google Scholar]
INGEMMET Geocatmin Digital Repository. 2023. Available online: https://geocatmin.ingemmet.gob.pe/geocatmin/ (accessed on 12 November 2024).
Donkin, R.A. Agricultural Terracing in the Aboriginal New World; Viking Fund Publications in Antopology (56); University of Arizona Press: Tucson, AZ, USA, 1979. [Google Scholar]
Cruden, D.M.; Varnes, D.J. Landslides: Investigation and mitigation. Chapter 3-Landslide types and processes. Transp. Res. Board Spec. Rep. 1996, 247, 36–75. [Google Scholar]
Highland, L.M.; Bobrowsky, P. The Landslide Handbook—A Guide to Understanding Landslides (No. 1325); US Geological Survey: Reston, VA, USA, 2008. [Google Scholar]
Yamagishi, H.; Moncada, R. TXT-tool 1.081-3.1 landslide recognition and mapping using aerial photographs and Google Earth. In Landslide Dynamics: ISDR-ICL Landslide Interactive Teaching Tools: Volume 1: Fundamentals, Mapping and Monitoring; Springer International AG: Berlin/Heidelberg, Germany, 2018; pp. 67–82. [Google Scholar]
Comber, A.; Fisher, P.; Brunsdon, C.; Khmag, A. Spatial analysis of remote sensing image classification accuracy. Remote Sens. Environ. 2012, 127, 237–246. [Google Scholar] [CrossRef]
Sokolova, M.; Lapalme, G. A systematic analysis of performance measures for classification tasks. Inf. Process. Manag. 2009, 45, 427–437. [Google Scholar] [CrossRef]
Ayalew, L.; Yamagishi, H. The application of GIS-based logistic regression for landslide susceptibility mapping in the Kakuda-Yahiko Mountains, Central Japan. Geomorphology 2005, 65, 15–31. [Google Scholar] [CrossRef]
Lombardo, L.; Mai, P.M. Presenting logistic regression-based landslide susceptibility results. Eng. Geol. 2018, 244, 14–24. [Google Scholar] [CrossRef]
Regmi, N.R.; Giardino, J.R.; Vitek, J.D. Characteristics of landslides in western Colorado, USA. Landslides 2014, 11, 589–603. [Google Scholar] [CrossRef]
Southerland, L.; Zhou, W. Comparison of Two Logistic Regression Models for Landslide Susceptibility Analysis Through a Case Study. Environ. Eng. Geosci. 2021, 27, 471–486. [Google Scholar]
Crawford, M.M.; Dortch, J.M.; Koch, H.J.; Killen, A.A.; Zhu, J.; Zhu, Y.; Bryson, L.S.; Haneberg, W.C. Using landslide-inventory mapping for a combined bagged-trees and logistic-regression approach to determining landslide susceptibility in eastern Kentucky, USA. Q. J. Eng. Geol. Hydrogeol. 2021, 54, qjegh2020-177. [Google Scholar]
Killen, A. Using Logistic Regression to Generate and Compare Six Landslide Inventory-Based Susceptibility Maps in El Paso County, Colorado. Master’s Thesis, Colorado School of Mines, Golden, CO, USA, 2023. [Google Scholar]
Rocscience Inc. RS2 Finite Element Modeling Software.; Rocscience Inc.: Toronto, Canada, 2024. [Google Scholar]
Vezzoli, L.; Apuani, T.; Corazzato, C.; Uttini, A. Geological and geotechnical characterization of the debris avalanche and pyroclastic deposits of Cotopaxi Volcano (Ecuador). A contribute to instability-related hazard studies. J. Volcanol. Geotherm. Res. 2017; 332, 51–70. [Google Scholar]
Cervantes, J.; Martínez Valladares, W.; Romero Fernández, D.; Zapata Montes, A.A.; Navarro Colque, P.A. The Matalaque Formation of southern Peru: New stratigraphic and geochemical data. In Proceedings of the 6th International Symposium on Andean Geodynamics, Extended Abstracts, Barcelona, Spain, 12–14 September 2005. [Google Scholar]
Ronda, G. Linking Inca Terraces with Landslide Occurrence in the Moquegua Region, Peru. Master’s Thesis, Colorado School of Mines, Golden, CO, USA, 2024. [Google Scholar]
Bosco, G.; Dalpiaz, M.; Simeoni, L. Effects of modified irrigation procedures on the stability of cultivated slopes. In Landslides and Engineered Slopes. Experience, Theory and Practice; CRC Press: Boca Raton, FL, USA, 2018; pp. 483–489. [Google Scholar]

Figure 1. AoI of the project, the Ticsani valley located in the Moquegua region, Southern Peru. Ticsani is an active volcano. The red box outlined in the insert corresponds to the area of study, which includes the communities of Carumas, Sacuaya, San Cristóbal, and Challsahuaya.

Figure 2. (A) The main fluvial landforms in the AoI. The Carumas river can be seen here to be bounded by landslides affecting a fluvial terrace covered by agricultural terraces. (B) Landslides affecting a slope close to a community in the AoI. Note the construction of terraces on prior landslide deposits.

Figure 4. Typical profile of terraces built in Southern Peru. Wall rocks are directly piled over an excavated trench in bedrock without mortar. An upward decreasing gradation is used to fill the internal portion of the wall in compacted layers. Fertile arable soil is placed in the top 30 cm. Modified from [19].

Figure 5. (Left) Confusion matrix used to classify landslide and terraces polygons to quantify the spatial coherence between both. (Right) Examples showing the classification procedure conducted visually, where terraces are green polygons and landslides are yellow polygons (indicated by white arrows). True positives (top left) corresponded to landslides that had, at least partially, terraces upstream, up to a distance of 500 m. False positives (top right) corresponded to terraces or irrigated areas without any associated landslides downstream. False negatives (bottom) corresponded to landslides that did not have any associated terraces or that had terraces that were built over the landslides deposits at a later date.

Figure 6. Sketch showing the collection of variables to perform logistic regression. The analysis was performed for the whole AoI, only a portion is shown as an example in the sketch. Continuous variables were obtained by averaging raster values within polygons. Distance to drainage, representative of fluvial erosion potential, was calculated by nearest distance to polygons tools. Both the dependent variable representing landslide presence and the independent variable representing terrace presence were quantified by a binary categorical response (1 or 0).

Figure 7. (A) A 3D view of the San Cristobal landslide. (B) A 1:50,000 geological map of the landslide area, modified from [33]. Ki-mat = Matalaque Fm. P-Pi = Puno Fm. Q-pl = colluvial deposits, debris avalanche. Qp-vl-pi = pyroclastic deposits. Red line marks the topographic cross-section shown in A. (C) Cross-section of the San Cristobal landslide with estimated base groundwater conditions (blue) and materials boundaries (green) to be used for modeling and interpreted potential failure surfaces (dashed red).

Figure 8. Polygons used to conduct a confusion-matrix based analysis of areas surrounding main rivers and communities. Terraces and flat irrigated areas were integrated into one class for this analysis.

Figure 9. Multi-variable regression-calculated probability of landslide presence (fits) versus: (A) mean geology; (B) longitudinal curvature; (C) perpendicular curvature; (D) distance to rivers. Data are separated into observations with and without terraces. A quadratic regression was included, which shows higher probabilities of landslide presence when terraces are present, independently of the mean geology value.

Figure 10. (A) Topographic model of the San Cristobal landslide and interpreted potential slip surfaces. (B) Maximum shear strain marking potential slip surfaces in the base model performed. Note the spatial correlation between interpreted initial failure surfaces and the modeling results. (C) Seepage modeling results. Most models yielded a 2 piezometric surface (in purple) solution with a main aquifer within the debris avalanche deposits at 2700 m elevation.

Figure 11. Sensitivity analysis results. Increased sensitivity was obtained for friction angle (PHI) and groundwater depth (GW), whereas models were relatively insensitive to cohesion (C), Young’s modulus (E), and Poisson ratio (P).

Figure 12. Infiltration vs. permeability plot displaying the obtained FoS from the seepage analysis. Gray curves are lower limits for FoS fields.

Table 1. Variables and data sources used for logistic regression of landslides inventories in published research and this study.

Parameter	Description	Source	References
Elevation	Elevation (masl) from available DEM of the AoI	Alos Palsar DEM (JAXA 2013), 12.5 m from original 30 m resolution (ASF 2023)	[41] [42] [43] [45] [46]
Slope	First derivative of the elevation raster	ArcGIS PRO toolbox	[41] [42] [43] [44] [45] [46]
Curvature (axial and perpendicular)	Second derivative of the elevation raster	ArcGIS PRO toolbox	[42] [43] [45] [46]
Geology	Four-category raster of decreasing quality	Geological Map	[41] [43] [44] [46]
TWI	Topographic wetness index	$l n (\frac{u p s l o p e c o n t r i b u t i n g a r e a}{t a n t a n (s l o p e)})$	[42] [43] [46]
Fluvial erosion	Used as a proxy of fluvial erosion	Nearest distance to drainage ArcGIS PRO toolbox	[42] [43] (stream power index)

Table 2. FEM material properties. Debris avalanche and pyroclastic deposits properties were taken from [48]. Matalaque Fm. properties were taken from Swiss Standard SN 670.

Material	Unit Weight (kN/m³)	Friction Angle (°) Peak/Residual	Cohesion (MPa)	Tensile Strength (MPa)	Young´s Modulus (MPa)	Poisson Ratio	Permeability (m/s)
Debris avalanche (Q-Pl)	25.3	42.5/42	0	0	30	0.25	1 × 10⁻⁴
Pyroclastic deposits (Qp-vl-pi)	16.5	48/43	0	0	40	0.25	1 × 10⁻⁷
Matalaque Fm. (Ki-mat)	26.5	39/25	11.5/3.1	1	40	0.25	1 × 10⁻⁸

Table 3. Spatial correlation results for the entire AoI. Complete dataset included in [50]. Four scenarios were tested using different polygons as predictors including only terraces within the AoI, irrigated areas (<25° slope) within the AoI, terraces and irrigated areas within the AoI, and terraces and irrigated areas within the main rivers buffer (Figure 8).

AoI Analysis	SQRT Area Normalized Classes			Confusion Matrix Statistics
	TP	FN	FP	PR	RE	F1
Terraces	0.572	0.428	0.595	0.490	0.572	0.528
Irrigated	0.361	0.639	0.616	0.369	0.361	0.365
Terraces + irrigated	0.630	0.370	0.602	0.511	0.630	0.565
Terraces + irrigated (main rivers buffer)	0.809	0.191	0.446	0.644	0.809	0.717

Table 4. Logistic regression results for the models tested, including different variables. AUC is the area under curve, and AIC is Akaike’s information criterion, both statistics that reflect predictive power.

	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7	Model 8
Variables	Slope, elevation, curvature, TWI, geology, distance to rivers	Slope, elevation, curvature, TWI, geology, distance to rivers, terrace presence	Curvature, geology, distance to rivers	Curvature, geology, distance to rivers, terrace presence	Curvature, distance to rivers	Curvature, distance to rivers, terrace presence	Distance to rivers	Distance to rivers, terrace presence
AUC	0.81	0.82	0.80	0.81	0.78	0.80	0.64	0.69
AIC	289	284	292	288	297	291	348	339

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Ronda, G.; Santi, P.; Pope, I.E.; Vargas Luque, A.L.; Paria, C.J.B. Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru. Geosciences 2024, 14, 315. https://doi.org/10.3390/geosciences14110315

AMA Style

Ronda G, Santi P, Pope IE, Vargas Luque AL, Paria CJB. Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru. Geosciences. 2024; 14(11):315. https://doi.org/10.3390/geosciences14110315

Chicago/Turabian Style

Ronda, Gonzalo, Paul Santi, Isaac E. Pope, Arquímedes L. Vargas Luque, and Christ Jesus Barriga Paria. 2024. "Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru" Geosciences 14, no. 11: 315. https://doi.org/10.3390/geosciences14110315

APA Style

Ronda, G., Santi, P., Pope, I. E., Vargas Luque, A. L., & Paria, C. J. B. (2024). Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru. Geosciences, 14(11), 315. https://doi.org/10.3390/geosciences14110315

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Linking Inca Terraces with Landslide Occurrence in the Ticsani Valley, Peru

Abstract

1. Introduction

2. Background

2.1. Geomorphology

2.2. Geology, Volcanism, and Seismicity

2.3. Landslide Geohazards

2.4. Terraces and Irrigation Techniques

3. Methods

3.1. Landslide Terraces Spatial Coherence Assessment

3.2. Logistic Regression Analysis of the Landslide Inventory

3.3. Finite Element Modeling of the San Cristobal Landslide

3.4. Seepage Analysis

4. Results

4.1. Spatial Coherence Analysis Results

4.2. Logistic Regression Results

4.3. Finite Element Modeling Results

4.4. Seepage Analysis Results

5. Discussion

5.1. Spatial Coherence Analysis

5.2. Logistic Regression

5.3. Finite Element Numerical Modeling

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI