Arequipa, Peru’s second-largest city (1,081,000 inhabitants as per the 2017 national census; INEI, 2018), is one of the fastest-growing regions in the country, at the contact between the Andean highlands and the coast. Arequipa is situated at c. 2,300 masl in a pull-apart basin between the west flank of the Western Cordillera of the Central Andes and the Coastal batholith to the south. This basin was filled by volcanic deposits, mostly Neogene and Pleistocene ignimbrites, covered in part by volcanoclastic deposits that make up the foothills of the Nevado Chachani, El Misti, and Pichu volcanoes. Chachani (6,080 masl) and Pichu (5,600 masl) are extinct volcanic clusters, but El Misti active composite stratocone currently exhibits seismic and fumarolic activity (Aguilar Contreras et al., 2022). Arequipa is the largest South American city exposed to a large variety of natural hazards from El Misti volcano, earthquakes, flash floods, and mass flows. Between 50,000 and 100,000 inhabitants living <1 km from active river or ravine channels are likely to be affected by flash floods and/or lahars, and many more if El Misti volcano, located 17 km NE of the city of Arequipa, produced tephra fall and PDCs. Socio-economic hardships are added to environmental hazards because poor urban development and planning have led to informal settlements and marginal contributions to improving living conditions (INEI, 2018).

The multidisciplinary project “Impacts in the City of Arequipa” aims to address hazards, exposure, and vulnerability of buildings and infrastructure that face frequent flash floods, hyperconcentrated flows (HCFs hereafter), and debris flows (DFs hereafter) along five ravines or quebradas that cross the city. The five ravines are tributaries to Río Chili (inset map, Figure 1), whose discharge is artificially controlled for electric supply and water intake. Impacts of potential tephra fallout have also been addressed here so that probabilities of roof collapse in case of eruption may enable stakeholders to undertake retrofitting projects and improve risk management in urban planning.

The research objectives are threefold: 1) to assess the extent of impacts of potential tephra fall on building roofs based on earlier classification of the building stock; 2) to characterize and simulate mass flows, that is, HCFs and DFs in two ravines often affected by floods in the recent past: Quebrada (Qda.) San Lazaro and Qda. Huarangal-Mariano Melgar (Huarangal-MM thereafter), draining the SSW and South slopes of El Misti, respectively (Figure 1); and 3) to assess impacts of these flows on the building stock in order to contribute to potential loss assessment.

Flash floods are local floods triggered within 6 h and often within 3 h, generally resulting from heavy rainfall over a catchment (National Weather Service, 2022). Flash floods frequently occur in arid or semi-arid environments where rainfall is concentrated over a short period of time, and the infiltration capacity of the soil is limited so that runoff is triggered shortly after the rainfall onset on relatively steep slopes. Flash floods can occur in cities where construction and pavement made the ground impervious. The environment of the city of Arequipa is propitious to flash floods and mass flows as the result of the semi-arid climate (c. 150 mm a year), the concentration of rainfall as thunderstorms often delivering more rainfall than the average monthly amount. A series of factors favor FFs, HCFs, and DFs across Arequipa (See Table 1): (i) short-lasting (3–4 h) and intense rainstorms (26–42 mm/h, from five meteorological stations across the basin of Arequipa over the past 60 years; Fuse and Benites, 2003) typical of a semiarid climate (≤150 mm/year at 2,300 masl); (ii) rapid runoff (30–60 min) in relatively small mountainous catchments (18–40 km²); (iii) the impervious network of roads and waterproofing of the built-up area that increase runoff; and (iv) the steep, bare slopes of El Misti’s foothills carved by ravine channels, which are filled by loose volcanic debris mixed with household and industrial waste (Thouret et al., 2013, 2014; Ettinger et al., 2015; Thouret et al., 2018; Arapa et al., 2019).

Lahars (volcanic debris flows, LHs) and rain-triggered debris flows (DFs) are potentially destructive water-saturated mass flows acting in mountainous regions and areas where steep slopes cut down in loose sediment. These flows can be divided into two categories using sediment concentration, grain-size distribution, and bulk density (Vallance, 2000; Hungr and Jakob, 2005; Manville et al., 2013; Vallance and Iverson, 2015; Thouret et al., 2020):

1. Debris flows (DFs) are mixtures of debris and water with high sediment concentrations that move downslope due to gravity as surging sediment slurries (Vallance, 2000). DFs comprise a solid phase of at least 60 vol% (>80 wt%), thoroughly mixed with water. The solid component includes mostly gravel and boulders with sand, silt, and clay proportions remaining low. A threshold of 3 wt% of silt and clay helps distinguish non-cohesive from cohesive DFs (Scott, 1988; Scott et al., 2005). The density of a DF ranges between 1,800 and 2,400 kg/m³, twice as much as the muddy water that flows in stream channels during floods (Pierson, 1980).

2. Hyperconcentrated flows (HCFs) are two-phase flows intermediate in sediment concentration between “normal” streamflows and DFs, with densities between 1,300 and 1,800 kg/m³. HCFs transport between 20 and 60 vol% (40 and 80 wt%) of sediment (Beverage and Culbertson, 1964; Pierson, 2005).

1. We used probabilistic maps of tephra fall extent and load from El Misti on the city and surrounding areas produced by Sandri et al. (2014) . Taking into account past events according to the recorded eruptive styles, the authors ran many simulations using TEPHRA2 (Connor and Connor, 2011) for each type of eruption and volumes of tephra fall between 10⁶ and 10⁹ m³. Using a Bayesian event tree and elicitation method, the authors measured a 0.05% tephra-fall probability of occurrence over 1 year for Arequipa in case of Vulcanian, SubPlinian, or Plinian eruption, particularly during the wet season (December–March). Their study provided four maps showing the areal distribution of probabilities of tephra fall with two mass loading thresholds, 100 and 250 kg/m², of tephra loads exceedance on roofs for both rainy and dry seasons and two, 10 km and 20 km columna heights (Figures 6, 7 in Sandri et al., 2014).

2. We assessed the vulnerability of roofs with respect to tephra load, including two steps: i) field surveys helped us recognize nine types of constructions and roofs based on several criteria collected by means of field surveys across 300 city blocks and ii) classification of about 2,000 buildings was achieved by means of field surveys at the building scale, completed by remote sensing using high-spatial-resolution Pléïades satellite images and drone images (50 and 25 cm pixel, respectively) and Google Street View (Figures 3A–D).

1. We elaborated a database of flow events over the period 1915 to 2020 to highlight hydrological input parameters that will be used in numerical simulations (See Table 1; Supplementary Table S1).

2. We chose two catchments, namely, the Qda. San Lazaro and the Qda. Huarangal-MM, which initiate >3,300 masl on the SW flank of El Misti volcano, for two reasons: they have been often affected by flash floods and mass flows, and their channels run across a wide range of buildings representing different periods of construction from suburbs to the historic center of Arequipa in the medial and distal reaches, respectively. All ravines, usually dry, become active during the rainy season (December to March) when they convey flash floods and mass flows to Río Chili. Table 2 displays a few morphometric characteristics of both catchments.

Qda. San Lazaro forms high (3,680 masl) on the SW flank of El Misti volcano and follows a NE-SW direction over a length of 14.1 km with an average slope of c. 8° from its source to its confluence with Río Chili near the Grau Bridge (2,372 masl) (Figure 2A). The catchment area is 18.80 km². The ravine has incised past lahar and pyroclastic-flow deposits a few tens of meters thick, forming two steep (9%) and widespread (5.50 km × 1.80 km) fans, which have been entirely built over the past 50 years. Qda. San Lazaro has been diverted to the SW, 1 km upstream of the confluence with Río Chili, by the older, eastern fan, which merges down valley with the high terrace supporting the historic center (Thouret et al., 2014). The longitudinal profile of San Lazaro shows three reaches and a gradient of 0.10 (Figure 2A). Nine bridges were built over Qda. San Lazaro, the most vulnerable being José Olaya, Parque Selva Alegre, and Ingreso Universidad Católica San Pablo (Figure 2A; Thouret et al., 2014).

The source of Qda. Huarangal-MM is located at 3,270 masl on the South flank of El Misti volcano, and the ravine follows an ENE–WSW direction, over a length of 18.9 km with an average slope of 5.3° from its source to the confluence with Río Chili near the suburb Tingo (2,259 masl; Figure 2B). Several tributaries join Qda. Huarangal-MM, but the c. 50-m-wide channel narrows down valley near the Geronimo–Paucarpata suburb. As a result, the Huarangal-MM catchment is more extensive than that of San Lazaro (Figure 2C; Table 1), while its longitudinal profile with a gradient of 0.05 is half as steep as its counterpart (Figure 2B). Five amongst ten bridges that cross the ravine are more vulnerable: 8 de Octubre, Tupac Amaru, Jorge Chavez, Santa Rosa, and Gran Unidad Escolar (Figure 2B; Thouret et al., 2014).

3. We utilized three numerical codes, Titan2F, VolcFlow, and FLO-2D, to simulate mass flows. Our purpose is to test the abilities of these numerical codes to simulate DF and HCF events and detect their limitations compared with real cases derived from field data or previous simulations described in the literature.

4. We assessed mass-flow impacts on habitat and drew fragility curves involving hydrodynamic and hydrostatic pressure (via flow depth) from DF and HCF, resulting from simulations and earlier studies (Mead et al., 2017; see Thouret et al., 2020 for a review of lahar characteristics and impacts). This assessment was completed by in situ measurements on home and bridge construction material, both on walls and common construction material, using devices such as sclerometers for hardness and compressive strength, an ultrasonic analyzer for crash resistance, Young’s modulus, and surface velocity (Chehade, 2021).

5. Using results from fragility curves for two categories of building wall thickness, we computed loss fractions for the most vulnerable classes A0, A, and B according to the position of buildings adjacent to two ravine channels and the delineated area of simulated dynamic pressure near both ravine channels. From the building loss fractions, we estimated the cost of reconstruction or retrofitting for each vulnerable class of A0, A, and B according to the percentage of building loss fraction in case of DF and HCF.

Potential volcanic hazards from the active El Misti are acute, as the recorded eruptive activity exhibited a variety of styles and deposits (Thouret et al., 2001; Harpel et al., 2011; Charbonnier et al., 2020). Among the volcanic hazards (Supplementary Figure S1), we focus on tephra fallout, while modeling of PDCs has been the target of recent work (Charbonnier et al., 2020). Tephra fall is considered the most common hazard from volcanic activity. It has rarely caused human life loss, but the high mobility of small particles is responsible for health problems and has economic impacts in cities around active volcanoes. Globally, the impacts of tephra fall include building roof collapse, damage to aircraft, disruption of lifelines, and health and environmental effects. Tephra fall from El Misti occurred in the mid-15th century in Arequipa and fine ash from Nevado Sabancaya in 2016. However, explosive activity over the past 50,000 years has shown a succession of Vulcanian to Plinian eruptions and a frequency of about 500 years for VEI 2 events and 2,000–5,000 years for VEI 3-4 events (Thouret et al., 2001; Sandri et al., 2014).

Holocene and historical Misti’s activity has encompassed a range of eruptions with VEI 2-4 magnitudes. Two recent eruptions (i.e., the 1440–1470 CE Vulcanian event and the ca. 2070-year BP Plinian eruption) are working references for probable scenarios (e.g., Sandri et al., 2014). Tephra-fall deposits at least 10 cm thick were emplaced by the mid-15th-century event (Chávez, 1992; Thouret et al., 2001), while 2070-year BP pumice-fall deposits exceeding 50 cm in thickness have been observed across the city area (Cobeñas et al., 2012; Cobeñas et al., 2014; Harpel et al., 2011; Harpel et al., 2013). Dispersal of tephra fallout toward the city depends principally on the eruptive style, column height, and patterns of wind directions and velocities over southern Peru.

Mass flows including DFs and HCFs threaten the city area, to which flash floods (FF) are added in the semi-arid environment of southern Peru. We selected 39 FFs and mass-flow events over the period 1915–2020 from reports, newspapers, and previous studies (e.g., Fernandez-Davila and Benites Montufar, 2003) that we completed with INDECI (National Institute of Civil Defense) reports (2019) on disasters (Table 1; Supplementary Table S1). As a result, FFs and lahars (HCFs or DFs) occur every 2.7 years on average and can be lethal and destructive every 5 years on average (Table 1; Nagata, 1999; Thouret et al., 2013; Thouret et al., 2014; Xue, 2016). The disrupting flows affected about 33,000 inhabitants and damaged more than 7,600 buildings over the past 105 years (Supplementary Table S1; Regional de Arequipa, 2014; INDECI, 2019). One of the destructive and lethal events occurred on 8 February 2013 when heavy rainfall (124.5 mm during 3 h) triggered an HCF and extensive flood along the Qda. Venezuela, affecting 280 buildings and 23 bridges in the city (Ettinger et al., 2015).

Structural classification of nine roof types (Table 4), based on Blong et al. (2017) classification and field surveys, encompasses mostly habitat buildings and administrative, educational and religious edifices, historic monuments, and sport facilities. Types 1, 2, 3, 4, and 6 are houses and residential buildings, while types 5, 7, 8, and 9 correspond to commerce and other edifices (Supplementary Table S4; Supplementary Figure S3). We defined five types of residential habitats ranging from houses with poor quality of construction material (e.g., type 4) to residential, 2 + story buildings for which design, structure, and material offer better resistance to load and impacts. Poor quality and vulnerable houses either have roofs made up of light material (plastic, wood, zinc, or a combination of those) or are often un-finished constructions without a roof. Type 2 average quality houses and commerce or industrial buildings have metal roofs and better-quality materials. Types 1 and 3 are high-quality and resistant buildings, while type 3 edifices show a better design and roof resistance due to pitch and reinforced concrete. Type 6 includes high-rise edifices with reinforced concrete, and potential impacts from tephra fallout have been estimated in terms of probabilities. The impacts would be more important during the rainy season when wind patterns would disperse more tephra toward the city area (Figures 6, 7inSandri et al., 2014). We assessed the extent of roof damage by comparing the tephra load (in particular 250 kg/m²) with the least resistant buildings and roofs that belong to vulnerability classes A0 and A. The latter classes prevail along the upper and middle reaches of both quebradas (Supplementary Figure S2).

We refer to Sandri et al. (2014) for the probability maps of tephra fall with 10 or 20 km high eruption column heights, wind patterns, and velocities during two dry and wet seasons in southern Peru and for two hazard metrics (tephra loads of 100 and 250 kg/m²). We evaluated the structural performance and susceptibility of roof collapse across built areas along both banks of Qda. San Lazaro and Qda. Huarangal-MM. Each building roof was classified based on seven criteria, including quality of the construction material, usage, and surface area. We obtained nine types of roofs mapped in the Q-GIS software.

The nine roof types have been converted into four vulnerability classes (Blong et al., 2017) according to the mechanical resistance of the material and based on two hazard metrics 100 and 250 kg/m². These tephra loads, in turn, yield pressures of 1–2 kPa and forces of 1,000 and 2,500 N/m² (inset graphics, Table 5). A substantial proportion (38% in area) of regular-quality roofs may not withstand 250 kg/m² or 2.5 kN/m², and almost all low-quality roofs (11%) may suffer structural damage or collapse if the tephra load reaches 100 kg/m² or 1 kN/m². Particular attention was paid to light and heterogeneous material covering low-quality houses, ranked as “hybrid” roofs, which can be affected by minor tephra loads ≤100 kg/m². Flat roofs or upper ceilings without proper cover may be more vulnerable than pitched roofs, a category that represents 3%-4% (in area) only of the surveyed buildings located in high-income neighborhoods only. Both values of hazard metrics may increase in case of rainfall concomitant with tephra fall because the probability of ashfall toward the city increases during the rainy season (Sandri et al., 2014). Both values of forces (1–2.5 kN/m²) are 2.5–6 times lower than forces exerted by tephra thickness between 20 and 50 cm in case of Plinian tephra fall (inset graphics, Table 5), to be compared to the 50-cm-thick 2070 BP pumice-fall deposit observed in the north and NW areas of the city.

The vulnerable roof classes rely on structural description and mechanical resistance (Table 4), which has led to impacts according to damage level in case of tephra fall exceeding the tephra loads and forces defined above (Figures 3C,D). From maps of vulnerable roofs and vulnerability index (Figures 3A,B; Table 5), we estimate three damage levels in case of the higher hazard metrics (load of 250 kg/m² defined above): 1) heavy/severe (i.e., roof collapse and burnt timber, failure of trusses and supporting structure, severe to moderate damage to rest of the edifice) for roof types 4, 7, and 9 (c. 11% in area); 2) significant damage to the roof structure and some damage to walls, at least one wall damaged/misaligned, collapse of part of ceiling for roof type 2 (c. 13%); and 3) slight (bending, damage to roofs overhangs and ceilings without roof, required cleaning, repainting or material overhaul) for roof types 1 and 5 (c. 14%). If the tephra load exceeded 250 kg/m² and reached about 6 kN/m², which would be the case with a 50-cm-thick tephra fall comparable to the c. 2070 BP Plinian pumice-fall (i.e., between 5 and 6 kPa), the aforementioned roof types (categories A, B, and C; Table 4) would all experience collapse, while moderate-to-severe damage would affect the rest of edifices. If we consider the probability of roof collapse triggered by a combination of tephra load (250 kg/m²) and rainfall, the proportion will increase to the majority of roof categories beyond classes A, B, and C, except for the RC (reinforced concrete) and pitched types of roofs, category D.

We estimated the cost of roof construction in case of damage due to tephra mass loading (inset graphics, Table 5), considering at least 425 building roofs along the medial reach of Qda. San Lazaro together with the proportion and estimated cost range of each roof type (Table 5). Severe damage of 75% of low-quality roof types 4, 7, and 9 would amount to a range of 18,000–20,000 US$. Significant damage of 50% of regular-quality roof type 2 would reach 6,000–8,000 US$, while slight damage of 25% of average-quality roof types 1 and 5 would translate into a cost range of 75,000–80,000 US$.

Flash floods were the most frequent events (55%), while HCF represented 33% and DF only 12% of the total events. Because solid concentration was not measured in situ, DFs with at least 40% vol. solid concentration may actually be more frequent. All flow events occur during the rainy, austral summer, with the majority (54%) in February, 33% in January, and 13% in March. Over 105 years, one event occurred every 2.70 years on average. A significant flow event (i.e., with >1 fatality and damaged buildings and/or bridges) occurred every 5 years on average.

Ettinger et al. (2015) reported some characteristics of the 8 February 2013 HCF. Total volume of the flow ranged between 50,000 and 100,000 m³. Peak flow discharge increased from 124 m³/s in the upstream reach to 425 m³/s in the middle reach, and the average discharge was seven times that of average discharge (19,50 m³/s) of streamflows (maximum 43 m³/s in San Lazaro and 69 m³/s in Huarangal-MM; Fuse and Benites, 2003) occurring in the rainy season. Velocity ranged between 5.9 and 10.9 m/s and flow depth between 3 and 6 m. The return period of discharges in Quebradas using the Gumbel distribution was compared to the case study of Río Chili (Martelli, 2011) and estimated on the basis of discharges reported from the 2013 event (Table 6).

In recent decades, various computational models have been developed that have allowed simulating the physics of mass flows from the resolution of equations by numerical methods (Iverson, 2014). Computational models or numerical codes have become one of the most important tools for HCF/DF hazard assessment. Thouret et al. (2020), their Table 4 elaborated on a list of models used for simulating lahars or HCFs/DFs, such as LAHARZ-py (Schilling, 2014), Titan2D (Pitman and Le, 2005; Williams et al., 2008) and Titan2F (Córdoba et al., 2015; Córdoba et al., 2018), FLO-2D (O’Brien et al., 1993), RAMMS (Cesca and D’Agostino, 2008), and VolcFlow (Kelfoun and Druitt, 2005). We tested three models, Titan2F, VolcFlow, and FLO-2D, all shallow-water depth-averaged models for modeling HCFs and DFs in the Qdas. San Lazaro and Huarangal. Assuming that HCF or DF thickness is much smaller than its length, it is possible to integrate the 3D mass and momentum balance equations over depth to obtain the depth-averaged continuum flow equations (Savage and Hutter, 1989). All simulations were performed over the 2013 DSM of the El Misti-Arequipa area (Oehler et al., 2014; Charbonnier et al., 2020) after smoothing and resampling the grid to 4 m spatial resolution and removal of all bridges crossing both Quebradas to avoid artificial blocking of the channeled flows. Charbonnier et al. (2020) used another couple of Pléaïdes images to enlarge the DSM to the SW flank and a TanDEM-X (12 m) to include the NE of the city, Huarangal-MM catchment, and El Misti’s summit. This has led to a DSM where bridges and other obstacles have been removed to simulate “free” flows, avoiding artificial damming and thickening of flows in case of artificial obstacles. We also ran simulations using the original DTM in which bridges and obstacles were present: the artificial obstacles (i.e., bridges digitally represented as vertical walls) dammed channeled flow deposits and induced artificial overbank flows with high velocities (5–7.5 m/s). Based on the outputs of the computational models, we defined mass-flow inundation zones, flow depths, and dynamic pressure and compared the results from the three codes despite their different purposes and input parameters.

The computational models require several input parameters (either from measurements in river channels and flumes, empirical models, or literature) and initial and boundary conditions, which are different for the three models (Table 7). In order to apply reliable rheological models to the modeled flows, we defined the density, viscosity, and yield strength parameters for each type of flow. Among required input parameters, flow rate and sedimentation concentration are critical for simulating sediment-water flows. Despite the lack of recorded flow discharges (except for streamflows and a few recent events, e.g., the 8 February 2013 HCF), we estimated discharges for both Quebradas based on the available discharge data of Río Chili and the Gumbel law distribution (Section 5.1; Table 6). The paucity of sedimentation concentration values makes it difficult to apply rheological properties and determine flow behavior, representing critical input information for modeling. Pallares et al. (2015) found high sand and gravel contents and a minimum amount of fine particles in historical lahar deposits across the city of Arequipa, although it is known that the grain-size distribution of lahar deposits does not reflect the wide flow spectrum that emplaced them (Dumaisnil et al., 2010). In case of heavy rainstorms (≥41 mm/h) over the small catchments under study, erosion and incorporation of coarse-grained volcaniclastic deposits to lahars may produce a limited increase in fine-grained particles (<8% silt and clay; Pallares et al., 2015). A small proportion of clay and silt in flows would influence the behavior of less cohesive HCFs and DFs (Bardou et al., 2007). Thus, we derived additional physical and rheological parameters of the flows to be simulated from previous studies of lahar characteristics at Sarno, Italy (Zanchetta et al., 2004); Mt. Ruapehu, New Zealand (Manville et al., 1998); Semeru, Indonesia (Dumaisnil et al., 2010; Doyle et al., 2010, 2011); and Panabaj, Guatemala (Charbonnier et al., 2018). HCFs prevailed among lahar flows and deposits measured around Semeru, which provided characteristics similar to recent flows observed in Arequipa, such as the 8 February 2013 HCF (Ettinger et al., 2015).

We mapped scenario-based inundation zones and evaluated impacts on the building stock using simulation runs of real and potential mass flows obtained from the three numerical codes. Three scenarios for HCFs and DFs have been derived from the recorded small-to-moderate events with 1 ̶ 1.5 × 10⁵ m³ (e.g., the 8 February 2013 HCF) and moderate-to-large historical lahars (c. 2.5 m³ ̶ 3.5 × 10⁵ m³ to 5 × 10⁵ m³).

Titan2F was developed at SUNY University at Buffalo for the modeling of biphasic (sediment-water) gravity-driven flows. The Titan2F model uses a similar framework to that developed by Pitman and Le (2005) for Titan2D. However, Titan2F has a complete set model of equations for both the granular phase and fluid phase (Córdoba et al., 2015; Córdoba et al., 2018). Titan2F is based on a depth-averaged, shallow layer model, as the fluid phase of the flow is modeled using a typical shallow water layer approach, while the frictional behavior of the solid phase is accounted for during the propagation of the simulated flow by assuming the Mohr–Coulomb model (Savage and Hutter, 1989).

VolcFlow was developed at the Laboratoire Magmas et Volcans to simulate the complex behavior of volcanic flows (Kelfoun and Druitt, 2005; Kelfoun et al., 2009; Kelfoun, 2017). It is based on the shallow water depth-averaged equations of mass and momentum balance, and the code allows the use of different rheological laws and source conditions. The equations are solved in a MATLAB^® environment using a numerical method based on an Eulerian scheme by finite difference. The input parameters are the volume rate with time, and the flow is defined as a viscoplastic Bingham fluid (Rickenmann, 1991; Phillips and Davies, 1991; Julien and Lan, 1991; Major and Pierson, 1992) with a yield strength and viscosity (Table 9), combined with a turbulent-dispersive law that considers the effects of turbulence and dispersive stresses caused by sediment collisions. In the depth-averaged form, it is expressed as follows: τ = τ 0 × u ‖ u ‖ + 3 × μ × u h + c × ρ × u × ‖ u ‖ , (1) where τ ₀ is the yield strength, μ is the viscosity, ρ is the flow density, h its thickness, and u = ( u x , u y ) its velocity. The coefficient c is a dimensionless number that can be linked to the coefficient ξ of Voellmy ( c = g ξ ) and that rules the turbulent dispersive pressure. The VolcFlow version used here is the single-phase version, in contrast to the two-phase regime used by Titan2F and FLO-2D codes.

Other input parameters used in VolcFlow are as follows (Table 9): 1) the same 4 m resolution DSM of the Arequipa area used for other codes and 2) an inflow hydrograph (inset, Figures 5A,B) imposed as a boundary condition and based on the measured rainfall amount and values of water and sediment discharge (yielding a flow volume). For this study, the volume rate Q of the hydrograph is approximated by the following law: Q = k 1 × t k 2 e − t / k 3 , (2) where k ₂ = 2.5, k ₃ = 200, and k ₁ is computed to obtain the total flow volumes ranging from 1.5 × 10⁵ m³ to 5 × 10⁵ m³. Given these parameters, the volume rate increases to reach a peak discharge at 500 m³/s and then decreases exponentially during ∼2,000 s. The physical and rheological values, shown in Table 9, have been selected from empirical and experimental values calculated by Rickenmann (1991), Phillips and Davies (1991), Pierson (2005), Iverson et al. (2010), Iverson et al. (2011), and Manville et al. (2013).

Modeling scenarios were established for a set of HCFs and DFs in both Qda. San Lazaro and Huarangal (Table 9). In total, three scenarios, each including two types of mass flows, were defined for both ravines with a similar range of volumes (1.5–5 × 10⁵ m³). We assume the worst-case scenario WF3 to be the result of maximum volumes (4.5 to 5 × 10⁵ m³) with a density of 2,200 kg/m³, yield strength between 20 and 200 Pa, and viscosity between 1 and 20 Pas (Table 9). This scenario represents HCF and DF events triggered by heavy, durable rainstorm remobilizing abundant, loose material emplaced by post-eruptive lahars.

The results of numerical simulations corresponding to HCF showed a realistic extent of the impacted areas (Figure 5A) and reasonable values of velocity and flow depth compared to the observed flows from the database. However, overbanks were observed at the large bend, followed by channel constrictions near the J. Olaya bridge in the middle reach of Qda. San Lazaro (Figure 2A) when the hydrograph reaches a maximum flow rate of 500 m³/s (inset Figures 5A,B).

The distribution of hydrodynamic pressures is between 35 and 100 kPa (Figure 5C) across the channel and low terraces 2–4 m high, where the most vulnerable houses (classes A0, A, and B) are located. For volumes in the order of 2.5–3.5 × 10⁵ m³, the large-volume scenario shows that major overbanks occurred from a sharp bend located just upstream of the José Olaya bridge and another just upstream El Progreso bridge (Figures 2B, 5B).

FLO-2D (FLO-2D Software Inc, 2021; O’Brien et al., 1993) is based on depth-averaged continuum flow equations. This code models lahar rheology using a shear stress relationship, and this rheology depends on the sediment concentration and flow mass, which vary during the flow propagation. Details about the code itself can be found in O’Brien et al. (1993), while the FLO-2D code was applied, for instance, by Charbonnier et al. (2018) to simulate the 2005 lahars at Panabaj (Guatemala). Two separate case studies are investigated here for the emplacement of 1) three DFs in the San Lazaro and Huarangal catchments and 2) three HCFs in the same catchments. The DFs and HCFs separately exhibited the same rheology and sediment concentration (i.e., all DFs were the same and all HCFs were the same). However, they each had three different volumes/discharges corresponding to small-, moderate-, and large-volume events.

FLO-2D includes various input parameters: 1) the DTM/DSM of the study area; 2) an inflow hydrograph based on measured rainfall accumulation and values of water and sediment discharge (i.e., flow volume); 3) a range of Manning roughness coefficients (0.03–0.08 according to the channel geometry and material); 4) yield strength and dynamic viscosities expressed as a function of the volumetric sediment concentration: τ y   =   α i e β i C v , (3) η   = α j e β j C v , (4) where τ _y is the yield strength, η is the dynamic viscosity, C _v is the volumetric sediment concentration, and α _i,j and β _i,j are two empirical coefficients obtained from O’Brien and Julien (1988) for HCFs and DFs with similar characteristics as the Misti ones; (5) the dimensionless laminar flow resistance K (Table 7); and (6) a Froude number, a dimensionless value defined as the ratio of kinetic energy to the potential energy that also accounts for the effect of gravity (Table 7 for best-fit input values). Inflow hydrographs for each simulation are shown in Figures 6A,B. Peak discharge rates (water only) vary from 45–60 m³/s (small-volume events) to 160–210 m³/s (large-volume events) for both DF and HCF simulations, while peak sediment volumetric concentrations (attained after ∼25 min) are set at 60 and 40 vol% for DF and HCF runs, respectively. Flow volumes vary between 1.5 × 10⁵ m³ (small-volume events) and 5 × 10⁵ m³ (large-volume events). Following Charbonnier et al. (2018), the Manning roughness coefficients vary with channel geometry and roughness between 0.03 and 0.08, the laminar flow resistance is set at 400, and the extended Froude number at 0.9 and two different rheological parameters (α _i,j and β _i,j ) were used to simulate the differences in yield strengths and viscosities of DFs versus HCFs (Table 7). Inundation maps for the DFs and HCFs are shown in Figures 6A,B, respectively.

FLO-2D model results (Figures 6A,B) highlight the contrasting behavior of HCFs compared to DFs. Runout distances of the three simulated DFs (Figure 6B) do not exceed 16 km, and flows are confined to the main river channels only, except for a flow overspill that occurred in a sharp bend located just before the J. Olaya bridge (Miraflores district) in the Qda. San Lazaro (dashed black rectangle area in Figure 6B). In contrast, runout distances of the three simulated HCFs (Figure 6A) vary between 17 km (small-volume event) and >20 km (moderate-to-large-volume events), and major overspills from both sides of the San Lazaro and Huarangal-MM channels occurred between 15 and 20 km runout distances (dashed black rectangle areas in Figures 6A,B), with the moderate- and large-volume scenarios (volume ≥3.5 × 10⁵ m³) flooding the neighboring streets next to both Quebradas with <4 m thick water-rich HCFs. Along the San Lazaro channel section passing through the city of Arequipa, the major HCF overspills (up to 300 m in width) occurred at runout distances between 16.5 and 17.5 km (overbank area two in Figure 6A), in areas of channel constriction, first at the El Progreso bridge and in the Selva Alegre and Historic Center districts. Along the Huarangal-MM channel in Arequipa, a minor overspill (<130 m width) occurred first along an area of channel constriction due to the presence of the Tupac Amaru bridge (overbank area one in Figure 6A). Then, a major overspill occurred ∼1 km downstream of the Santa Rosa bridge, flooding the streets on the NW area up to 800 m away from the Qda. Huarangal-MM (overbank area 2 in Figure 6A). Flow overspills located at the distal end of both Quebradas (before reaching the Río Chili in the San Lazaro and in the Explanada Dolores (a meeting square for public events) in the Huarangal-MM, are due to the poor resolution of the topography on the DSM in these areas, mostly due to the presence of trees.

These first FLO-2D modeling results at Misti highlight that while DFs are mostly confined in the main river channels, impacted areas of moderate-to-large-volume (>3.5 × 10⁵ m³) HCFs in the city of Arequipa can be extensive and dangerous. HCFs can spill over from the shallow and narrow riverbeds and inundate highly populated streets and critical city roads and bridges.

Uncertainties and limitations associated with tephra-fall scenarios and modeling have already been addressed in Sandri et al. (2014). The authors considered the natural variability of the eruptive process in terms of eruptive sizes and vent location. Many possible scenarios were combined, each one weighted with its own probability of occurrence. This allowed the authors to quantify how the uncertainty in the eruptive size and vent location affected the probability for different areas around the volcano to be impacted in a given time window: the associated hazard. Herein, uncertainties are linked to the choice of two hazard metrics: the tephra load of 100 and 250 kg/m². The first choice was driven by the poor-to-moderate resistance of the majority of the roofs in the surveyed neighborhoods along the San Lazaro ravine. The majority of flat roofs and the paucity of pitched roofs across Arequipa reduces the roof resistance. As the probability of ashfall toward the city increases during the rainy season, one of the study limitations stems from the fact that the rainfall was not included in the tephra-fall scenario. Wet tephra will have a considerable effect on the metal sheet and hybrid roofs if the latter are not cleaned during the eruption.

Input parameters of the three numerical codes cannot be used for benchmarking purposes. However, we can compare outputs (extent, paths, overbank areas, velocity, and flow depth) with characteristics reported from the most recent flow events. All three codes did a good job in reproducing observed flows, particularly recent HCFs across two basin ranges in Arequipa, which are the most common events from the database 1915–2020. Simulation runs show the expected path and extent along the entrenched ravine channels, particularly Qda. San Lazaro. All runs show overbank areas down valley at sharp bends, break in slope between volcaniclastic fans, and upstream of narrow bridges. In the Qda. Huarangal-MM, the most voluminous DFs and HCFs only can reach the densely populated neighborhoods along the middle and lower reaches. This is due to the long (16 km) path and gentle slope down valley, and Titan2F simulated flows are shorter than other flows simulated with VolcFlow and FLO-2D. The shorter inundated areas may be due to at least three factors: 1) the Mohr–Coulomb model prevents the flow from propagating on slopes below the lowest basal friction angles chosen as input; 2) viscosity, allowing the flow to propagate on gentle slopes, is not considered in Titan2F; and 3) the initiation process is not a mass flux, but a pile of material which collapses due to gravity.

Titan2F-simulated flows display a narrow extent, shorter paths, particularly in the Qda. Huarangal-MM, and a limited number of overbanks with respect to documented channel sinuosity and constrictions. VolcFlow- and FLO-2D-simulated flows exhibit longer and/or wider paths, higher velocity, and overbanks in the same areas at bends and down valley, which pinpoint reaches where attention must be drawn for mitigation measures to be implemented (containment wall, check dams, or enlargement of bridges). VolcFlow and FLO-2D reproduce overbank at the Olaya bridge bend in the middle reach of San Lazaro that does not appear with Titan2F runs.

The outputs of all three codes included flow depth, velocities, and dynamic pressure. Ranges of flow depth (0.1–4.6 m) and velocities (1–9 m/s) are comparable between codes and real flows and with simulated flows in the literature, whereas the distribution of dynamic pressure differs between the code outputs. Dynamic pressure has never been measured during or immediately after flows in Arequipa. By considering the same middle reach of the Qda. San Lazaro, the range of dynamic pressure is narrower with Titan2F runs (0–50 kPa), widens to 75 kPa with VolcFlow runs, and expands beyond 100 kPa with FLO-2D runs. VolcFlow- and FLO-2F-simulated HCF and DF display a wider area impacted by dynamic pressure on both sides of the ravine channel, extending onto the low terraces or banks on which low-quality buildings are located (Figures 5C, 7). Usually, HCFs exert a weaker dynamic pressure range than DFs, but FLO-2D-simulated HCFs show a higher dynamic pressure range (Figure 7).

By using FLO-2D, density is automatically computed according to the specific weight of the flow, which varies progressively as the flow propagates down valley. Density values can be much higher in HCFs than DFs in channels; in contrast, density is an initial parameter within Titan2F and VolcFlow and does not change down valley. In terms of rheology, HCFs mostly displayed the highest dynamic pressures acting on parallel walls (Mead et al., 2017). The higher density (compared to Newtonian flows) is responsible for the larger dynamic pressures (Jenkins et al., 2015). This effect is moderated by the yield strength of HCFs, causing the velocity to be lower than Newtonian flows near perpendicular walls. DF pressures are much lower than both Newtonian and HCFs as the yield strength and dilatant rheology components limit overbank flow velocities.

From estimated values of lahar impact pressure in Sarno, Italy, Zanchetta et al. (2004) argued that the complete destruction of a building required a pressure exceeding 90 kPa. Below this value, structures were generally heavily (35–90 kPa) or moderately (<35 kPa) damaged. The required pressure to break a large single-glazed window was as low as 1–2 kPa, and it was little more than 3 kPa for a wooden door. Therefore, an important point is that openings are the weakest points of buildings and are a critical parameter to quantify in buildings surveys. The dynamic pressure threshold of 35 kPa (Zanchetta et al., 2004) for moderate damage would be the case along both ravines for flows simulated with Titan2F and VolcFlow. The dynamic pressure threshold of 90 kPa for severe damage has been reached along both channels, as shown by HCFs and DFs simulated with FLO-2D. Heavy damage triggered by dynamic pressure between 35 and 90 kPa can be attained in the channel center and along sharp curves and/or narrow sections (e.g., upstream and near the J Olaya bridge; Figures 5C, 7). In all cases, the high (30%–50%) proportion of openings, the prevailing single-glazed windows in surveyed buildings in Arequipa, represent the weakest point for HCF and DF impacts in Arequipa.

The evaluation of tephra load on roofs has shown that at least 50% of building roofs of classes A0, A, and B are highly vulnerable to tephra load of 250 kg/m², a force of 2 kN/m². This is equivalent to a 10-cm-thick tephra that would fall in the rainy season over the city. As future El Misti eruptions would likely produce 10- to 50-cm-thick tephra-fall deposits, the probabilities of roof collapse would affect at least four among the nine categories of roofs.

Considering the cost of roof collapse, the low quality of roofs, and the low-income households in the surveyed neighborhoods, retrofitting and construction need to be undertaken by the municipality, the regional government, and NGOs. Specific construction rules and prevention procedures should be implemented as soon as El Misti awakens and before an imminent eruption.

Damage from water-rich flows is induced by a combination of three principal forces: 1) hydrodynamic pressure, 2) hydrostatic pressure, and 3) the collisional force of boulders (Thouret et al., 2020, and references therein). Hydrodynamic pressure exerted by the flow front/head, having horizontal and vertical components, is the kinetic energy per unit volume of the flow, which changes with flow density. Hydrostatic pressure depends on the depth and weight of the wet flow until it solidifies and acts toward the end of the event once the deposit is stabilized. Herein, we will not consider boulders, which act like missiles, hitting structures, and obstacles. Hydrodynamic pressure is proportional to the square of the frontal velocity as follows: P d y = k ρ d f v 2 , (5) where Pdy is the dynamic pressure; k is a constant, which depends on the density of the granular material, flow dynamics, and flow homogeneity and constituents; ρdf is the mean density of the mixture; and v is the velocity of the flow front.

Zeng et al. (2015) distinguished three principal damage categories for reinforced concrete and masonry structures by linking flow processes with hazard intensity: 1) inundation or burial, including damage to the ground floor or external walls of a building together with debris entering rooms, but without significant damage to structural components; 2) serious structural damage owing to flow front impact or boulder impacts triggering the failure of single structural elements, or the collapse of the whole structure; and 3) undercutting where soil erosion and/or liquefaction deforms foundations and subsequently buildings tilt or collapse. Faella and Nigro (2003) and Nigro and Faella (2008) stated that the severity of LH/DF impacts largely depended on the orientation and structural types of edifices and the kinetic energy of flows (a function of velocity, density, and Froude number).

From the evaluation of fragility curves by Mead et al. (2017) in the Qda. Dahlia, a tributary of Qda. Huarangal, we elaborated fragility curves for 425 buildings surveyed in 2018 and 2019 along the Qda San Lazaro and Huarangal. Fragility curves (Supplementary Figure S4 taken from Mead et al., 2017) were adapted and expanded to all surveyed buildings of classes A0, A, and B along both quebradas. Fragility curves based on hydrodynamic pressure and flow depth indicate to which extent buildings of classes A0, A, and B with wall bricks of 25 and 15 cm width would be affected by the ultimate bending moment before collapsing. The ultimate bending moment M u = ( f t + f d ) w b 2 / 6 , (6) where ft is the tensile strength of the masonry wall, fd is the design compressive stress acting on the wall, w is the width of the wall facing the flow, and b is the thickness of the wall.

Brick walls (and “bloquetas,” i.e., lapilli-rich rubble stones) represent about 60% of the buildings in these areas. Figures 5C, 7 and Supplementary Figure S4 suggest that moderate dynamic pressure in the range of 20–30 kPa and flow depth <1.5 m would suffice to bend the 15-cm-thick and 25-cm-thick brick walls, respectively.

Figure 8A, adapted and expanded from Mead et al. (2017), shows the loss fractions of buildings having brick walls (15 and 25 cm thick, respectively) due to HCF and DF with flow rates as low as 25–100 m³/s. Three out of four flow rate values (25, 50, and 75 m³/s) considered in the loss evaluation of buildings of classes A0, A, and B correspond to the 5- to 10-year return period of flow events in Arequipa (Table 6), to be compared with flow rates in the order of 125 m³/s during the 8 February 2013 HCF in Qda. Venezuela, which represents a 50-year return period flow according to Table 6. The loss proportion of all class buildings is 30% higher for the brick wall of 15 cm thickness when HCFs are considered with respect to DFs. For the 15-cm-thick brick walls, the loss proportion of class A0 buildings is high (75%–100%) in all HCF and DF cases, the loss proportion of class A buildings is relatively high (85%–100%) in case of HCF but less in case of DF (10–85%). The loss is relatively minor (0%–40%) for class B buildings in case of DF except for flow rates ≥100 m³/s and substantial (45%–100%) in case of HCF. As a result, a substantial loss of low-quality buildings on both sides of the ravines in Arequipa would result from HCFs within the range of 5- to 10-year return period of low-to-moderate (25–75 m³/s) flow rates.

The monetary cost of loss fractions of buildings has been estimated according to flow scenarios to contribute to planning emergency procedures and civil defense works. The range of cost has been estimated for each structural type of building (Table 3; Figure 8C). Evaluating cost is challenging in Arequipa because financial issues are a sensitive subject, particularly in city suburbs, where the declaration of assets is biased and tax evasion is omnipresent. Figure 8C displays the estimated cost range of construction/retrofitting of each building type, which is not the market value, depending on construction age and location, use, maintenance, and other factors. The estimated cost range (Figure 8C) is based on damage level proportions of building classes (Figure 8B). Class A0 buildings show the highest (16%–21%) damage levels in all categories (heavy, significant, and slight). Class A buildings display moderate (7%–13%) damage levels in almost all categories, while class B buildings exhibit low (0.5%–5%) damage levels. In terms of monetary cost, the situation is the opposite, as the most vulnerable class is low quality and “cheap” housing. Thus, the estimated cost of reconstruction or retrofitting of damaged class A0 buildings is quite low (0.7%–2.2%). The estimated cost of damaged class A buildings is moderate to substantial (9.7%–32.9%), while the cost of less damaged but regular-quality class B buildings covers a wide range between low (0.6%) and high values (34.1%), which is explained by the fact that a small proportion (<5%) of class B buildings would suffer severe or significant damage.

Should a prevention policy be implemented in the near future, we suggest that houses of classes A0 and A bordering the ravines could be relocated on high terraces at a minimum distance of 250 m away from the channel. Alternatively, the financial and logistical support of the Municipality, Regional Government, Ministry of Housing, and NGOs would allow the reconstruction of houses of classes A0 and A, for example, in the SE basin of Arequipa away from the ravines and El Misti (Municipalidad Provincial de Arequipa, 2011). A large proportion of severely and significantly damaged class A buildings should be reconstructed, while the slightly damaged proportion should be retrofitted. Class B buildings, probably slightly damaged in small proportion, can be retrofitted. Wall containment at sharp bends and enlargement of bridge arches in constricted channels are recommended where overbank flows were simulated and actually occurred in the past.