Introduction

Agriculture is the main water-consuming industry in the world and accounts for 70% of the global water withdrawals from ground and surface freshwater resources (United Nations, 2021). In arid and semi-arid regions, withdrawals for agriculture exceed 90% of ground and surface water usage due to the exclusive use of irrigation (Frenken and Gillet 2012). Irrigation is a way to combat desertification in many countries that simultaneously increase food security (Williams 2011; Zhang and Huisingh 2018). However, irrigation in these regions can adversely affect the landscape by triggering landslides (Garcia-Chevesich et al. 2021).

Intense irrigation can cause a rise in the water table level and increase the water saturation of the soil (Lee et al. 2008). This results in higher hydrostatic pressure at the bottom of the slope, and some cases can result in soil liquefaction and deformation of the slope (Yan et al. 2019). Irrigation-triggered landslides occur in many different lithologies and soil types such as loess terrain, sedimentary and alluvial deposits, and volcanic materials (Garcia-Chevesich et al. 2021).

Irrigation-induced landslides have been studied worldwide for decades. Over the past few decades, the Chinese Loess Plateau (Heifangtai, Gansu province) has experienced thousands of landslides close to irrigated areas (Gu et al. 2020; Guo et al. 2015; Xu et al. 2012). The 2018 Palu landslide in Indonesia, which was originally thought to have been initiated due to an earthquake, was later found to have been caused by the excessive irrigation of rice crops (Bradley et al. 2019; Cummins 2019; Watkinson and Hall 2019). Irrigation-induced landslides have also been documented in California due to the irrigation of avocado fields (Knott 2008). Historical slides, like an ancient landslide in British Columbia, have also been attributed to the effects of irrigation (Clague and Evans 2003).

One area that has been highly affected by irrigation-triggered landslides is the agricultural Majes area in the Siguas Valley in southern Peru (Fig. 1). The Majes area is a critically important for the production of agricultural goods in the Arequipa region, and contributes over US-$800 million to the Peruvian GDP per year (Efrain and Churata 2018). Development for agriculture started along the northwest side of the Siguas River Valley in the early 1980s, and irrigation started in 1983. In total, the development has ~ 160 km2 of irrigated land in an arid environment (Majes I) (Wei et al. 2021). Expansion of the irrigation into Majes II could add 460 km2 of irrigated land in the southeastern part of the Siguas River Valley (Zeballos Arivilca 2020). While irrigation has brought benefits, it has also affected the local hydrogeology and is contributing to active landslide development along the Siguas River (Araujo et al. 2016); similar irrigation-induced landslides also affect other regions worldwide (Garcia-Chevesich et al. 2021).

Fig. 1
figure 1

Location maps of the Majes area in Peru. A Regional map of Peru showing the location of our study area; B Map of Majes region showing its agricultural areas. The purple box represents the area of our model domain; C zoom in to the landslide area (marked in red dashed rectangular) which also shows the proximity of the landslide to local irrigation-supported agriculture

Full size image

One of the most impactful landslides affecting the region is El Zarzal. The first water seepage from the river valley walls in the area was observed in 1996. Seepage at the lower toe of the valley walls started in 2004, and the failure that created the El Zarzal landslide started in 2005. Since then, retrogressive failure has continued to expand the landslide size and its impacts on the region. Graber et al. (2021) used landslide modeling to conclude that irrigation water in Majes I is raising the local water table and helping to drive El Zarzal’s continued failure.

While prior studies assumed water table change (Graber et al. 2021), we developed a groundwater modeling study of the Majes region and the surrounding non-irrigated land to better understand the relationship between irrigation since 1983 and landslide evolution. We use water table interpretations and the irrigation history data to constrain the hydrogeological properties of the sediments in the region. Then, using our calibrated flow models, we examined whether groundwater pumping could be used to lower the elevation of the water table in the area, and to reclaim water that could be reused for irrigation. We focused on the area near Majes I, and specifically on potential mitigation measures near El Zarzal. The management approaches outlined in this paper could be used to prevent future failure in areas like Majes II where simultaneous irrigation and pumping could minimize water table rise.

Methods

Groundwater flow numerical model

We built a 3D groundwater flow model domain of the Majes area and the El Zarzal landslide using DEM data with a resolution of 12 m. The model was built using FEFLOW (Diersch 2014), which uses the finite-element method to solve the groundwater flow equation. The model simulates the initial state of the water table level before irrigation started, after 40 years of irrigation, and after pumping near the landslide head scarp. The 3D model included 20 horizontal layers that created a model domain with dimensions of 14 km in length, 7 km in width, and a depth below the land surface of 138–618 m (Fig. 2). The model mesh contains 1,164,135 nodes and 2,202,480 triangular prismatic elements.

Fig. 2
figure 2

3D groundwater flow model domain and hydrogeological layers. The colors represent the hydraulic conductivity values of the different layers. The model domain location is in Fig. 1b. The yellow markers represent the location of the pumping wells simulated for this study

Full size image

Based on the stratigraphic description of the Majes geology (Flamme et al. 2022), our model included four hydrogeological layers with different hydraulic conductivities. The upper layer (the Millo Formation) is a conglomerate with a thickness of 20 m. Beneath the Millo Formation lie the Upper and Lower Moquegua Formations with thicknesses estimated to be about 120 m and 80 m, respectively. The Upper Moquegua Formation consists of sandstone and limonitic gravel, whereas the Lower Moquegua Formation is made up of claystone and siltstone. The stratigraphic layering of the model was fit to the geological cross-section next to the El Zarzal landslide and was extrapolated through the model domain. The Millo Formation hydraulic conductivity was set at 50 m d−1, a representative value for conglomerates. The Upper Moquegua Formation, which is the second hydrogeological layer in the model, was investigated with different hydraulic conductivity values. We focused on variability here as the water table location and most of the flow dynamics occur within this formation. The hydraulic conductivity of the 3rd layer (the Lower Moquegua Formation) was set with a value of 0.0001 m d−1 which is typical for clay, while the bottom layer was set with a hydraulic conductivity value of 0.00001 m d−1 which is typical for basement rocks (Fig. 2). The porosity used in the model is 0.3, which is within the sandstone porosity range.

In our simulations, we used a hydraulic head boundary condition on the north-eastern and south-western sides of the model, coordinated with the natural flow path of the region. The only hydrological information available is the interpretation of the water table depth of 80 m based on interpretations from transient electromagnetics (TEM) data (Flamme et al. 2022). Considering the assumptions we put on the hydraulic conductivity of the different layers, and because of the lack of hydrological information, our calibration process was by trial and error, changing the hydraulic head boundary condition such that the water table was at ~ 80 m depth, consistent with geophysical interpretations and seepage observations (Flamme et al. 2022). During our calibration process, we assumed an anisotropy ratio (Kz/Kx) of 0.5. Our approach resulted in a regional hydraulic gradient \(\left(\frac{dh}{dx}\right)\) of 0.0168. A hydraulic head boundary condition was also set for the riverbed with the head being set equal to the riverbed elevation to allow local groundwater flow to the river. The bottom, the north-western, and the southeastern sides of the model were set as no-flow boundaries.

Sensitivity analysis

To examine the sensitivity of the model to irrigation, the hydraulic conductivity of the 2nd hydrogeological layer (the Upper Moquegua Formation) was changed in multiple scenarios, and we evaluated the irrigation rate required to raise the water table by 40 m over 40 years, consistent with the irrigation history for Majes I. We also evaluated the water table response for this layer with the only water input being precipitation (17 mm year−1) (Wei et al. 2021) and no irrigation. After an irrigation-induced water table rise of 40 m was achieved, different pumping rates were tested to identify the rates required for reducing the water table back to a pre-irrigation level near the cliff’s edge where landside activity is prevalent. That was done to examine the option of mitigating the landslide and recovering groundwater for irrigation. Nine boreholes with depths of 70–150 m were simulated near the head scarp of the landslide (Fig. 2). The 3D approach enables an accurate calculation of the water table drawdown due to pumping. Another set of simulations was conducted in transient mode testing the response time for a water table rise of 10 m, and thereafter, the water table drop of 10 m due to pumping of 1 million cubic meters per year (MCMY). These scenarios were conducted for future water management in other places in the Majes region, presenting a sustainable option for irrigation sources without causing landslides by managing water table rise. However, reclaimed water could require treatment before reuse.

Results and discussion

Water table height optimization

Water table height analysis in the groundwater model was first conducted to match the interpreted water table next to the El Zarzal landslide (Flamme et al. 2022). Flamme et al. (2022) estimated the top of the water table through DC electrical resistivity and time-domain electromagnetic surveys next to the landslide and on the other side of the Siguas Valley which has not been irrigated and has a low water table. We interpret that the water table level of the non-irrigated region reflects the regional water table before irrigation started in Majes. A simple approach for model calibration is 1D where the water table rose from 80 m depth to 40 m depth over 40 years, which we use as the baseline for the simulations in this study.

Steady-state modeling results demonstrate that for the water table to be raised by 40 m, an irrigation rate of ~ 1400 mm year−1 is required. Intense irrigation in Majes started in 1983 and it continues until this day. Therefore, to understand if the water table level today is at a steady state, the same scenario was run in transient mode for 40 years and the results were compared to the steady-state results. We found less than 0.5% error between the water table level predicted the steady-state and the transient simulates, which implies that the hydrological system near El Zarzal is in a steady state with irrigation.

This irrigation rate value, however, is valid only for a specific range of hydraulic conductivities. Figure 3 shows the steady-state modeling results of the irrigation rate that is needed to raise the water table by 40 m in 40 years with different hydraulic conductivities. The mid-range values that were tested (10–3–1 m d−1) show similar results with an average of 1377 mm year−1 and an error of less than 3%. The lowest (10–4 m d−1) and highest (10 m d−1) hydraulic conductivity values show ~ 13% and ~ 50% error from the mid-values average, respectively. Based on this, we used a representative value of 0.01 m d−1 for simulations with pumping.

Fig. 3
figure 3

Relation between the Upper Moquegua Formation (2nd model layer) hydraulic conductivity and the irrigation rate required to raise the pre-irrigation water table by 40 m

Full size image

To fully understand our model sensitivity, mid-range hydraulic conductivity values were explored to investigate the impact of only precipitation (17 mm year−1) and no irrigation on the water table within the Upper Moquegua Formation (Fig. 4). As the hydraulic conductivity value increases, water table elevation increases. It can be noted, however, that the difference between the water table elevation with the lowest hydraulic conductivity value (10–3 m d−1) compared to that with the highest value (10–1 m d−1) is only ~ 1.5 m. Therefore, we conclude that the primary reason for the 40 m water table rise is irrigation.

Fig. 4
figure 4

Water table elevation at steady state with no irrigation for different hydraulic conductivities of the Upper Moquegua Formation (2nd model layer)

Full size image

Water table rise increases pore pressure and decreases frictional resistance, which can trigger landslides (Bogaard and Greco 2016). In heavily irrigated areas where landslides may be triggered, water table management may be an effective mitigation strategy. One way to prevent water table rise due to irrigation is the use of drip irrigation systems that provide the exact amount of water needed for each crop at different stages of growth (Hanson and May 2004; Smajstrla et al. 2000). In areas with more traditional irrigation methods, pumping wells in strategic locations can capture the percolated water and lower the water table. With this pumping approach, the water table can be lowered or managed, the likelihood of landslides onset can be decreased, and the pumped water can be reused which will lower the stress of the traditional water sources. Extracting less water from the surface reservoir (e.g., rivers, lakes) will be beneficial for the water body ecosystem and regional water sustainability. Re-use of the pumped groundwater could be problematic as the percolated water quality deteriorates with time. Salinization of the groundwater beneath an agricultural field may occur (Pulido-Bosch et al. 2018) and enrichment in other contaminants such as uranium may hinder the use of the water for agriculture (Jurgens et al. 2010). Therefore, the pumped water may need designated treatment such as ultra/nano-filtration or desalination depends on the specific contaminant before reuse for irrigation.

Pumping optimization for mitigation of further land failure

We evaluated how pumping from nine wells along the head scarp of the landslide could be used to manage the local water table. Given the limited hydrogeological data that are available for the area, we investigated a range of hydraulic conductivity values (0.001–0.1 m d−1) of the Upper Moquegua Formation to estimate the pumping rate that would be needed to lower the water table by 40 m along the cliff.

We specifically focused on the pumping rate that would be needed to reduce the water table elevation by 40 m to counter the 40 m rise of the water table that was caused by irrigation, such as interpreted at Majes I (Fig. 5). Lower hydraulic conductivities required lower pumping rates to reduce the water table and higher hydraulic conductivities required higher pumping rates to achieve the same water table reduction. This trend was expected as pumping from a low conductivity layer creates a larger drawdown cone than pumping from a higher conductivity layer which allows for the capture of water from a larger area with less drawdown. That happens due to the larger hydraulic gradients required to drive flow in lower conductivity layers.

Fig. 5
figure 5

Pumping rate that is needed to lower the water table by 40 m near the cliff edge for different hydraulic conductivities of the Upper Moquegua Formation (2nd model layer)

Full size image

This pressure difference between the pumping well and its surroundings is critical to understand and evaluate for each location where pressure reduction through groundwater pumping is suggested as a mean for mitigating landslide onset. Figure 6 shows the simulated pressure distribution after irrigation but before pumping in the upper model layer and 150 m below the land surface where pumping wells screens are producing water. Figure 7 shows the simulated pressure distribution on the surface and 150 m below the surface after pumping to change the pressure. The differences between the figures show the pressure reduction next to the El Zarzal landslide and the potential mitigation of it. In the case where K = 0.01 m d−1 for the Upper Moquegua layer, a pumping rate of 0.36 MCMY would be needed for the water table elevation to be lowered by 40 m. Using a cross-section of the pressure distribution, we can also see how the pressure changes near the stream, the cliff, and the headward of the scarp (Fig. 8), which clearly shows that the water table was lowered near the pumping zone and near the active head scarp.

Fig. 6
figure 6

Pressure maps of the model domain after irrigation. a 3D representation of the topography (the upper slice of the model); b 2D representation of the black square in a, showing the landslide area. The white line represents the water table location; c 3D representation of the subsurface 150 m below the topography; d 2D representation of the landslide area 150 m below the topography (black square in c)

Full size image
Fig. 7
figure 7

Pressure maps of the model domain after pumping. a 3D representation of the topography (the upper slice of the model); b 2D representation of the black square in a, showing the landslide area. The white line represents the water table location, and the red dots represent the pumping well locations.; c 3D representation of the subsurface 150 m below the topography; d 2D representation of the landslide area 150 m below the topography (black square in c)

Full size image
Fig. 8
figure 8

Pressure cross-sections after irrigation (a) and after pumping (b). c The location of the cross-section which is also shown in Fig. 7b (black line). The black line (in a and b) represents the water table level. The vertical black lines in (b) represent the pumping well location

Full size image

Lowering the height of the water table by 40 m along the cliff results in a pressure reduction of ~ 200 kPa at the elevation of the pre-irrigation water table. Despite the reduction in pore pressure, in some instances, such as in the El Zarzal landslide, this will not stop the landslide propagation but might slow its rate (Graber et al. 2021). Xu et al. (2012) showed that landslides may result in the formation of a concave topographic shape, resulting in increased groundwater accumulation and therefore a higher water table level. Therefore, it is more likely that a secondary landslide will occur there. Thus, greater attention must be paid to lowering the pressure within the El Zarzal landslide than in other nearby areas.

While pumping groundwater along the head scarp of the landslide will reduce the pore pressure and help slow down the landslide propagation, it could possibly provide another water source for irrigation, depending on the quality and treatment of the reclaimed water. This methodology of pumping groundwater for irrigation near a failing cliff could be adopted in relevant areas where the water quality for irrigation is sufficient. Furthermore, it may act as a preventive measure for the onset of a future landslide in some areas.

Pumping as a preventive measure for future landslides

To fully understand the benefit of pumping from a well field to prevent landslides in a newly irrigated area, we modeled several scenarios for a water table rise of 10 m due to irrigation and then, pumping groundwater to lower the water table by 10 m. That simulates a scenario where water table rise is detected before any landslide initiation, and pumping is employed to prevent any landslide activity. Given the uncertainty about the aquifer properties, we tested different hydraulic conductivities of the Upper Moquegua layer and pumped 1 MCMY from 9 wells along the cliff. The response time of the water table rise due to irrigation \({(\tau }_{i})\) and drop due to pumping \({(\tau }_{p})\) was monitored for each scenario.

Figure 9 shows the results of the water table response time to irrigation and pumping for our threshold of 10 m of change. The hydraulic conductivity value changes the response time of the water table for constant irrigation and pumping. When the hydraulic conductivity is low, the response times due to irrigation and pumping are longer than when hydraulic conductivity is higher. Even though the total pumping rate (1 MCMY) is significantly lower than the total irrigation rate (80 MCMY = 1400 mm year−1 over the entire irrigation area of 57 km2), it effectively lowers the water table near the cliff face. In addition, the ratio between the water table response time due to pumping and irrigation (τp τi−1) increases as hydraulic conductivity increases (Fig. 9 inset). From this, we interpret that as hydraulic conductivity increases, the pumping effect becomes more dominant than the irrigation effect on the water table response time. Moreover, this means that the water table in aquifers with higher hydraulic conductivity will be easier to moderate and there will be a lower risk of landslides if the water table rise is detected early.

Fig. 9
figure 9

The water table response time for raising it by 10 m due to irrigation (τi, blue squares), and the water table response time to lower it 10 m back due to pumping of 1 MCMY (τp, red circles) for different hydraulic conductivities of the 2nd layer. The inset represents the ratio of τp to τi

Full size image

While we focused on the Majes region of southern Peru, this approach could also be applied to other arid-to-semi-arid environments that have heavy irrigation and landslide risks. Our generalized model approach for the Majes region shows promise, however, it is limited due to the lack of water table and hydrogeological parameter. To improve the model application, which could lead to more precise information on number and location of pumping wells, monitoring the water table in newly developed irrigation regions would be beneficial. Furthermore, we propose that groundwater flow modeling approach such as presented in this study would be highly beneficial before developing agricultural fields next to cliffs and that local regulators should consider requiring such work before irrigation. With that, data on hydraulic conductivity of the different geological strata in the area, along with groundwater heads at multiple locations in any study region must be obtained for sufficient model calibration and validation. Therefore, monitoring wells should be situated in strategic locations for groundwater head monitoring and conducting pumping tests for obtaining hydrological information prior to model development.

Conclusions

We built a numerical 3D groundwater flow model of the Majes region around the El Zarzal landslide in southern Peru. We simulated the historic water table depth and its rise due to irrigation to estimate the hydraulic conductivity of the system where the water table is changing and the El Zarzal landslide has been active. After constraining hydraulic conductivity from historical data, we simulated different pumping approaches for lowering the water table for landslide mitigation. We show that pumping in the vicinity of the landslide lowers the pore pressure along the cliff which may reduce the propagation velocity of the landslide. In unirrigated areas where agriculture is planned, placing pumping wells along the cliff may significantly decrease landslide risk. Several simulations show that if the water table rises due to intense irrigation, pumping can successfully reduce the water table in a shorter time than the water table rise due to irrigation. As an example, in lower hydraulic conductivity (10–3 m d−1), the time it takes for the water table to rise 10 m due to irrigation is ~ 1000 days, while the time it takes for the water table to be lowered by 10 m due to pumping is ~ 250 days. We also demonstrate that formations with higher hydraulic conductivities have a faster response time of the water table level to pumping and thus, are easier to manage. Water recovered from pumping can be reapplied as irrigation water (pending appropriate treatment) which will reduce stress on the traditional water sources (e.g., rivers, lakes) and maintain their cultural and ecological importance in the area.