Abstract
In arid regions, mudflows are triggered by extreme precipitation in hydrologically sensitive basins. In the district of Coronel Gregorio Albarracín Lanchipa (Tacna, Peru), these processes mainly affect areas adjacent to ravines and steep slopes. This study delineates zones susceptible to flooding and mudflows using sequential hydrological–hydraulic modelling. Hydrometeorological records, high-resolution topography (0.1 m), and soil and land-cover data were integrated. Simulations were performed using HEC-HMS and HEC-RAS, and peak discharges were estimated for return periods (RP) of 25, 50, 100, and 200 years. Hydrographs were analyzed to assess basin response under extreme rainfall conditions. Model performance was evaluated through event-based validation using the March 2020 mudflow event. Estimated peak discharges were 11.9, 17.3, 25.3, and 122.4 m³/s for RP25, RP50, RP100, and RP200, respectively. The RP200 value is consistent with a historical event (~200 m³/s) reported in the 1920s. Results indicate a marked increase between RP100 and RP200, suggesting a nonlinear basin response under extreme conditions. Validation yielded an F1-score of 0.92, indicating strong agreement between observed and simulated extents. Simulations identified impacts on two population centers, one educational facility, 1.5 km of roads, and 44.6 ha of agricultural land, highlighting the importance of incorporating extreme scenarios in risk management and land-use planning.
1 Introduction
Mudflows are flow-type mass movements composed of highly concentrated mixtures of water and fine-grained sediments, whose rheological behavior ranges from viscous fluids to materials exhibiting significant yield stress (O’Brien et al., 1993; Vallejo, 1979). These processes represent one of the most destructive natural hazards in arid and mountainous regions due to their rapid propagation, high erosive capacity, and severe impacts on infrastructure, agricultural land, and human settlements (Gaprindashvili et al., 2021; Khamidullaev et al., 2025; Tsereteli et al., 2022a, 2022b). In arid and semi-arid environments, mudflows are commonly triggered by short-duration, high-intensity rainfall events that generate infiltration-excess (Hortonian) runoff in steep catchments with sparse vegetation cover (Goodrich et al., 1997; Highland, 2008; Iverson, 1997). Under these conditions, relatively small increases in rainfall intensity or duration can produce disproportionately large increases in runoff generation and sediment transport, resulting in highly nonlinear hydrological responses and the rapid mobilization of dense flows (Huang et al., 2020; Hungr et al., 2013; Reid et al., 2016).
Over recent decades, numerical modeling of mudflows and debris flows has advanced substantially through the incorporation of non-Newtonian rheological formulations (RAMSS, FLO 2D, HEC-RAS) and coupled hydrological–hydraulic modeling frameworks (Bruno et al., 2022; El-Bagoury, 2024; Flores et al., 2026). Rheological models such as Bingham and Herschel–Bulkley provide more realistic representations of sediment–water mixtures than conventional Newtonian approaches, particularly for hyperconcentrated or partially cohesive flows (O’Brien et al., 1993). Recent studies have demonstrated that rheological selection significantly influences flow depth, velocity, and runout distance, highlighting the importance of incorporating non-Newtonian behavior into hazard assessments (Hdeib et al., 2018; Lima et al., 2025). In parallel, coupled HEC-HMS (USACE-HEC, 2025a) and HEC-RAS (USACE-HEC, 2025b) approaches have gained increasing relevance for linking runoff generation with the hydraulic propagation of dense sediment-laden flows, especially in data-scarce catchments (Lee et al., 2025; Peker et al., 2024; Wei et al., 2024). Nevertheless, integrated hydrological–hydraulic modeling incorporating non-Newtonian formulations remains limited in the arid basins of southern Peru.
In southern Peru, mudflows and debris flows have repeatedly affected urban areas developed over alluvial fans and sectors adjacent to ephemeral channels (Del Savio et al., 2019). In the district of Coronel Gregorio Albarracín Lanchipa (CGALD) in Tacna, Peru, local reports and technical studies document recurrent events associated with short-duration, high-intensity storms despite the predominantly arid climatic conditions of the region (INGEMMET, 2021). This problem is further exacerbated by urban expansion into areas showing geomorphological evidence of previous flow activity and by the limited availability of continuous hydrometeorological records. Although recent studies in Peru have applied hydrological and hydraulic tools to evaluate hazards associated with floods and mass movements (García, 2016; Goyburo et al., 2025; Millán-Arancibia and Lavado-Casimiro, 2023), assessments integrating hydrometeorological data, high-resolution topography, and non-Newtonian hydraulic modeling to delineate mudflow-prone areas in arid basins of southern Peru remain scarce.
A particularly underexplored aspect in these catchments is the nonlinear hydrological response to extreme rainfall events. Unlike humid environments, where subsurface flow processes dominate, arid basins are primarily controlled by Hortonian runoff, producing highly nonlinear relationships between rainfall intensity, storm duration, and peak discharge (Goodrich et al., 1997). Consequently, relatively small storms may generate limited hydrological responses, whereas events exceeding critical intensity–duration thresholds can trigger abrupt increases in runoff and sediment entrainment, ultimately favoring mudflow initiation (Reid et al., 2016). Despite its relevance for hazard assessment, this nonlinear response remains poorly documented in Tacna and, more broadly, in arid basins of southern Peru.
Within this context, the present study evaluates mudflow susceptibility in the CGALD using a coupled hydrological–hydraulic modeling framework based on HEC-HMS and HEC-RAS, incorporating a non-Newtonian Bingham rheological formulation. Using hydrometeorological records, high-resolution topography, and basin-scale physical and environmental characterization, peak discharges associated with different return periods were estimated and subsequently used to simulate mudflow propagation and inundation patterns. This study additionally provides one of the first integrated non-Newtonian mudflow hazard assessments in arid basins of southern Peru, highlighting the role of nonlinear hydrological response in the generation and propagation of extreme sediment-laden flow events.
2 Materials and methods
2.1 Study area
The study area is located in the district of Coronel Gregorio Albarracín Lanchipa (CGALD), within the province and department of Tacna, southern Peru (Figure 1). The hydrological basin draining through CGALD toward the coastal plain covers approximately 162 km2 and was delineated using the 30 m spatial resolution FABDEM digital elevation model (Hawker et al., 2022). The basin forms part of the Caplina River system, which originates in the western Andes cordillera and discharges into the Pacific Ocean. The basin exhibits a pronounced altitudinal gradient, ranging from 400 to 1,045 m a.s.l. (INGEMMET, 2021), favoring short concentration times and rapid hydrological responses to intense rainfall events. Geomorphologically, the area is composed of coastal hills, piedmont slopes, and alluvial fans dissected by ephemeral channels draining toward the coastal plain. The principal drainage axis corresponds to the Caplina River.
Figure 1
Geological assessments conducted by INGEMMET identify extensive sectors underlain by unconsolidated Quaternary alluvial deposits (Qal unit), characterized by low cohesion and high erodibility, and classified as having “high” to “very high” susceptibility to debris flows and mudflows under intense precipitation conditions (INGEMMET, 2021). These geomorphological conditions are coupled with an arid to semi-arid climate (BWh/BSh), where mean annual precipitation ranges from 20 to 50 mm in lowland sectors to 150–250 mm in the upper basin (Aybar et al., 2020). More than 80% of annual precipitation occurs during the austral summer (January–March), typically as short-duration, high-intensity convective storms. Owing to the low infiltration capacity of desert soils, these events commonly generate infiltration-excess (Hortonian) runoff and abrupt hydrological responses (Goodrich et al., 1997; Reid et al., 2016).
Since its establishment in 2001, CGALD has undergone rapid and largely unplanned urban expansion, currently exceeding 15,000 inhabitants (INEI, 2024). Urban growth has progressively encroached onto alluvial fans and sectors adjacent to ephemeral channels, increasing exposure to hydromorphological hazards (Gebremicael et al., 2022; Reid et al., 2016). Several debris-flow and mudflow events have been documented through satellite imagery and technical reports. Among these, a notable event occurred in March 2020, when extreme rainfall (~30 mm in 24 h) triggered a mudflow affecting areas downstream of the Arunta Bridge and nearby quarry sectors. This event provides representative evidence of the basin’s susceptibility to short-duration, high-intensity storms and was used as a reference event for evaluating the consistency of the hydraulic modeling conducted in this study.
2.2 Method and data preprocessing
Figure 2 summarizes the methodological framework adopted in this study, which is organized into three coupled components: (i) hydrological modeling, (ii) two-dimensional hydraulic modeling, and (iii) hazard scenario generation. The approach integrates hydrometeorological forcing with the simulation of mudflow propagation in arid basins characterized by limited data availability.
Figure 2
The hydrological component began with basin delineation and morphometric characterization using the 30 m spatial resolution FABDEM digital elevation model (Hawker et al., 2022). Geomorphological and hydrological parameters, including slope, drainage network, and concentration time, were derived from this dataset. Rainfall–runoff transformation was performed using HEC-HMS v4.12 (USACE-HEC, 2025b), applying the SCS-CN method for loss estimation and a unit hydrograph approach for discharge generation.
The hydraulic component was developed using HEC-RAS v6.5 (USACE-HEC, 2025b) in a two-dimensional (2D) configuration, using the generated hydrographs as upstream boundary conditions. High-resolution topography was obtained through RPAS photogrammetry, enabling detailed representation of the microtopography controlling flow propagation across the channel and alluvial fan. Flow behavior was represented using a non-Newtonian Bingham rheological formulation, with rheological parameters adopted from the literature for unconsolidated alluvial materials comparable to those present in the study area (Bilal and Kumar, 2025; Moreiras et al., 2021; Ortega et al., 2022).
Hydraulic model validation was conducted through spatial comparison between the simulated inundation extent and the observed footprint of the March 2020 event, using the F1-score as the performance metric. Finally, model outputs were integrated to generate flow depth and velocity maps, delineate potentially affected areas, and classify hazard levels based on hydrodynamic criteria. The following subsections describe each methodological component in detail.
2.2.1 Digital elevation model for hydrology and hydraulic modeling
In this study, two digital elevation models (DEMs) with different spatial resolutions were employed according to the scale requirements of hydrological and hydraulic modeling. For basin-scale hydrological characterization, the FABDEM product with 30 m spatial resolution was used (Hawker et al., 2022). This bare-earth DEM is suitable for delineating drainage networks and estimating geomorphological parameters such as slope, contributing area, and time of concentration, which are controlled by mesoscale topographic features (Lindsay, 2016). For hydraulic modeling, a high-resolution DEM (0.17 × 0.17 m) was generated using RPAS photogrammetry, covering approximately 6.1 km2 along the main channel and alluvial fan (Figure 3). This fine-scale dataset captures microtopographic features such as channel banks, terraces, and anthropogenic modifications (quarry excavations), which exert strong control on flow routing, depth, and lateral spreading in mudflow events (Bates et al., 2010). The high-resolution DEM was integrated into the two-dimensional (2D) HEC-RAS model, where it defines the computational domain and spatial discretization of the flow equations. The selected resolution is consistent with hydraulic modeling practices that require terrain data capable of resolving dominant flow features to minimize numerical diffusion and improve simulation accuracy (Neal et al., 2012). Despite its high resolution, the RPAS-derived DEM is subject to photogrammetric uncertainties, including vertical errors and surface noise, which may influence hydraulic results, particularly in low-gradient areas. These limitations are considered in the interpretation of model outputs.
Figure 3
2.2.2 Hydrometeorological data
The study area is monitored by pluviometric stations operated by National Service of Meteorology and Hydrology of Peru (SENAMHI) and National Water Authority of Peru (ANA), with records dating back to 1960 and 1990, respectively, complemented by local stations installed for this study (Table 1). The network covers an altitudinal gradient from 58 m a.s.l. (La Yarada) to 3,100 m a.s.l. (Palca). All records were homogenized to daily scale. The Jorge Basadre station (560 m a.s.l.) was selected as the primary reference gauge because: (i) its 31-year continuous record (1993–2024) (OMM, 2025); (ii) proximity to the CGALD (<1 km); and (iii) representativeness of precipitation conditions in the CGALD. To extend the series to the 1981–2023 period, two gap-filling methods were applied: (i) simple linear regression with nearby stations (Calana and Magollo, calibration period 1993–2014); and (ii) simple linear regression with the PISCO gridded product (Aybar Camacho et al., 2017). Coefficients were estimated using ordinary least squares. Validation combined leave-one-out cross-validation and an independent period (2015–2016), with performance evaluated using R2 and RMSE. The procedure was verified not to introduce bias in extreme values. From the resulting continuous series, the Annual Maximum Series (AMS) of 24-h precipitation was constructed. For each year of the precipitation series, the maximum daily value was extracted. Years with more than 10% missing daily data were excluded from the analysis. The resulting series (1981–2023, 43 years) forms the basis for extreme event frequency analysis and design storm derivation for different return periods.
Table 1
| Station | Latitude (°) | Longitude (°) | Altitude (masl) | Type | Source | Period | Temporal resolution |
|---|---|---|---|---|---|---|---|
| La Yarada | −18.21 | −70.52 | 58 | Conventional | Senamhi | 2017–2024 | Daily |
| Jorge Basadre | −18.03 | −70.25 | 560 | Conventional | Senamhi | 1993–2024 | Daily |
| Palca | −17.77 | −69.97 | 3,100 | Conventional | Senamhi | 1965–2014 | Daily |
| Calientes | −17.88 | −70.14 | 1,200 | Conventional | Senamhi | 2018–2024 | Daily |
| Calana | −17.94 | −70.18 | 848 | Conventional | Senamhi | 1963–2024 | Daily |
| Magollo | −18.12 | −70.33 | 288 | Conventional | Senamhi | 1965–1995 | Daily |
Precipitation stations used for hydrological modeling in the CGALD (Tacna, Peru).
Figure 4a shows the comparison of daily precipitation data between the Jorge Basadre rain gauge (red) and the corresponding PISCO gridded cell (blue) at the same location. The time series reveals gaps in the historical record during the periods 1981–1993 and 2014–2016, which were filled using a combined approach based on nearby stations, daily PISCO data, and linear regression. Figure 4b presents the scatter plot between the PISCO product and the Jorge Basadre precipitation station. The results indicate strong agreement between both datasets, with correlation coefficients of 0.92 at the monthly scale and 0.97 at the daily scale, demonstrating a high level of consistency and an adequate fit for gap-filling purposes.
Figure 4
Figure 5 presents the annual time series of maximum 24-h precipitation (PPmax24h) for the period 1981–2023. The bars represent the observed annual maxima, and the red line indicates the fitted linear trend. The series shows pronounced interannual variability, with predominantly low to moderate values during the 1980s and early 1990s (<10 mm), followed by higher-magnitude events from the mid-1990s onward. Notable peaks occur in 1997 (≈10 mm) and especially in 2020 (≈26 mm), the latter representing the maximum value of the entire record. Subsequently, values return to moderate levels during 2021–2023. The positive linear trend suggests a gradual increase in the intensity of extreme precipitation events over the study period; however, this tendency is strongly influenced by the occurrence of isolated extreme episodes. Overall, the series reflects non-stationary behavior characterized by sporadic high-intensity events that increase variability and contribute to the observed upward trend in PPmax24h.
Figure 5
2.2.3 Hydrological modeling
Peak discharges associated with different return periods (RP) were estimated using HEC-HMS v4.12 (USACE-HEC, 2025a), configured under a semi-distributed framework. Rainfall–runoff transformation was performed using the SCS-CN loss method, which is appropriate for infiltration-excess runoff conditions characteristic of arid basins. Curve Number (CN) values were assigned based on land use and hydrologic soil groups derived from GIS analysis and thematic mapping. Runoff generation was represented using the SCS unit hydrograph, whereas flow routing was simulated using the kinematic wave method (Meselhe et al., 2004; Orellana, 2021; Salazar-Briones et al., 2018).
The concentration time (Tc) was estimated from geomorphological parameters derived from the 30 m DEM, including slope, flow length, and drainage area. Owing to the absence of continuous discharge records, conventional hydrological calibration was not feasible. Consequently, model parameters were constrained within physically plausible ranges reported for comparable arid basins (Cardich Motta, 2017; Roque, 2022).
Hydrological consistency was evaluated through indirect validation. Simulated peak discharges for the RP50 scenario were compared with the maximum available records at the Caplina gauging station, showing agreement in order of magnitude. Additionally, the RP50 hydrograph used as upstream boundary condition in the hydraulic model reproduced an inundation extent spatially consistent with the March 2020 event. Although this procedure does not constitute a formal calibration, it provides a reasonable assessment of model plausibility under data-scarce conditions.
Design storms were derived through precipitation frequency analysis for return periods of 25, 50, 100, and 200 years. The resulting hydrographs were subsequently used as upstream boundary conditions in the hydraulic model.
2.2.4 Hydraulic modeling
The propagation of sediment-laden flows was simulated using the two-dimensional (2D) unsteady flow module of HEC-RAS version 6.5 (USACE-HEC, 2025b), employing its integrated non-Newtonian formulation based on a Bingham rheological model (Vallejo, 1979). This approach represents the flow as an equivalent continuum characterized by a finite yield stress (τᵧ) and a viscous component (η), which is suitable for dense, sediment-rich flows commonly observed in arid alluvial environments (INGEMMET, 2021; Pino, 2013; Vilcanqui-Alarcón et al., 2022).
The model solves the depth-averaged Saint-Venant equations for mass and momentum conservation. Within this framework, flow resistance is internally modified to account for non-Newtonian behavior through the specification of τᵧ and η(Equation 1). A conceptual expression of the modified friction slope for a Bingham fluid can be written as:Where h is flow depth (m), U is the velocity magnitude (m·s−1), ρ is fluid density (kg·m−3), and g is gravitational acceleration (9.81 m·s−2). This expression is provided for interpretative purposes, as the full solution is handled internally by the HEC-RAS numerical solver.
The computational domain was discretized using a high-resolution Digital Elevation Model (0.17 m spatial resolution) derived from RPAS photogrammetry. This dataset captures detailed channel geometry, overbank areas, and anthropogenic features (e.g., quarry excavations), which exert strong control on flow routing and lateral spreading across the alluvial fan.
Upstream boundary conditions were defined using hydrographs generated by the hydrological model for return periods of 25, 50, 100, and 200 years, while downstream conditions were specified as normal depth based on local channel slope. Simulations were performed under unsteady flow conditions to reproduce the temporal evolution of the flow wave.
The rheological parameters required for the non-Newtonian Bingham model in HEC-RAS were estimated based on specialized literature and previous studies conducted in similar geomorphological settings (Goyburo et al., 2025; INGEMMET, 2021; Keaton et al., 2019; Pino, 2013; USACE-HEC, 2025b). Rheological parameters (Table 2) were estimated using empirical relationships (Equation 2) that express τᵧ and η as exponential functions of the volumetric sediment concentration (Cᵥ), following established formulations (Iverson, 1997; O’Brien et al., 1988):Where α₁, β₁, α₂, and β₂ are empirical coefficients dependent on sediment characteristics. A representative range of Cᵥ = 0.30–0.35 was adopted, consistent with dense sediment-laden flows in arid basins, where steep slopes, short-duration high-intensity rainfall, and abundant unconsolidated deposits promote high sediment availability. In the absence of site-specific rheometric data, these parameters should be interpreted as first-order approximations.
Table 2
| Parameter | Symbol | Value | Units |
|---|---|---|---|
| Volumetric sediment concentration | Cv | 0.30–0.35 | – |
| Plastic viscosity | η | 22.0 | Pa·s |
| Yield stress | τy | 210 | Pa |
| Specific gravity | Gs | 2.65 | – |
Rheological parameters used as input in the HEC-RAS non-Newtonian (Bingham) model for mudflow simulation in the CGALD (Tacna, Peru).
The modeling approach assumes a single-phase equivalent fluid and does not explicitly represent processes such as grain-size segregation, pore-pressure dynamics, or phase separation (Hungr et al., 2013). Consequently, model outputs are sensitive to the selected rheological parameters, particularly Cᵥ, which exerts a primary control on both τᵧ and η.
Model performance was evaluated through indirect, event-based validation by comparing the simulated inundation extent for the RP50 scenario with the observed footprint of the March 2020 event derived from Sentinel-2 imagery. The spatial agreement supports the plausibility of the adopted parameterization, although it does not constitute a formal calibration.
3 Results
3.1 Hydrological modeling
3.1.1 Construction of design hyetographs
The goodness-of-fit of five probability distributions (Normal, Log-Normal, Pearson Type III, Log-Pearson Type III, and Gumbel) was evaluated against the hydrological data series using the Kolmogorov–Smirnov test, with a significance level (α) of 0.05 and a sample size (n) of 43. The results, based on the maximum absolute deviation statistic (Dmax), indicate that only the Log-Normal distribution does not exceed the critical value (Dcritical = 0.20323) and therefore provides an adequate fit to the observed data. The Log-Normal distribution exhibits the best performance (Dmax = 0.19926), while the Gumbel (Dmax = 0.22476), Normal (Dmax = 0.23598), Pearson Type III (Dmax = 0.93903), and Log-Pearson Type III (Dmax = 0.9494) distributions exceed the critical threshold, leading to their rejection for modeling the analyzed series (Figure 6). Accordingly, the Log-Normal distribution is recommended for subsequent frequency analyses. The fitted distribution functions further confirm that the Log-Normal model provides the best representation of the maximum precipitation series at the Jorge Basadre station (Table 3).
Figure 6
Table 3
| Statistics | Normal | Log-Normal | Pearson lll | Log Pearson lll | Gumbel | ||
|---|---|---|---|---|---|---|---|
| n | 43 | Dmax | 0.23598 | 0.19926 | 0.93903 | 0.9494 | 0.22476 |
| a | 0.05 | Dcritical > Dmax | Does not fit | Fits | Does not fit | Does not fit | Does not fit |
| Dcritical | 0.20323 | Best fit | 3 | 1 | 4 | 5 | 2 |
Kolmogorov–Smirnov goodness-of-fit test.
Figure 7 presents the estimated maximum annual 24-h precipitation for different return periods at the Jorge Basadre station. Based on the observed data, projected precipitation depths were derived for various return periods using several probabilistic distributions. As the return period increases, the estimated precipitation values rise progressively, reaching up to 72.5 mm for a return period of 1,000 years. This highlights the potential occurrence of low-frequency but high-impact extreme rainfall events. Following the application of the Kolmogorov–Smirnov goodness-of-fit test, the Log-Normal distribution was identified as the model that most accurately reproduces the statistical behavior of the observed data. Therefore, it is considered the most suitable distribution for design and planning purposes related to extreme hydrological scenarios.
Figure 7
Table 4 presents the estimated design precipitation mudflow depths for different storm durations of less than 24 h at a rain gauge representative of the study area. A progressive increase in precipitation intensity is observed as both storm duration and return period increase, ranging from values close to 11.5 mm for 10-min events with a 2-year return period to more than 158 mm for 24-h storms associated with a 1,000-year return period. This pattern reflects the direct relationship between event duration and expected precipitation magnitude, underscoring the importance of probabilistic analysis for the design of hydraulic infrastructure and the evaluation of hydrological risk scenarios.
Table 4
| Duration | Return period (years) | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Hr | min | 2 | 5 | 10 | 20 | 25 | 50 | 100 | 200 | 300 | 500 | 1,000 |
| 0.17 | 10.00 | 0.74 | 1.40 | 2.11 | 3.09 | 3.47 | 4.96 | 7.00 | 9.79 | 11.89 | 15.14 | 20.94 |
| 0.33 | 20.00 | 0.88 | 1.67 | 2.51 | 3.67 | 4.13 | 5.90 | 8.32 | 11.65 | 14.14 | 18.00 | 24.90 |
| 0.50 | 30.00 | 0.97 | 1.85 | 2.78 | 4.06 | 4.57 | 6.52 | 9.21 | 12.89 | 15.65 | 19.92 | 27.55 |
| 0.67 | 40.00 | 1.04 | 1.99 | 2.99 | 4.37 | 4.91 | 7.01 | 9.89 | 13.85 | 16.81 | 21.41 | 29.61 |
| 0.83 | 50.00 | 1.10 | 2.10 | 3.16 | 4.62 | 5.19 | 7.41 | 10.46 | 14.65 | 17.78 | 22.64 | 31.31 |
| 1.00 | 60.00 | 1.16 | 2.20 | 3.31 | 4.83 | 5.44 | 7.76 | 10.95 | 15.33 | 18.61 | 23.69 | 32.77 |
| 1.50 | 90.00 | 1.28 | 2.43 | 3.66 | 5.35 | 6.02 | 8.59 | 12.12 | 16.96 | 20.59 | 26.22 | 36.26 |
| 2.00 | 120.00 | 1.37 | 2.61 | 3.93 | 5.75 | 6.47 | 9.23 | 13.02 | 18.23 | 22.13 | 28.18 | 38.97 |
| 4.00 | 240.00 | 1.63 | 3.11 | 4.68 | 6.84 | 7.69 | 10.97 | 15.49 | 21.68 | 26.31 | 33.51 | 46.34 |
| 6.00 | 360.00 | 1.81 | 3.44 | 5.18 | 7.56 | 8.51 | 12.14 | 17.14 | 23.99 | 29.12 | 37.08 | 51.28 |
| 7.00 | 420.00 | 1.88 | 3.57 | 5.38 | 7.86 | 8.84 | 12.62 | 17.81 | 24.93 | 30.26 | 38.54 | 53.30 |
| 8.00 | 480.00 | 1.94 | 3.70 | 5.56 | 8.13 | 9.14 | 13.05 | 18.42 | 25.78 | 31.29 | 39.85 | 55.11 |
| 10.00 | 600.00 | 2.05 | 3.91 | 5.88 | 8.59 | 9.67 | 13.80 | 19.47 | 27.26 | 33.09 | 42.13 | 58.27 |
| 11.00 | 660.00 | 2.10 | 4.00 | 6.03 | 8.80 | 9.90 | 14.13 | 19.94 | 27.92 | 33.89 | 43.15 | 59.68 |
| 12.00 | 720.00 | 2.15 | 4.09 | 6.16 | 9.00 | 10.12 | 14.44 | 20.38 | 28.53 | 34.63 | 44.10 | 60.99 |
| 24.00 | 1440.00 | 2.56 | 4.86 | 7.32 | 10.70 | 12.03 | 17.17 | 24.24 | 33.93 | 41.18 | 52.44 | 72.53 |
Estimation of sub-daily precipitation for different durations using the Dick-Peschke method.
Table 5 presents the design precipitation intensities associated with short-duration storms, estimated for different return periods at a representative rain gauge within the study area. The results show a progressive decrease in intensity as storm duration increases, ranging from values exceeding 120 mm/h for 10-min intervals with a 10-year return period to intensities on the order of 2–3 mm/h for 24-h events. Likewise, for any given duration, precipitation intensity increases with the return period, reflecting the more extreme and less frequent nature of high-recurrence events. This information is fundamental for hydraulic design and urban drainage analysis, as it enables the estimation of expected peak discharges and the establishment of design criteria for high-intensity precipitation events occurring over short time intervals.
Table 5
| Duration | Return period (years) | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Hr | min | 2 | 5 | 10 | 20 | 25 | 50 | 100 | 200 | 500 | 1,000 |
| 0.17 | 10.00 | 4.43 | 8.43 | 12.69 | 18.53 | 20.84 | 29.75 | 41.98 | 58.76 | 90.83 | 125.62 |
| 0.33 | 20.00 | 2.63 | 5.01 | 7.54 | 11.02 | 12.39 | 17.69 | 24.96 | 54.01 | 54.01 | 74.69 |
| 0.50 | 30.00 | 1.94 | 3.70 | 5.56 | 8.13 | 9.14 | 13.05 | 18.42 | 39.85 | 39.85 | 55.11 |
| 0.67 | 40.00 | 1.57 | 2.98 | 4.48 | 6.55 | 7.37 | 10.52 | 14.84 | 32.11 | 32.11 | 44.41 |
| 0.83 | 50.00 | 1.32 | 2.52 | 3.79 | 5.54 | 6.23 | 8.90 | 12.55 | 27.16 | 27.16 | 37.57 |
| 1.00 | 60.00 | 1.16 | 2.20 | 3.31 | 4.83 | 5.44 | 7.76 | 10.95 | 23.69 | 23.69 | 32.77 |
| 1.50 | 90.00 | 0.85 | 1.62 | 2.44 | 3.57 | 4.01 | 5.72 | 8.08 | 17.48 | 17.48 | 24.18 |
| 2.00 | 120.00 | 0.69 | 1.31 | 1.97 | 2.87 | 3.23 | 4.61 | 6.51 | 14.09 | 14.09 | 19.48 |
| 4.00 | 240.00 | 0.41 | 0.78 | 1.17 | 1.71 | 1.92 | 2.74 | 3.87 | 8.38 | 8.38 | 11.59 |
| 6.00 | 360.00 | 0.30 | 0.57 | 0.86 | 1.26 | 1.42 | 2.02 | 2.86 | 6.18 | 6.18 | 8.55 |
| 7.00 | 420.00 | 0.27 | 0.51 | 0.77 | 1.12 | 1.26 | 1.80 | 2.54 | 5.51 | 5.51 | 7.61 |
| 8.00 | 480.00 | 0.24 | 0.46 | 0.70 | 1.02 | 1.14 | 1.63 | 2.30 | 4.98 | 4.98 | 6.89 |
| 10.00 | 600.00 | 0.21 | 0.39 | 0.59 | 0.86 | 0.97 | 1.38 | 1.95 | 4.21 | 4.21 | 5.83 |
| 11.00 | 660.00 | 0.19 | 0.36 | 0.55 | 0.80 | 0.90 | 1.28 | 1.81 | 3.92 | 3.92 | 5.43 |
| 12.00 | 720.00 | 0.18 | 0.34 | 0.51 | 0.75 | 0.84 | 1.20 | 1.70 | 3.67 | 3.67 | 5.08 |
| 24.00 | 1440.00 | 0.11 | 0.20 | 0.31 | 0.45 | 0.50 | 0.72 | 1.01 | 1.41 | 2.19 | 3.02 |
Estimated design intensities for sub-daily precipitation durations.
Figure 8 presents the Intensity–Duration–Frequency (IDF) curves, which constitute a fundamental tool for characterizing the precipitation response to events of different magnitudes and recurrence levels. The results show that precipitation intensity decreases as storm duration increases, ranging from values exceeding 100 mm/h for 10-min events with return periods between 5 and 10 years to intensities of approximately 2–4 mm/h for 24-h storms. Conversely, for a given duration, intensity increases with the return period, reflecting the greater severity associated with low-frequency events, such as those with 500- or 1,000-year recurrences. This typical IDF structure enables the estimation of design precipitation for hydraulic infrastructure, urban drainage systems, and hydrological modeling applications. Consequently, it provides a key basis for defining expected peak discharges and safety margins in planning and risk management under extreme precipitation scenarios.
Figure 8
Figure 9 shows the sub-daily precipitation intensities (mm/h) for different return periods (RP25 to RP200), derived using the Dick-Peschke method from records at the Jorge Basadre station. The results reveal a marked concentration of the storm core around 780 min (13 h), where the absolute maximum intensity reaches 16.21 mm/h for RP200, representing a + 209% increase relative to RP100 and a + 409% increase compared to RP25. Before the peak, intensities increase progressively and moderately (e.g., at 600 min, RP200 records 1.13 mm/h versus 0.37 mm/h at RP25); after the peak, a sharp decline is observed, with values dropping below 2 mm/h from 840 min onward even for RP200. This pattern reflects the structure of the SCS Type II 12-h hyetograph, where approximately 80% of the precipitation is concentrated within a 3-h window (between 660 and 840 min), consistent with high-intensity convective storms typical of arid regions. The nonlinear response to increasing return periods, particularly the sharp jump between RP100 and RP200, suggests that once an intensity threshold is exceeded, the risk of mudflow generation in the study area increases dramatically.
Figure 9
3.1.2 Simulated hydrographs for different return periods
Table 6 summarizes the peak discharges derived from the synthetic hydrographs generated for different return periods, allowing the interpretation of the basin’s hydrological response to events of increasing severity. For a 25-year return period (RP25), a peak discharge of approximately 11.9 m3/s is obtained, which increases progressively to 17.3 m3/s for RP50 and 25.3 m3/s for RP100, reflecting greater runoff concentration and flow velocity under more intense precipitation events. This represents a 45% increase from RP25 to RP50 and a 113% increase from RP25 to RP100, indicating a nonlinear yet gradual amplification of the hydrological response as return periods lengthen. In contrast, the value associated with RP200 (122.4 m3/s) exhibits a marked increase in the hydrograph peak, 929% higher than the RP25 baseline, representing an extreme scenario in which the natural or structural conveyance capacity of the channel system may be exceeded. This sharp escalation underscores a potential threshold effect in the basin’s response, where the combination of precipitation intensity and antecedent saturation conditions leads to disproportionately high runoff generation. This result is particularly relevant for the design of hydraulic infrastructure and for flood and mudflow risk assessments (Figure 10).
Table 6
| Statistic | RP25 | RP50 | RP100 | RP200 |
|---|---|---|---|---|
| Peak Flow (m3/s) | 11.9 | 17.3 | 25.3 | 122.4 |
| Peak Flow Increase vs. RP 25 | – | +45% | +113% | +929% |
Peak flow estimation using the SCS unit hydrograph method (HEC-HMS) for different return periods based on precipitation data from the Jorge Basadre station.
Figure 10
3.2 Hydraulic modeling
3.2.1 Model validation
Figure 11 presents a comparison between the observed mudflow extent recorded in 2020, delineated from satellite imagery, and the simulated flood footprint corresponding to a 50-year return period (RP50). The spatial analysis indicates a high level of agreement between both extents, with an estimated overlap of approximately 95% of the affected area. The spatial overlay shows particularly strong correspondence along the main channel and in downstream lateral expansion sectors. Spatial agreement was further quantified using areal overlap metrics and the F1-score, yielding a value of 0.92, which indicates excellent correspondence between observed and simulated extents and supports the robustness of the adopted model parameterization. The main discrepancies are concentrated within the quarry area, where recurrent excavation and material removal activities have produced geomorphological changes after the documented event. These anthropogenic modifications likely altered the local microtopography and, consequently, the preferential flow paths.
Figure 11
3.2.2 Hydraulic modeling of different return periods
The hydrographs for the different return periods (RPs), previously derived in the maximum-flood hydrological analysis, were used as upstream boundary conditions for the hydraulic modeling. Based on these inputs, the fluvial system response was simulated under scenarios of increasing magnitude; the resulting mudflow depths and inundation extents are presented in Figure 12. The comparative analysis of the modeled scenarios reveals a progressive increase in both inundated area and mudflow depth as the return period increases. For the least severe scenario (Figure 12a, RP25), mudflow depths along the main channel are predominantly shallow, generally below 0.25 m. However, localized sectors with intensified hydraulic behavior are identified, most notably at the Arunta Bridge, where mudflow depths of approximately 1.5 m are reached. Likewise, the quarry area, whose morphology has been altered by extractive activities, acts as a zone of flow retention and concentration, exhibiting mudflow depths between 1.5 and 3.0 m, with a maximum of about 4.25 m at its outlet. In the RP50 scenario (Figure 12b), a moderate intensification of flooding is observed. Mudflow depths in the main channel increase to around 0.5 m, while values at the Arunta Bridge reach approximately 1.75 m. Within the quarry, mudflow depths range from 2.0 to 3.25 m and increase to about 4.5 m at the outlet, confirming its role as a critical zone of hydraulic accumulation. The RP100 scenario (Figure 12c) exhibits a more severe hydraulic response. Main-channel mudflow depths approach 1.0 m, while mudflow depths at the bridge increase to roughly 2.0 m. In the quarry sector, mudflow depths vary between 2.25 and 3.5 m, with a maximum near 4.75 m at the outlet, indicating a sustained increase in flow energy and stored volume. Finally, the most extreme scenario (RP200; Figure 12d) represents a critical flooding condition. In this case, mudflow depths in the main channel locally exceed 4.0 m, reaching values close to 5.0 m at the bridge. The most severe impacts are again concentrated in the quarry area, where the combination of high discharges and strongly modified geomorphology produces exceptional mudflow depths, modeled between 4.0 and 5.0 m. The maximum value, on the order of 5.5 m, is again recorded at the quarry outlet, identifying this sector as the zone of highest hydraulic hazard within the study area.
Figure 12
The hydraulic modeling reveals a well-defined spatial pattern in flow velocity distribution for the analyzed return periods (Figure 13). In general, velocities in both the main channel and adjacent floodplain remain within moderate ranges, with most values below 0.5 m/s across the modeled domain. However, the results consistently identify a localized zone of flow acceleration present in all simulated scenarios. This acceleration zone is concentrated in the mid-reach of the quarry, where the combined effects of channel geometry, depressions generated by extractive activities, and flow confinement promote a significant increase in hydraulic energy. The comparison among scenarios indicates a progressive rise in maximum velocities as event magnitude increases. In the lower-magnitude scenarios, peak velocities reach approximately 2.25 m/s (Figure 13a) and 2.5 m/s (Figure 13b), increasing to around 2.75 m/s in the intermediate scenario (Figure 13c). The most extreme case, associated with the highest return period (Figure 13d), exhibits markedly more critical behavior, with maximum velocities reaching up to 5.5 m/s in the same sector. Outside this localized zone, flow velocities remain within low to moderate ranges, and no other acceleration cores of comparable magnitude are identified along the analyzed reach. This pattern underscores the decisive role of local geomorphological alterations in intensifying flow dynamics and confirms the quarry sector as the area of highest hydraulic hazard in terms of erosion potential, sediment transport capacity, and potential damage to nearby infrastructure.
Figure 13
3.2.3 Mudflow hazard map for historical events
Figure 14a illustrates the spatial distribution of mudflow depth associated with the simulated historical extreme flood event (200 m3/s discharge). Mudflow depths are classified into four ranges (≤1 m, 1–2 m, 2–3 m, and >3 m), with the highest values concentrated along the main channel and in topographically depressed areas. Zones with mudflow depths exceeding 2 m appear discontinuously, following paleo-flow paths and natural flow expansion sectors, which reflects the strong control exerted by local geomorphology. In contrast, areas located farther from the channel axis predominantly exhibit shallower mudflow depths, although localized impacts are still identified within urbanized zones (parks, residential areas, among others). This spatial distribution enables the identification of sectors with potential for significant structural damage, particularly in built-up areas located near the active channel margins. Figure 14b depicts the flow velocities reached during the same extreme event, classified into intervals of ≤1 m/s, 1–2 m/s, 2–3 m/s, and >3 m/s. The highest velocities are mainly concentrated along the main channel and at bends, where geomorphological confinement promotes higher hydraulic energy. In lateral overbank areas, velocities decrease progressively, with low to moderate values predominating in more open urban sectors. Nevertheless, localized high-velocity corridors are observed crossing populated areas, indicating a considerable potential for erosion and structural washout in zones adjacent to housing. The superposition of the local road network and populated areas makes it possible to identify infrastructure exposed not only to inundation but also to dynamically hazardous flow conditions. The simulated flood extent initiates notably upstream of the Arunta hill bend; within this affected corridor, approximately 300 houses belonging to the Alfonso Ugarte and San Antonio population centers, the educational institution I. E. 42255 Santa Teresita del Niño Jesús, about 1.5 km of roads, and 44.6 ha of agricultural land are directly impacted.
Figure 14
4 Discussion
Unlike studies addressing other extreme hydrometeorological hazards, research on mudflows in Peru remains scarce. This study identified areas exposed to mudflows through the integrated application of hydrological and hydraulic modeling in the CGALD. The following discussion examines the main impacts, methodological limitations, and key factors that should guide future research on mudflow hazards in arid and semi-arid environments of southern Peru.
4.1 Nonlinear hydrological response under extreme precipitation
The hydrological simulations reveal that the CGALD basin exhibits a markedly nonlinear response to extreme precipitation. While peak discharges increase progressively between RP25 and RP100, the RP200 scenario produces a disproportionate increase in discharge, rising from 25.3 m3/s (RP100) to 122.4 m3/s (RP200). This behavior suggests the existence of a hydrological threshold beyond which runoff generation and sediment mobilization become highly amplified.
The identified threshold response is consistent with the dominant hydrological processes in arid and semi-arid catchments (Montes-Pajuelo et al., 2024), where infiltration-excess (Hortonian) runoff governs storm response dynamics (Wang et al., 2017). Under such conditions, relatively small increases in rainfall intensity or duration may exceed infiltration capacity, rapidly increasing surface runoff connectivity and sediment entrainment. Similar nonlinear flood responses have been documented in steep arid basins where hydrological extremes are strongly controlled by short-duration convective storms and rapid concentration times (Tabari, 2020).
The precipitation record associated with the March 2020 event further supports this interpretation. The extreme rainfall value (~30 mm) identified in the frequency analysis lies within the upper tail of the distribution and corresponds to a documented mudflow occurrence in the study area. Its inclusion in the Log-Normal distribution fitting strongly influences the estimated magnitude of low-frequency events, emphasizing the sensitivity of flood estimation to rare but physically plausible extremes in data-scarce arid environments.
These findings indicate that conventional hazard assessments based exclusively on RP100 scenarios may underestimate the magnitude of exceptional events in southern Peruvian arid basins. The abrupt increase observed at RP200 highlights the need to incorporate extreme precipitation scenarios into hydraulic design, land-use planning, and risk management frameworks, particularly in rapidly urbanizing alluvial fan environments (Gioia et al., 2021).
4.2 Hydraulic behavior and geomorphological controls on mudflow propagation
The hydraulic simulations indicate that channel morphology exerts a primary control on mudflow propagation within the CGALD. For the RP25–RP100 scenarios, the main channel remains largely confined due to its incised geometry, the presence of lateral levees near the Arunta Bridge, and the relatively high conveyance capacity of the downstream reach (Pino, 2013).
A key geomorphological control identified in this study corresponds to the quarry sector located downstream of the bridge. The excavated depressions generated by extractive activities function as temporary storage zones that attenuate lateral flow expansion under moderate flood scenarios. However, their buffering capacity decreases substantially during extreme events characterized by high sediment concentration and elevated flow volumes.
Historical observations associated with flows exceeding approximately 200 m3/s indicate that overtopping and lateral spreading can occur in urban sectors adjacent to the active channel (Gómez Alanoca et al., 2025). This behavior is consistent with previous studies conducted in southern Peru, where unconsolidated alluvial materials, steep slopes, and sudden runoff inputs promote rapid increases in flow mobility during exceptional storms (Romero et al., 2013).
The simulations also demonstrate the importance of high-resolution topography for reproducing localized hydraulic controls. The RPAS-derived DEM enabled the identification of preferential flow paths, artificial depressions, and confinement zones that would likely remain unresolved using conventional topographic datasets. This represents a particularly relevant contribution for hazard assessment in data-scarce arid basins, where local microtopography strongly influences flow routing and inundation extent.
Model validation against the March 2020 event yielded an F1-score of 0.92, indicating excellent agreement between simulated and observed inundation extents. Although this validation does not constitute a full calibration, it supports the physical plausibility of the adopted hydrological and rheological parameterization.
4.3 Methodological implications and limitations
This study provides one of the first coupled hydrological–hydraulic assessments of mudflow hazard in southern Peru under data-scarce arid conditions. Nevertheless, several limitations associated with data availability and model parameterization should be acknowledged. The main source of uncertainty arises from the hydrometeorological input data. In arid environments, precipitation exhibits strong spatial and temporal variability that is difficult to capture using sparse rain-gauge networks and discontinuous historical records. As a result, the characterization of short-duration extreme storms remains uncertain, particularly for localized convective events that control rapid runoff generation and mudflow initiation.
Additional uncertainty is associated with the hydraulic parameterization. Due to the absence of site-specific rheometric and multi-temporal soil data, the non-Newtonian simulations relied on rheological parameters derived from the literature and calibrated within physically plausible ranges. Although this approach is acceptable for first-order hazard assessment in data-scarce regions, it may not fully represent the spatial variability of sediment concentration, yield stress, and flow mobility within the basin. Consequently, the proposed hydrological thresholds and inundation patterns should be interpreted as model-conditioned estimates rather than deterministic predictions.
Despite these limitations, the coupled HEC-HMS/HEC-RAS framework produced physically coherent results and satisfactory event-based validation against the March 2020 mudflow. This demonstrates the potential of integrated hydrological–hydraulic approaches for hazard assessment in arid Andean basins where direct observations are limited. From a methodological perspective, future studies should evaluate the sensitivity of the hydrological response using alternative semi-distributed platforms such as RS MINERVE (Astorayme, 2017), while precipitation gap-filling procedures should be contrasted with gridded climate products and alternative interpolation approaches to reduce uncertainty propagation in extreme-event analyses (Qquenta et al., 2023).
Future improvements should also incorporate field-based geotechnical characterization, rheometric testing of alluvial sediments, and high-resolution LiDAR topography to better constrain non-Newtonian parameters and improve the representation of flow propagation. Likewise, the implementation of specialized debris-flow models such as FLO-2D or RAMMS could provide more detailed simulations of velocity, runout distance, and depositional behavior. Strengthening meteorological monitoring networks, particularly at middle and upper elevations of the basin, together with systematic post-event documentation, would further improve future calibration and validation efforts.
Finally, future research should integrate complementary approaches focused on operational hazard assessment, including rainfall intensity–duration thresholds for mudflow initiation (Goyburo et al., 2025), probabilistic scenario analysis, and scenario libraries for rapid emergency response. These developments would strengthen the predictive capability and transferability of mudflow hazard assessments in southern Peru and other arid mountainous regions with limited observational data.
5 Conclusion
This article presents one of the first coupled hydrological–hydraulic assessments of mudflow hazard in a data-scarce arid basin of southern Peru. Integrating HEC-HMS and HEC-RAS with a non-Newtonian rheological approach, it provides new evidence on the hydrological controls and hazard dynamics governing extreme mudflow events in the Coronel Gregorio Albarracín Lanchipa District (CGALD). The principal conclusions are:
The basin exhibits a pronounced nonlinear hydrological response under extreme precipitation conditions. Peak discharges increase progressively from 11.9 m3/s (RP25) to 25.3 m3/s (RP100), whereas the RP200 scenario generates a disproportionate increase to 122.4 m3/s. This threshold-like behavior indicates that mudflow hazard in small arid basins cannot be adequately characterized through conventional linear extrapolation of design events.
The identified hydrological threshold suggests that assessments based exclusively on RP100 scenarios may substantially underestimate extreme-event magnitude and associated impacts. Relatively small increases in extreme precipitation can produce disproportionately large increases in runoff generation and flow discharge once critical conditions are exceeded.
Channel morphology exerts a dominant control on flow propagation. Although the incised channel and existing containment structures largely confine flows under RP25–RP100 scenarios, quarry-related depressions act as preferential accumulation zones that locally increase hazard levels.
Under the extreme historical scenario (~200 m3/s), potential impacts include approximately 300 dwellings, two populated sectors, one educational institution, 1.5 km of road infrastructure, and 44.6 ha of agricultural land. These findings highlight the increasing exposure associated with urban expansion over alluvial and geomorphologically active areas.
The coupled HEC-HMS/HEC-RAS framework, incorporating non-Newtonian Bingham rheology, achieved satisfactory agreement with the March 2020 event (F1-score = 0.92). The methodology therefore constitutes a robust first-order approach for mudflow hazard assessment in poorly monitored arid basins.
Overall, the results indicate that return periods ≥100 years represent a critical hazard threshold in the CGALD and potentially in similar arid basins of southern Peru. Consequently, hazard assessment, land-use planning, and risk-management strategies should explicitly incorporate extreme scenarios beyond conventional RP100 standards to avoid substantial underestimation of future mudflow impacts.
Statements
Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.
Author contributions
JQ: Conceptualization, Formal analysis, Investigation, Methodology, Software, Writing – original draft, Writing – review & editing. AG: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. EP-V: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing. VP: Conceptualization, Formal analysis, Writing – review & editing. JE-M: Formal analysis, Resources, Writing – review & editing, Funding acquisition, Project administration. KA-C: Formal analysis, Project administration, Writing – review & editing. GH: Formal analysis, Writing – review & editing. ET-A: Formal analysis, Writing – review & editing. WL-C: Conceptualization, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This work was funded by the Jorge Basadre Grohmann National University, Tacna, Peru, Vice-Rectorate of Research and Research Institute for providing the canon, canon and mining royalties funds for the development of the project “Mitigation of the risk of overflow and flooding based on a proposal of assessment of sustainable public space in arid areas”, approved with R. R. N°13626-2024-UNJBG.
Acknowledgments
The authors thank Jorge Basadre Grohmann University for its institutional and technical support, Senamhi, and the Water Research Group (H2O’UNJBG) for their collaboration, technical discussions, and support throughout this research.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that Generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
Astorayme (2017). Análisis y evaluación comparativa de modelos hidrológicos agrupados y semidistribuidos aplicados al pronóstico de caudales diarios del río Chillón. Lima: Universidad Nacional Mayor de San Marcos.
2
Aybar CamachoC. L.Lavado-CasimiroW.HuertaA.Fernández PalominoC.Vega-JácomeF.Sabino RojasE.et al. (2017). Uso del producto grillado PISCO de precipitación en estudios, investigaciones y sistemas operacionales de monitoreo y pronóstico hidrometeorológico. Nota Técnica No 001 SENAMHI-DHI-2017. Lima: Repositorio Institucional - SENAMHI.
3
AybarC.FernándezC.HuertaA.LavadoW.VegaF.Felipe-ObandoO. (2020). Construction of a high-resolution gridded rainfall dataset for Peru from 1981 to the present day. Hydrol. Sci. J.65, 770–785. doi: 10.1080/02626667.2019.1649411
4
BatesP. D.HorrittM. S.FewtrellT. J. (2010). A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling. J. Hydrol.387, 33–45. doi: 10.1016/J.JHYDROL.2010.03.027,
5
BilalM.KumarR. (2025). Physical modelling to map the rheological and morphological dynamics along with entrainment mechanics of debris flow: a comprehensive review of the state of the art. Geoenviron. Disasters 202512, 45–51. doi: 10.1186/s40677-025-00349-1
6
BrunoL. S.MattosT. S.OliveiraP. T. S.AlmagroA.RodriguesD. B. B. (2022). Hydrological and hydraulic modeling applied to flash flood events in a small urban stream. Hydrology9:223. doi: 10.3390/hydrology9120223
7
Cardich MottaK. A. (2017). Modelación De Máximas Avenidas En La Cuenca Del Río Lurín Utilizando Modelos Hidrológico E Hidráulico, 1–186. Available online at: https://hdl.handle.net/20.500.12996/3732
8
Del SavioA. A.AstocahuanaS. I. Q.NavarroL. F. C. (2019). Numerical simulation of debris flows of the catastrophic event of February 2019 in Mirave - Peru. Ambiente Agua14:1. doi: 10.4136/AMBI-AGUA.2437
9
El-BagouryH. A.G. (2024). Flood prediction, and mitigation using meteorological. doi: 10.3390/w16020356
10
FloresH.MedinaK.Castillo-VergaraF.IribarrenP.AzócarG.SalazarC.et al. (2026). Numerical simulation of rainfall-induced debris flows triggered by cyclone Yaku 2023 in Chasquitambo, Peru. Hydrology13, 83–21. doi: 10.3390/hydrology13030083
11
GaprindashviliM.TsereteliE.GaprindashviliG.KurtsikidzeO. (2021). NLandslide and mudflow hazard assessment in Georgia. Dordrecht: ATO Science for Peace and Security Series C: Environmental Security, 265–279. doi: 10.1007/978-94-024-2046-3_14
12
García (2016). RS MINERVE – Technical Manual. Technical Report. Sion: RS MINERVE Group.
13
GebremicaelT. G.DeitchM. J.GancelH. N.CroteauA. C.HaileG. G.BeyeneA. N.et al. (2022). Satellite-based rainfall estimates evaluation using a parsimonious hydrological model in the complex climate and topography of the Nile River catchments. Atmos. Res.266:105939. doi: 10.1016/J.ATMOSRES.2021.105939
14
GioiaA.LioiB.TotaroV.MolfettaM. G.ApollonioC.BisantinoT.et al. (2021). Estimation of peak discharges under different rainfall depth–duration–frequency formulations. Hydrology8:150. doi: 10.3390/HYDROLOGY8040150
15
Gómez AlanocaC.Quispe CáceresT.Carpio ObertiD.Rejas CéspedesK.Huanca RamosS. (2025). Determinación de zonas vulnerables a inundación con HEC-RAS en el tramo poblado del río Caplina, Tacna. doi: 10.5281/zenodo.16517291
16
GoodrichD. C.LaneL. J.ShillitoR. M.MillerS. N.SyedK. H.WoolhiserD. A. (1997). Linearity of basin response as a function of scale in a semiarid watershed. Water Resour. Res.33, 2951–2965. doi: 10.1029/97WR01422
17
GoyburoA.GutierrezL.RauP.Lavado-CasimiroW. (2025). Empirical rainfall thresholds for mudflow events in an arid basin of the Peruvian coast. Front. Water7:1637115. doi: 10.3389/frwa.2025.1637115
18
HawkerL.UheP.PauloL.SosaJ.SavageJ.SampsonC.et al. (2022). A 30 m global map of elevation with forests and buildings removed. Environ. Res. Lett.17:024016. doi: 10.1088/1748-9326/AC4D4F
19
HdeibR.AbdallahC.ColinF.BroccaL.MoussaR. (2018). Constraining coupled hydrological-hydraulic flood model by past storm events and post-event measurements in data-sparse regions. J. Hydrol.565, 160–176. doi: 10.1016/J.JHYDROL.2018.08.008
20
HighlandL. M. (2008). The Landslide Handbook—A Guide to Understanding Landslides. Available online at: https://pubs.usgs.gov/circ/1325/
21
HuangM. L.SunD. A.WangC. H.KeletaY. (2020). Reliability analysis of unsaturated soil slope stability using spatial random field-based Bayesian method. Landslides18, 1177–1189. doi: 10.1007/S10346-020-01525-0
22
HungrO.LeroueilS.PicarelliL. (2013). The Varnes classification of landslide types, an update. Landslides11, 167–194. doi: 10.1007/S10346-013-0436-Y
23
INEI. (2024). Tacna Compendio estadístico 2024. Available online at: https://www.gob.pe/institucion/inei/informes-publicaciones/6502390-compendio-estadistico-tacna-2024
24
INGEMMET (2021). Peligro geológico en la región Tacna, 1–238. Available online at: https://hdl.handle.net/20.500.12544/3161
25
IversonR. M. (1997). The physics of debris flows. Rev. Geophys.35, 245–296. doi: 10.1029/97RG00426
26
KeatonJ. R. (2019). Review of Contemporary Terminology for Damaging Surficial Processes: stream flow, Hyperconcentrated Sediment flow, Debris flow, mud flow, mud Flood, Mudslide. doi: 10.25676/11124/173147
27
KhamidullaevS.OymatovR.JulievM. (2025). Global trend and conceptual structure of mudflow hazard research: a bibliometric analysis using Biblioshiny. J. Geol. Geogr. Geoecol.34, 783–795. doi: 10.15421/112566
28
LeeD. H.LeeS. R.ParkJ. Y. (2025). Coupled modeling framework for proactive design of debris-flow barrier placements. Sci. Rep.15:290. doi: 10.1038/s41598-025-15290-4,
29
LimaD. C. (2025). Hydraulic modeling of Newtonian and non-Newtonian debris flows in alluvial fans: a case study in the Peruvian Andes. Water17:2150. doi: 10.3390/W17142150
30
LindsayJ. B. (2016). Whitebox GAT: a case study in geomorphometric analysis. Comput. Geosci.95, 75–84. doi: 10.1016/J.CAGEO.2016.07.003
31
MeselheE. A.HabibE.OcheO. C.GautamS. (2004). “Performance evaluation of physically based distributed hydrologic models and lumped hydrologic models,” in Proceedings of the 2004 World Water and Environmetal Resources Congress: Critical Transitions in Water and Environmetal Resources Management, 1–10. doi: 10.1061/40737(2004)211
32
Millán-ArancibiaC.Lavado-CasimiroW. (2023). Rainfall thresholds estimation for shallow landslides in Peru from gridded daily data. Nat. Hazards Earth Syst. Sci.23, 1191–1206. doi: 10.5194/NHESS-23-1191-2023
33
Montes-PajueloR.Rodríguez-PérezÁ. M.LópezR.RodríguezC. A.Montes-PajueloR.Rodríguez-PérezÁ. M.et al. (2024). Analysis of probability distributions for modelling extreme rainfall events and detecting climate change: insights from mathematical and statistical methods. Mathematics12:1093. doi: 10.3390/MATH12071093
34
MoreirasS. M.SepúlvedaS. A.Correas-GonzálezM.LauroC.VergaraI.JeanneretP.et al. (2021). Debris flows occurrence in the semiarid Central Andes under climate change scenario. Geosciences11, 43–27. doi: 10.3390/GEOSCIENCES11020043
35
NealJ.SchumannG.BatesP. (2012). A subgrid channel model for simulating river hydraulics and floodplain inundation over large and data sparse areas. Water Resour. Res.48:2514. doi: 10.1029/2012WR012514
36
O’BrienJ. S.JulienP. Y. (1988). Laboratory analysis of mudflow properties. J. Hydraul. Eng.114, 877–887. doi: 10.1061/(ASCE)0733-9429(1988)114:8(877)
37
O’BrienJ. S.JulienP. Y.FullertonW. T. (1993). Two-dimensional water flood and mudflow simulation. J. Hydraul. Eng.119, 244–261. doi: 10.1061/(ASCE)0733-9429(1993)119:2(244)
38
OMM. (2025). Guía de prácticas hidrológicas. Available online at: https://library.wmo.int/idurl/4/32737
39
OrellanaR. E. (2021). Modelamiento hidrológico e hidráulico para el análisis de inundaciones en la ciudad de Piura utilizando HEC-HMS y HEC-RAS. Lima: Pontifia Universidad Católica Del Perú.
40
OrtegaJ. C. B.BendezuM. A. L.Del SavioA. A.CanalesF. A. (2022). Effect of lithological and geotechnical characteristics on the generation of debris flows in the arid basin of Mirave, Peru. Ambiente Agua17, 1–18. doi: 10.4136/AMBI-AGUA.2785
41
Pekerİ. B.GülbazS.DemirV.OrhanO.BedenN. (2024). Integration of HEC-RAS and HEC-HMS with GIS in flood modeling and flood hazard mapping. Sustainability16:1226. doi: 10.3390/su16031226
42
PinoT. C. A. (2013). Caracterización hidrogeomorfológica de la cuenca del Río Caplina – Tacna. Tacna: Universidad Nacional Jorge Basadre Grohmann.
43
QquentaJ.RauP.BourrelL.FrappartF.Lavado-CasimiroW.QquentaJ.et al. (2023). Assessment of bottom-up satellite precipitation products on river streamflow estimations in the Peruvian Pacific drainage. Remote Sens.16:11. doi: 10.3390/RS16010011
44
ReidM. E.CoeJ. A.BrienD. L. (2016). Forecasting inundation from debris flows that grow volumetrically during travel, with application to the Oregon coast range, USA. Geomorphology273, 396–411. doi: 10.1016/J.GEOMORPH.2016.07.039
45
RomeroH.SmithP.MendonçaM.MéndezM. (2013). Macro y mesoclimas del altiplano andino y desierto de Atacama: desafíos y estrategias de adaptación social ante su variabilidad. Rev. Geogr. Norte Gd.55, 19–41. doi: 10.4067/S0718-34022013000200003
46
RoqueG. (2022). Riesgo de inundaciones fluviales por máximas avenidas en la cuenca baja de Rio Lurin. Available online at: https://hdl.handle.net/20.500.13084/6031
47
Salazar-BrionesC.Hallack-AlegríaM.Mungaray-MoctezumaA.LomelíM. A.Lopez-LambrañoA.Salcedo-PerediaA. (2018). Modelación hidrológica e hidráulica de un río intraurbano en una cuenca transfronteriza con el apoyo del análisis regional de frecuencias. Tecnol. Cienc. Agua9, 48–74. doi: 10.24850/j-tyca-2018-04-03
48
TabariH. (2020). Climate change impact on flood and extreme precipitation increases with water availability. Sci. Rep.10:13768. doi: 10.1038/s41598-020-70816-2,
49
TsereteliE.BolashviliN.GaprindashviliG.GaprindashviliM. (2022a). Mudflow processes in Georgia. Geography Water Res.2, 28–34. doi: 10.55764/2957-9856/2022-2-28-34.10
50
TsereteliE.GaprindashviliG.GaprindashviliM. (2022b). Natural Disasters (Mudflow, Landslide, Etc.), 55–69. doi: 10.1007/978-3-030-90753-2_7
51
USACE-HEC. (2025a). HEC-HMS User’s Manual. Available online at: https://www.hec.usace.army.mil/confluence/hmsdocs/hmsum/latest
52
USACE-HEC. (2025b). HEC-RAS User’s Manual. Available online at: https://www.hec.usace.army.mil/confluence/rasdocs/rasum/latest
53
VallejoL. E. (1979). An explanation for mudflows. Géotechnique29, 351–354. doi: 10.1680/GEOT.1979.29.3.351
54
Vilcanqui-AlarcónA.Pino-VargasE.Vargas-BernuyJ. (2022). Geomorphological alteration in relation to anthropic actions in the Caplina riverbed, Tacna, Peru. Agroindustrial Sci.12, 47–58. doi: 10.17268/agroind.sci.2022.01.06
55
WangW.LiH. Y.LeungL. R.YigzawW.ZhaoJ.LuH.et al. (2017). Nonlinear filtering effects of reservoirs on flood frequency curves at the regional scale. Water Resour. Res.53, 8277–8292. doi: 10.1002/2017WR020871
56
WeiZ. L.ShangY. Q.LiangQ. H.XiaX. L. (2024). A coupled hydrological and hydrodynamic modeling approach for estimating rainfall thresholds of debris-flow occurrence. Nat. Hazards Earth Syst. Sci.24, 3357–3379. doi: 10.5194/nhess-24-3357-2024
Summary
Keywords
arid basins, couple modeling, hazard mapping, mudflow, nonlinear hydrological response
Citation
Qquenta J, Goyburo A, Pino-Vargas E, Pocco V, Espinoza-Molina J, Acosta-Caipa K, Huayna G, Taya-Acosta E and Lavado-Casimiro W (2026) Mudflow hazard mapping and nonlinear hydrological response in arid basins of southern Peru using coupled hydrological–hydraulic modeling. Front. Water 8:1822397. doi: 10.3389/frwa.2026.1822397
Received
03 March 2026
Revised
02 June 2026
Accepted
10 June 2026
Published
01 July 2026
Volume
8 - 2026
Edited by
Nanditha J. S., Indian Institute of Technology Kanpur, India
Reviewed by
Bidhan Kumar Sahu, Indian Institute of Technology Gandhinagar, India
Elmer Calizaya, Universidad Nacional del Altiplano, Peru
Updates
Copyright
© 2026 Qquenta, Goyburo, Pino-Vargas, Pocco, Espinoza-Molina, Acosta-Caipa, Huayna, Taya-Acosta and Lavado-Casimiro.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jonathan Qquenta, jonathan.qquenta@unmsm.edu.pe
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.