Examining the impact of lahars on buildings using numerical modelling

Abstract

Lahars are volcanic flows containing a mixture of fluid and sediment which have the potential to cause significant damage to buildings, critical infrastructure and human life. The extent of this damage is controlled by properties of the lahar, location of elements at risk and susceptibility of these elements to the lahar. Here we focus on understanding lahar-induced building damage. Quantification of building damage can be difficult due to the complexity of lahar behaviour (hazard), varying number and type of buildings exposed to the lahar (exposure) and the uncertain susceptibility of buildings to lahar impacts (vulnerability). In this paper, we quantify and examine the importance of lahar hazard, exposure and vulnerability in determining building damage with reference to a case study in the city of Arequipa, Peru. Numerical modelling is used to investigate lahar properties that are important in determining the inundation area and forces applied to buildings. Building vulnerability is quantified through the development of critical depth–pressure curves based on the ultimate bending moment of masonry structures. In the case study area, results suggest that building strength plays a minor role in determining overall building losses in comparison to the effects of building exposure and hydraulic characteristics of the lahar.

Introduction

Lahars, defined as gravity-driven flows containing a mixture of volcanic sediment and water (Vallance and Iverson, 2015), have caused severe damage to infrastructure and buildings (e.g. de Bélizal et al., 2013; Pierson et al., 2013; Ettinger et al., 2015; Jenkins et al., 2015) in addition to being responsible for a large proportion of volcanic fatalities (Auker et al., 2013). Assessing the extent of potential lahar damage can be difficult due to the complexity of flow behaviour, varying number of elements (e.g. buildings and bridges) exposed to lahars and a lack of knowledge on the structural capacity of these elements to withstand damage-causing components of the lahar flow. Using the common definitions of Varnes (1984), we define the damaging components of lahar flow (e.g. velocity, depth and pressure) as the hazard, environmental characteristics of exposed elements (e.g. building locations and orientations) as the exposure and the ability of exposed elements to withstand the hazard (e.g. building strength) as vulnerability. Lahar-induced damage is controlled by the interactions of these factors; however, the importance of each component can vary. Here we focus on quantifying and examining the role of hazard, exposure and vulnerability in determining lahar-induced building damage.

Post-event field assessments of building damage can elicit information relating lahar hazard to structural damage. However, these assessments tend to only record information on substantial damage, are affected by terrain changes during the event which alter exposure and often rely on a priori assumptions of building strength and vulnerability (Ettinger et al., 2015). Pre-event assessments are affected by the lack of reliable hazard intensity measures (van Westen et al., 2006; Ettinger et al., 2015), differences in spatial and temporal scales, uncertainty surrounding site-specific lahar triggers (Di Baldassarre and Montanari, 2009) and a lack of structural information on building stock (Ettinger et al., 2015). These issues are reflected in the relative lack of studies on hazard impact in urban areas (Jenkins et al., 2015) and often result in a reliance on expert judgement to develop vulnerability models for lahars and flash floods (Ettinger et al., 2015).

The physical vulnerability of buildings, defined as the susceptibly of a building to damage with respect to the hazard (Künzler et al., 2012), is a function of building characteristics such as size, shape, age, construction materials, structural integrity, maintenance and build quality (Martelli, 2011; Künzler et al., 2012; Ettinger et al., 2015). Information on these building properties is often lacking and hard to collect on a large scale. This commonly leads to the simplification of vulnerability into a measure that can provide a relative indication of vulnerability and consequent damage (Künzler et al., 2012). Studies simplifying vulnerability as a relative index use a combination of qualitative and quantitative metrics obtained through building surveys, interpretation of remote sensing data and GIS techniques to map and analyse vulnerability on a large scale (e.g. Lavigne, 1999; Künzler et al., 2012; Galderisi et al., 2013; Thouret et al., 2013, 2014; Ettinger et al., 2015). These methods can be applied to understand and highlight spatial patterns in vulnerability; however, as a relative measure, they cannot provide guidance on absolute damage for any specific event.

A direct estimation of damage caused by specific events requires quantified relationships describing a building's response to the hazard. Buildings can be damaged through a number of mechanisms including (i) direct damage resulting from static and dynamic forces imposed by the flow, (ii) damage to foundations through erosion and scour, (iii) buoyancy effects of the flow causing structures to float, (iv) direct damage from larger debris (missiles) within the flow and (v) indirect damage caused by chemical and biological actions such as seeping induced weakness of mortar (Kelman and Spence, 2004). All these actions, apart from chemical and biological effects, are related directly to lahar depth, velocity or a combination of depth and velocity. Thus, a common approach in determining building damage thresholds for a particular building type is to relate damage to hazard intensity measures of depth and/or velocity (e.g. Zanchetta et al., 2004; Custer and Nishijima, 2015; Jenkins et al., 2015). However, building typologies are affected by socio-economic, cultural and institutional conditions (Künzler et al., 2012), and hazard intensities (flow depth and velocity) are affected by building environmental factors such as local elevation, distance from main channels and orientation (Thouret et al., 2014). This results in complex interactions between hazard, exposure and vulnerability. These issues cause direct vulnerability relationships to be site-specific and require detailed investigation of the regions at risk to examine the relative effects and role of hazard, exposure and vulnerability on building loss.

We attempt to quantify and examine the components that determine building damage in a small area within the city of Arequipa, Peru. A relative index of vulnerability on a city-block scale was developed for Arequipa in Thouret et al. (2013, 2014). The studies by Thouret et al. (2013, 2014) highlighted two groups of vulnerability indicators (building characteristics and the physical setting) that play an important role in determining vulnerability within Arequipa. In this study, we separate the building characteristics (i.e. vulnerability) from the physical setting (i.e. exposure) to examine the effect of hazard (flow rate and rheology), exposure (building orientation and location) and vulnerability (building type) components on building loss within Arequipa. Physical vulnerability of buildings is explicitly separated from exposure through the development of a building damage model that is dependent on flow velocity, depth and sediment concentration. Simulations of lahar flow using smoothed particle hydrodynamics are used to examine how flow characteristics and the physical setting of city blocks affects forces on buildings and the consequent damage. While, for reasons explained earlier, damage functions presented here are necessarily specific to Arequipa, the hazard modelling approach and vulnerability model development are described in detail to support risk assessment in other regions affected by lahars.

Case study: Quebrada Dahlia, Arequipa, Peru

The central business district of Arequipa, the second largest city in Peru, is situated 17 km south-west of the summit of El Misti (Fig. 1), a steep stratovolcano with a history of explosive eruptions. Rapid population growth since 1960 has resulted in an expansion of the city towards the ring plain and steep slopes of El Misti (Thouret et al., 2013). Arequipa is drained by several ravines (locally called quebradas or torrenteras), shown in Fig. 1, that have been shaped by lahars and floods originating from the volcano on volcanoclastic fans north-east of the city. These quebradas are normally dry but carry water sporadically during the December to March rainy season (Vargas Franco et al., 2010; Martelli, 2011; Thouret et al., 2013; Sandri et al., 2014). Flash floods and hyperconcentrated flows occur relatively frequently in the quebradas, with return periods between 2 and 10 years (Vargas Franco et al., 2010; Thouret et al., 2013). Previous studies of lahar hazard and vulnerability for Arequipa identified seven alluvial terraces (T0, T1, T1${}^{\prime }$, T2, T2${}^{\prime }$, T3 and T4) based on stratigraphy and local elevation above the quebrada and the Rio Chili valley (Martelli, 2011; Thouret et al., 2013, 2014). The likelihood of inundation by a lahar or flash flood decreases with each terrace. Terrace levels T0 and T1 (up to 3 m above the quebrada) are frequently flooded (approx. every 2 to 10 years). The higher terraces (T1${}^{\prime }$ to T2${}^{\prime }$, 3 to 10 m above the quebrada) are rarely flooded (estimated 20 to 100 years) and the highest terraces (T3 and T4) are only likely to be inundated by lahars linked to large eruptions (Thouret et al., 2013, 2014). A city-wide vulnerability study by Thouret et al. (2014) identified that the city blocks most vulnerable to flash floods and lahars were on the lower terraces and typically within 100 m of a quebrada.

Location of Arequipa in relation to the El Misti volcano, showing the main quebradas and the location of the Quebrada Dahlia study area.

[Figure omitted. See PDF]

Individual building type and vulnerability classes for each block in the Quebrada Dahlia study area. Block ID increases from north to south.

Block	ID	Type	Structural class		Structural class	Type	ID	Block
	1	1A	A		A0	2A	1
	2	1A	A0		A0	2A	2	East 1
	3	4	B		A	3	3
West 1	4	4	A0		A0	1B	1
	5	2B	A0		A	3	2	East 2
	6	3	A		A0	1B	3
	7	4	B		A0	1A	4
	1	4	B		B	4	1
	2	4	B		A	3	2
West 2	3	4	B		A	3	3	East 3
	4	4	B
	5	3	A

To build on this study and investigate the vulnerability of the quebrada channel and banks in detail, simultaneous photogrammetry and building surveys were undertaken along short sections (approximately 200 m) of several quebradas during September 2013. Here we focus on one 150 m-long section of Quebrada Dahlia to examine lahar hazard and building damage. Quebrada Dahlia is a small tributary of Quebrada Mariano Melgar-Huarangal (Fig. 1), which is situated in the Mariano Melgar District on the north-easternmost fan of Arequipa, shown in detail in Fig. 2. The case study area was chosen for the following reasons:

• The quebrada channel is reasonably straight, reducing the effect of bends in the watercourse on lahar dynamics.

• Building quality varies from well-built reinforced masonry buildings to makeshift structures with little to no mortar. This allows for an investigation of the effects of building quality on damage caused by lahars.

• All buildings are situated on the lowest terraces (T0–T1${}^{\prime }$, 1–5 m above the channel), meaning they may be affected by even the smallest events identified in Vargas Franco et al. (2010) and Thouret et al. (2013).

Building types and simplified structural classes from Thouret et al. (2014).

Typology	Building description	Simplified
		structural class
1A	Unreinforced masonry of lapilli, ignimbrite or terracotta with no roof support structure (i.e. metal sheet roof).	A0
1B		A0
2A		A0
2B		A0
3	Terracotta masonry with reinforced concrete roof.	A
4	Terracotta masonry with reinforced concrete frame and roof.	B
5	Historical ignimbrite building with mortar.	A
6A	Ignimbrite masonry with reinforced concrete elements or modifications.	B
6B		B
6C		B

A three-dimensional reconstruction of the terrain and buildings along Quebrada Dahlia was created using the photogrammetry method described in Mead et al. (2015). The surface reconstruction, shown in Fig. 2b, contained 1.4 million points with a surface density of between 150 and 750 points per m${}^{\mathrm{2}}$. A GNSS-D survey undertaken in October 2014 enabled georeferencing of the terrain reconstruction for possible inclusion in future GIS applications. The surface reconstruction was smoothed and reduced to create the lower-resolution terrain model shown in Fig. 2c to be used in the lahar simulations. The terrain model in Fig. 2c contains 22 buildings identified during building surveys in 2013. Streets and the quebrada (shown in Fig. 2a) separate these buildings into five city blocks (labelled in Fig. 2c), referred to hereafter as “blocks”. The typology of each building was characterised through surveys undertaken in 2013 following the approach of Thouret et al. (2014). In this approach, buildings are classified as one of eight structural types (1A–8C) based on a visual inspection to determine construction material, roof type and structural support (see Table 1). These types were then grouped into larger simplified structural classes. Using this building classification system, the study area contains eight class A0 buildings, seven class A buildings and seven class B buildings (see Table 2 for a description of building types and structural classes).

Overview of the Quebrada Dahlia study area, Arequipa, (a) aerial image with black outline showing study area, dashed outline showing channel banks and transparent lines showing streets in the area, (b) photogrammetric reconstruction of the surface and (c) individual buildings and building blocks identified from building surveys.

[Figure omitted. See PDF]

Developing building vulnerability relationships

Buildings and infrastructure can be damaged through a variety of mechanisms brought upon by the actions of a lahar. Here, as in most other studies of lahar damage (Zanchetta et al., 2004; Toyos et al., 2008; Ettinger et al., 2015; Jenkins et al., 2015), we focus on the direct damage resulting from hydrostatic and hydrodynamic forces applied to buildings. We regard these actions as the most important, although scour and large debris missiles within the flow can also cause significant damage (Jenkins et al., 2015). Scour and debris actions are neglected here as they are currently too difficult to predict and incorporate into large-scale loss analyses (Kelman and Spence, 2004), particularly in regions with limited hazard and exposure information.

The building stock within Arequipa is characterised mostly by masonry structures of varying quality, with some reinforced concrete structures (Thouret et al., 2014). Therefore, we develop vulnerability relationships that are primarily focused on masonry buildings. A structural failure model similar to those employed by Roos (2003), Custer and Nishijima (2015) and Zeng et al. (2015) is implemented. In these models, masonry walls are presumed to fail when the lateral pressure imposed on the wall results in a bending moment or shear force greater than the wall's calculated ultimate bending moment or ultimate shear force.

The ultimate bending moment ($M_{\mathrm{u}})$ is calculated using the following equation (Roos, 2003): $$M_{\mathrm{u}}=(f_{\mathrm{t}}+f_{\mathrm{d}}){\displaystyle \frac{w{b}^{\mathrm{2}}}{\mathrm{6}}},$$ where $f_{\mathrm{t}}$ is the tensile strength of the masonry wall, $f_{\mathrm{d}}$ 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, which is assumed equal to the brick width. Tensile strength and design compressive stresses for buildings in Arequipa are calculated using the approach specified in Australian Standard (AS) 3700–2011 and summarised in Appendix A in addition to an approach calculating ultimate shear force. Preliminary investigations using these two approaches suggested that the force required to overcome the ultimate moment was consistently lower than the force required to overcome the ultimate shear force. Therefore, we chose to focus the remainder of this study on the ultimate bending moment only.

The use of a foreign standard to calculate the ultimate moment should still be valid for the study area if construction material properties from Arequipa are used as inputs. However, some specifications and assumptions of the standard may not be relevant. Notably, observations during the building survey suggest that construction methods and conformity to specifications within the standard differs substantially to those specified in AS 3700–2011. This difference will influence ultimate bending moments, particularly those for low-quality unreinforced building types (i.e. building types 1A–2B) due to the makeshift nature of construction. For these classes, calculated bending moments will represent a best case scenario in which masonry unit strength and quality have not been compromised by construction methods.

Critical depth–pressure curves

The range of design compressive stress for each building typology is shown in Fig. 3. The range was obtained by calculating the design compressive stress for every configuration of masonry compressive strength ($f_{\mathrm{c}})$, wall thickness ($b)$ and thickness coefficient ($k_{\mathrm{t}})$ in Appendix A. Buildings with reinforced frames (types 3, 4 and 6) can withstand much greater compressive stresses than non-reinforced buildings (types 1, 2, 5). The wall thickness has a large effect on building strength, which is consistent with observations by Jenkins et al. (2015). Notably, the design compressive stresses are similar for building types that share the same simplified structural class identified in Thouret et al. (2014), based on the structural classes of Zuccaro et al. (2008). Given these similarities, we also use the simplified structural classes (A0, types 1A–2B; A, types 3 and 5; and B, types 4 and 6A–6C, see Table 2) from Thouret et al. (2014).

Range of design compressive stress for building types 1A–6C defined in Thouret et al. (2014). Compressive stress capacity was calculated for every configuration of compressive strength ($f_{\mathrm{c}})$, bedded area ($A_{\mathrm{b}})$ and thickness coefficient ($k_{\mathrm{t}})$ at brick widths ($b)$ of 150 and 250 mm.

[Figure omitted. See PDF]

Critical depth and dynamic pressures for failure of structural classes A0, A and B for brick widths of 0.15 m (top) and 0.25 m (bottom). Shading of the lines indicate flow type and density, dotted lines and dashed lines represent the minimum and maximum forces required. Densities are for a Newtonian flow (NF, $\mathit{\rho }=$ 1000 kg m${}^{-\mathrm{3}})$, hyperconcentrated flow (HCF, $\mathit{\rho }=$ 1500 kg m${}^{-\mathrm{3}})$ and debris flow (DF, $\mathit{\rho }=$ 1915 kg m${}^{-\mathrm{3}})$.

[Figure omitted. See PDF]

The critical depth (i.e. hydrostatic pressure) and dynamic pressure required to overcome the ultimate bending moment (Eq. 1) for each structural class is shown in Fig. 4. These curves assume that both hydrostatic and dynamic pressure act on walls. Other studies (e.g. Jenkins et al., 2015) assume only dynamic pressure acts on walls due to an equalisation of lahar depths on the inside and outside of buildings. This equalisation can take a reasonable amount of time, which is likely to be much longer than the simulation duration studied here (see following sections). The curves in Fig. 4 indicate the structural limit of each class; combinations of depth and pressure that fall above the curves indicate an applied moment greater than the building can withstand. Conversely, combinations of depth and pressure that fall below the curves indicate an applied moment less than the maximum the building can withstand. Figure 4 shows that the critical depth decreases with the density of flows as the hydrostatic pressure gradient is much larger for sediment-rich lahars. The critical depths and pressures are also affected by the structural class, with A0 structures being much less resilient than A and B structures. However, wall thickness has the most dominant effect on determining the strength of buildings. Wider walls increase the section modulus ($w{b}^{\mathrm{2}}/\mathrm{6}$ in Eq. 1), resulting in stiffer walls that also have a higher compressive stress capacity.

Lahar numerical modelling and results

Lahar rheology and implementation in smoothed particle hydrodynamics

Lahar flow behaviour varies depending on the sediment concentration and composition of the flow. At very low concentrations of sediment, lahars will flow in a similar manner to water. At higher concentrations, interactions between the sediment and water cause a non-linear response to stresses applied to the flow. This non-linearity in the stress–strain relationship requires the use of rheology models that capture both the linear (i.e. water-like, called Newtonian) and non-linear (called non-Newtonian) shear response. Here we implement a generalised quadratic rheology model to simulate lahar flows along the case study area. The quadratic rheology model can be expressed as follows (Julien and Lan, 1991; O'Brien et al., 1993; Jan and Shen, 1997): $$\mathit{\tau }={\mathit{\tau }}_y+\mathit{\mu }\dot{\mathit{\gamma }}+\mathit{\alpha }{\dot{\mathit{\gamma }}}^{\mathrm{2}},$$ where $\mathit{\tau }$ is the shear stress, ${\mathit{\tau }}_y$ the yield strength, $\mathit{\mu }$ the viscosity, $\dot{\mathit{\gamma }}$ the shear rate, and $\mathit{\alpha }$ is the turbulent-dispersive parameter, a coefficient that combines the effects of turbulence and dispersive stresses caused by sediment collisions. This model follows the general form of the Herschel–Bulkley equation commonly used to describe non-Newtonian lahar behaviour (Manville et al., 2013).

Commonly used lahar models such as the Pitman and Le (2005) model in Titan2D (Patra et al., 2005) or laharZ (Iverson et al., 1998) are able to delineate hazard zones or lahar inundation areas on a large scale. However, the reduced dimensions of these models (e.g. through depth-averaging in Titan2D) means they are unsuitable for the detailed modelling of lahar flow in urban environments required for this study. Instead, we implement the quadratic rheology model using three-dimensional smoothed particle hydrodynamics (SPH) to simulate lahar flows along the case study area. SPH is a Lagrangian method that tracks the physical motion of interpolation points (commonly referred to as particles) through space. It is well suited to modelling free surface fluid flows, predicting and tracking the motion of dynamic objects within the flow (e.g. Cleary et al., 2012, 2015; Prakash et al., 2014) and modelling complex flooding scenarios involving interactions with buildings (e.g. Mead et al., 2015). The SPH method used here is described in Cleary and Prakash (2004) and Prakash et al. (2014). Non-Newtonian lahar rheology was implemented in SPH using an apparent Newtonian viscosity ($\mathit{\eta })$. Assuming the fluid is isotropic, constitutive equations for rheology can be written as a generalised Newtonian fluid in terms of the apparent viscosity: $$\mathit{\tau }=\mathit{\eta }\dot{\mathit{\gamma }}.$$ When the apparent viscosity is constant the fluid is Newtonian with a viscosity of $\mathit{\eta }$. Non-Newtonian fluids can be modelled using Eq. (3) by developing relationships for $\mathit{\eta }$ based on constitutive equations (Mitsoulis, 2007). Using this approach, the apparent viscosity for the quadratic rheology is $$\mathit{\eta }={\displaystyle \frac{{\mathit{\tau }}_y}{\dot{\mathit{\gamma }}}}+\mathit{\mu }+\mathit{\alpha }\dot{\mathit{\gamma }}.$$ To reduce computational time we use the viscosity regularisation approach of Papanastasiou (1987), described in Mitsoulis (2007) and Minatti and Paris (2015). Regularisation is required as the apparent viscosity approaches infinity at low strain rates when using Eq. (4), reducing the simulation time step and significantly increasing computational cost. Using the Papanastasiou (1987) approach, the regularised viscosity used in simulations is $$\widehat{\mathit{\eta }}={\displaystyle \frac{{\mathit{\tau }}_y}{\dot{\mathit{\gamma }}}}\left(\mathrm{1},-,e^{-c\dot{\mathit{\gamma }}}\right)+\mathit{\mu }+\mathit{\alpha }\dot{\mathit{\gamma }},$$ where $c$ is the viscosity scaling parameter. Larger values of $c$ result in a better approximation of the constitutive equation (Eq. 4), while smaller values result in smaller apparent viscosities and larger simulation time steps. Here we set $c=\mathrm{200}$, a value which yielded the best balance between simulation speed and accuracy in validation simulations that compared flow down an inclined plane with analytical solutions.

Lahar simulations

Static and dynamic pressures acting on the buildings in the Quebrada Dahlia study area were determined for 12 different inundation scenarios. We use the same SPH particle spacing (12.5 cm) of previous simulations by Mead et al. (2015). This resolution provided the best balance between computational time and resolution of fine-scale features that can affect the flows. Inundation scenarios were designed to explore a wide range of flow types and velocities in order to investigate the effect of rheology and velocities on flow dynamics and forces exerted on buildings. Simulations were run for three different flow types (Newtonian, hyperconcentrated streamflow and debris flow) at constant flow rates of 25, 50, 75 and 100 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$.

Density, particle concentration and rheology coefficients for hyperconcentrated streamflow and debris flow simulations, taken from Govier et al. (1957), Julien and Lan (1991).

Flow type	Density	Particle concentration	Yield strength	Viscosity	Dispersive stress
	(kg m${}^{-\mathrm{3}})$	by volume (%)	(${\mathit{\tau }}_y$, Pa)	($\mathit{\mu }$, Pa s)	coefficient ($\mathit{\alpha })$
Hyperconcentrated	1500	30.3	0.94	0.0137	1.28 $\times$ 10${}^{-\mathrm{5}}$
streamflow
Fine-grained, matrix	1915	55.5	0.672	0.0485	0.00224
-supported debris flow

The flow rates were chosen to produce scenarios ranging from minimal (25 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}})$ to extreme (100 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}})$ overbank flooding. The ratio between inertial and gravitational forces, expressed through the Froude number, was kept below 1 (subcritical flow) for each flow rate by varying the inflow area. Froude number consistency was used here as inertial and gravitational forces are dominant controls on environmental flows such as these. Flow types were selected to represent the characteristics of the most commonly occurring flows in Arequipa – flash flood, hyperconcentrated streamflow and fine-grained, matrix-supported debris flow (Thouret et al., 2013). Rheology of flash flood flows was considered to be completely Newtonian with a viscosity of water (i.e. ${\mathit{\tau }}_y,\mathit{\alpha }=\mathrm{0}$, $\mathit{\mu }=\mathrm{0.001}$ and density ($\mathit{\rho })=\mathrm{1000}$). Rheological parameters for hyperconcentrated and debris flows (Table 3) were chosen using the dimensionless ratio between dispersive and viscous stresses explained in Julien and Lan (1991). Values for yield strength (${\mathit{\tau }}_y)$, viscosity ($\mathit{\mu }$) and the turbulent-dispersive coefficient ($\mathit{\alpha })$ were taken from the experiments of Govier et al. (1957) and Bagnold (1954), reported in Julien and Lan (1991). For a hyperconcentrated streamflow, we presumed a particle concentration by volume ($C_{\mathrm{v}})$ of approximately 30 %, consisting mostly of finer particles, meaning viscous stresses are still relatively important. Debris flow scenarios were assumed to contain larger particles at a higher value $C_{\mathrm{v}}$ of approximately 55 %. The particle concentration acts to increase density, viscosity and the dispersive stress coefficient in hyperconcentrated and debris flow rheologies compared to a fully Newtonian water flow. The higher particle concentration of the debris flow (compared to a hyperconcentrated flow) also results in a much higher dispersive stress coefficient, meaning that dispersive stresses will have more importance in determining flow behaviour.

Computational cost limits the length of simulations to the first 45 s of lahar flow for each scenario. The flow was not established and constant by 45 s, so these simulations do not represent the forces exerted on buildings by a steady flow rate. Instead, the scenarios considered here are more representative of the higher velocity and depth surges or waves in a lahar.

Flow behaviour

Figure 5 displays snapshots of velocity and dynamic pressure magnitudes for each flow type at a flow rate of 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$. Snapshots were taken at 15 s intervals and dynamic pressure was calculated as $\mathit{\rho }{v}^{\mathrm{2}}/\mathrm{2}$, where $v$ is the velocity magnitude. Lahars mostly followed the developed channel of Quebrada Dahlia for the first 15 s before overtopping the bank and spreading outwards. Channel and overbank pressures and velocity profiles are similar for Newtonian and hyperconcentrated flows, but the velocity of overbank flow is much lower for the debris flow rheology. This lower velocity is presumably caused by increased friction in the debris flow due to the higher viscosity and dispersive coefficients. The dynamic pressure differs between each rheology as a result of the varied densities (and lower velocity for debris flows); however, the maximum pressure is still similar between rheologies as maximum velocities are mostly confined to the channel.

Evolution of dynamic pressure and velocity magnitudes for a 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ flow along Quebrada Dahlia for a Newtonian flow (NF), hyperconcentrated flow (HCF) and debris flow (DF). Arrows indicate northerly (green) and easterly (red) directions.

[Figure omitted. See PDF]

The highest dynamic pressures in Fig. 5 are present along the centre of the channel, with much lower pressures near the buildings. The velocity magnitude may therefore not accurately represent the pressure forces acting perpendicular to the walls of each building. The critical strength of a wall is determined from the forces acting normal (perpendicular) to the structure; therefore it is important to calculate dynamic pressure from velocity normal to the wall. The section of Quebrada Dahlia studied here runs in a north–south direction and the buildings have walls that are oriented either parallel or perpendicular to the channel, so an initial understanding of the perpendicular forces acting on walls can be interpreted from the north–south (N–S) and east–west (E–W) velocity components. Figure 6 shows the dynamic pressure calculated from directional velocity components at 40 s for a flow rate of 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$. Figure 6 shows a consistent pattern for all rheologies where the pressure is dominated by the streamwise (N–S) velocity. The pressure applied to walls facing the stream ($\sim$ E–W direction, third column of Fig. 6) is much lower than the pressure applied to walls perpendicular to the stream. Higher pressures for E–W velocities are observed along cross streets splitting each city block; however, the pressure that acts perpendicular (N–S) to these walls is minimal. These observations indicate that pressure calculated from the magnitude of velocity, which is often assumed to be acting perpendicular to walls (e.g. Zanchetta et al., 2004; Jenkins et al., 2015), can be much higher than actual pressure acting on walls and the use of velocity magnitudes could therefore lead to an overestimation of building damage.

Directional components of dynamic pressure for a 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ flow along Quebrada Dahlia for a Newtonian flow (NF), hyperconcentrated flow (HCF) and debris flow (DF). Maximum pressure is 25 kPa for magnitude and N–S pressures, 1 kPa for E–W pressure.

[Figure omitted. See PDF]

In order to accurately estimate normal forces on walls, we calculate pressures from the velocity normal to each block. The normal velocity of fluid near each block face (e.g. north and west faces of the East 3 block) is calculated using the dot product of simulated velocity vectors and the direction vector of the block face. This normal velocity ($v_{\mathrm{n}})$ is averaged across the face and used to calculate a “normal” pressure using $\mathit{\rho }{v}_{\mathrm{n}}^{\mathrm{2}}/\mathrm{2}$. Figure 7 compares the dynamic pressures calculated from velocity magnitude and normal velocity for the West 2 block (see Fig. 2c) for Newtonian, hyperconcentrated and debris flow types. The pressures are measured for walls oriented approximately parallel to the quebrada (labelled “parallel”) and north-facing walls that are oriented approximately perpendicular to the quebrada (labelled “perpendicular”). The pressures exerted on parallel walls by the normal velocity are up to five times lower than velocity magnitude pressures. The pressure applied to perpendicular walls also differs between normal velocity and magnitude, with the timing of peak pressure also affected. This further demonstrates the importance of considering normal velocity rather than velocity magnitude when estimating dynamic pressures (and consequently damage).

Comparison of mean pressure magnitude (grey lines) and mean normal pressure (black lines) on block West 2 in the parallel and perpendicular orientations for a 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ flow along Quebrada Dahlia.

[Figure omitted. See PDF]

Mean normal pressures applied to each city block in the perpendicular and parallel orientations for a 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ flow.

[Figure omitted. See PDF]

Critical depth–pressure curves for building classes A0, A and B subjected to Newtonian flow. Peak normal pressures and corresponding depths applied to each city block are plotted as points for each flow rate.

[Figure omitted. See PDF]

Critical depth–pressure curves for building classes A0, A and B subjected to a hyperconcentrated flow. Peak normal pressures and corresponding depths applied to each city block are plotted as points for each flow rate.

[Figure omitted. See PDF]

Critical depth–pressure curves for building classes A0, A and B subjected to a debris flow. Peak normal pressures and corresponding depths applied to each city block are plotted as points for each flow rate.

[Figure omitted. See PDF]

Pressures acting on each block in the study area, calculated using the technique explained in the previous paragraph, are shown in Fig. 8 for a flow rate of 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$. Blocks East 1 and West 1 do not have walls facing perpendicular to the flow and therefore have no pressures recorded in that orientation. The pressure for each block generally follows a similar pattern through time with a well-defined peak pressure and a lower, steady background pressure. The rise of pressure to its peak value and reduction to its background value occurs over the space of approximately 20 s for each block. This time frame is too short to allow for an equalisation of hydrostatic pressure between the inside and outside of buildings, confirming that both hydrostatic and dynamic pressures are acting on walls during lahar surges. The timing of the peak is delayed for downstream blocks and the magnitude of the peak for each block varies. The differences in peak pressure are caused by exposure effects such as orientation and elevation of each block relative to the quebrada. Walls facing perpendicular to the stream are generally exposed to higher dynamic pressures than parallel walls, but this effect appears to vary and could be dependent on cross-street elevations (cross streets leading away from Qda. Dahlia increase in elevation at different rates).

In terms of rheology, hyperconcentrated flows mostly displayed the highest dynamic pressures acting on parallel walls. The higher density (compared to Newtonian flows) is responsible for the larger dynamic pressures (see Jenkins et al., 2015). This effect is moderated by the yield strength of the hyperconcentrated flows which cause the velocity to be lower than Newtonian flows near perpendicular walls. Debris flow pressures are much lower than both Newtonian and hyperconcentrated flows as the yield strength and dilatant rheology components limit overbank flow velocities.

Application of critical depth–pressure curves

Depth at the maximum value of pressure acting on block walls for each scenario is used to determine whether individual buildings in the study area can withstand the bending moment applied by hydrostatic and dynamic pressure. Figures 9 to 11 plot the peak pressure and “surge depth” (depth at the time of peak pressure) for Newtonian, hyperconcentrated and debris flows alongside critical depth–pressure curves for vulnerability classes A0, A and B with a wall thickness of 150 mm (results for 250 mm wall thicknesses are provided in the Supplement). The hazard variables of flow rate and lahar rheology appear to influence building damage, although the size of the effect is difficult to determine since most scenarios place depth and pressure combinations well above the critical curves for each block. The flow depth, which affects hydrostatic pressure and bending moment location, generally increases with the flow rate while the dynamic pressure appears to be mostly controlled by the rheology in combination with flow rate. The forces applied to the West 2 block, containing one class A and 4 class B buildings, are lower than the other blocks. This is possibly due to the relative elevation and orientation of each block to the quebrada (i.e. exposure) affecting dynamic pressure and lahar depth. Debris flow scenarios at flow rates of 25, 50 and 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ indicate depths and pressures below the critical limit for this block's building classes.

The orientation of walls to the flow direction is another element of exposure that affects the normal pressure exerted on walls. In several scenarios, perpendicular walls are subjected to higher dynamic pressures and lower depths than parallel walls. However, this effect appears to be conditional to the rheology of the flow as the opposite is true for debris flow scenarios. These two effects demonstrate the importance of considering exposure elements separately to vulnerability.

The proportion of buildings with depths and pressures above the critical curve for each scenario is shown in Fig. 12 for 150 mm brick widths and Fig. 13 for 250 mm wall thicknesses. Assuming a binary damage state model where damage is complete for depths–pressure combinations above the curve, these proportions can be used to directly represent building loss. For the thinner walls, all class A0 buildings are above the curve for all scenarios apart from the 25 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ debris flow. The East 1 block is not inundated in this scenario, resulting in two undamaged class A0 buildings. Class A and B buildings are also mostly destroyed, with the exception of lower flow rate hyperconcentrated and debris flow scenarios where some blocks are on the edge of inundation and therefore subjected to much lower depth–pressure combinations. Slightly fewer building losses occur with larger brick widths (Fig. 13) as the larger section modulus results in a greater resistance to bending moments. However, most buildings are still destroyed in Newtonian and hyperconcentrated flow scenarios. An exception to this is the 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ Newtonian flow where the highest pressure on the East 1 block occurs early in the simulation when the surge depth is low, reducing the magnitude of hydrostatic pressure and lowering the size of the applied moment.

Building loss fraction for all flow scenarios where buildings are assumed to have a brick width of 0.15 m.

[Figure omitted. See PDF]

The building loss results indicate that class A0 buildings are most vulnerable, with class A buildings marginally stronger due to the roof support. Losses for type B buildings in this area are much lower; however, this appears to be more closely related to building exposure than structural strength as most type B buildings are in two blocks subjected to lower depth–pressure combinations for all scenarios. Overall, similarly to the observations in Jenkins et al. (2015), the data presented here suggest that building strength (i.e. the vulnerability component) has a minimal effect on losses, and building location (i.e. exposure) relative to flow rate and type (i.e. hazard) plays a much greater role.

Limitations and discussion

The losses shown in Figs. 12 and 13 are estimates based on several assumptions that, while necessary for the estimation of building loss, could limit the accuracy of results. Firstly, the depth–pressure curves are created using ultimate bending moments derived from a foreign standard and do not consider proportional losses, only assuming damage is complete for depths and pressures above the critical curve. Secondly, the flow scenarios modelled here are a subset of likely scenarios and do not replicate all damage-causing actions of lahar flow. Finally, the maximum total pressure was assumed to be the sum of hydrostatic and dynamic pressure and to occur when dynamic pressure was at its peak. These limitations are discussed and justified in this section to highlight areas of improvement necessary for robust, quantitative estimation of lahar damage and vulnerability.

Building loss fraction for all flow scenarios where buildings are assumed to have a brick width of 0.25 m.

[Figure omitted. See PDF]

Depth–pressure curves

The critical depth–pressure curve is the contour where the ratio of applied (pressure) moment equals the ultimate (failure) moment of a given masonry wall. The calculation of ultimate bending moments followed an Australian standard (AS 3700–2011). Although bending moment calculations are similar for all national standards and material properties from Arequipa were used as inputs, the standard inherently assumes conformance to construction and design standards. This is demonstrated through the assumption of a minimum mortar strength of 0.2 MPa. Page (1996) suggests this strength can be achieved with correct mortar composition and laying; however, lower strengths are possible if there is low conformance to design standards. The makeshift structures that characterise class A0 buildings are likely to have mortar bond strengths that are much lower than the implied minimum of 0.2 MPa. This mischaracterisation of mortar strengths for makeshift structures will result in an overestimation of building strength and critical depth–pressure curves. Additionally, the depth–pressure curves assume a binary damage state, where failure is total when the applied (pressure) moment equals the ultimate (failure) moment. This neglects incremental damage states that require building repair (e.g. to doors or windows) and can cause a reduction in the overall building strength.

Both of these assumptions will result in an underestimation of loss if violated. Most of the flow scenarios caused depths and pressures that exceeded critical curves by a large margin and resulted in an almost total loss; conversely, flow scenarios that did not result in total or near-total losses usually had depths and pressure values that were well below the critical curves. This suggests that these assumptions are not critical to the results shown here, but may be important to consider in other case study areas.

Flow scenarios

The 12 flow scenarios were chosen to understand the effect of hazard properties (flow rate and rheology) on total loss. These scenarios may not represent any specific lahar event for Qda. Dahlia. Rather, scenarios were chosen to be representative of the range of lahar rheologies and flow rates that can cause building damage in Arequipa. The chosen flows have similar characteristics to observed lahars and lahar deposits (Thouret et al., 2013) and are therefore reasonably representative of the lahar hazards expected in Arequipa.

However, damage caused by these hazards may not be representative as only the direct actions of hydrostatic and dynamic pressure were considered in this study. While direct actions are regarded as the most important source of damage, they are also favoured in risk assessment due to the large-scale predictability of hydrostatic and dynamic forces (Kelman and Spence, 2004). Damage is likely to also be caused by scour and large debris missiles within the flow (Jenkins et al., 2015). In particular, boulders are often carried by lahars at the flow front (Iverson, 1997; Doyle et al., 2011) and can lead to significant damage (e.g. Zeng et al., 2015). However, these actions are harder to predict and incorporate into large-scale loss analyses (Kelman and Spence, 2004). These unstudied actions are generally proportional to depth, pressure or velocity, indicating that there may be a relationship between the ratio of applied to ultimate moment and damage through other actions. Given that only direct actions are considered in this study, the curves likely form an upper bound to complete damage, and depth–pressure combinations below the curve may still result in complete building damage through other mechanisms.

Pressure actions

Both hydrostatic and dynamic pressures were considered in bending moment calculations. Slower increases in depth, buildings with many openings and the location of buildings relative to the channel can also affect the equalisation of lahar depths and reduce the effect of hydrostatic pressure. However, lahar depth would still be an important factor to consider in building damage estimation as it controls location of the bending moment and can cause damage through other actions (e.g. inundation damage, buoyancy, corrosion).

The applied depth at the time of maximum pressure was used here to create the depth–pressure combinations to determine building loss. This surge depth was not necessarily the maximum depth of the lahar during the simulation. Maximum depths generally occurred at later times in the simulations when hydrostatic pressure may have equalised inside and outside buildings. This assumption of surge depth was valid for most cases, although the losses for the 75 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ Newtonian flows indicate that this approach can be too simplistic at times. The complexity of lahar flows within urban environments with intricate geometry and obstacles similar to the case study area means that broad generalisations and assumptions about flow dynamics, such as the assumption of a surge depth, are often limited in their validity.

Discussion

The combination of pressures applied to each block in the study area created bending moments that, with few exceptions, were much higher than the maximum moment buildings could withstand. The limitations identified in previous sections generally overestimate building strength and resilience to lahars which would result in greater damage than predicted here. The estimated building losses (Figs. 12 and 13) therefore represent the minimum expected losses for each flow scenario with damage likely to be more severe due to additional damage actions (e.g. boulders impacting structures) and the overestimation of building quality, particularly for class A0 buildings.

When inundated, blocks in this study area are subjected to depths and pressures higher than the strongest structural class buildings can withstand. Specific improvements to reduce vulnerability, such as adding roof support and utilising reinforced frames comprised of equally spaced RC columns will increase the overall strength of buildings by reducing the slenderness ratio (Eq. A6). Wider masonry units (wall thickness) and stronger mortar joints will also increase the overall building strength by increasing wall stiffness and therefore resistance to bending moments. However, this increased structural strength appears to only reduce losses in very low flow rate scenarios where there is proportionally less inundation. This suggests that, while each component of risk has a role in determining overall building losses, the variability in individual losses appears to be predominantly caused by flow dynamics (i.e. lahar hazard) and building exposure (e.g. proportion of building types and orientation within blocks).

Conclusion

Development of fragility functions in the form of critical depth–pressure curves for building classes within Arequipa has helped to provide insight into possible building losses and their cause. Building vulnerability is largely controlled by social, cultural and institutional factors (Künzler et al., 2012), so the depth–pressure curves are necessarily specific to Arequipa building typologies. However, given sufficient data on building strength, depth–pressure curves can be generated through the same approach as in Appendix A and used to quantify masonry building loss in terms of flow depth and pressure in other regions.

Estimated building losses (Figs. 12 and 13) are caused by the intersection of lahar hazard (flow rate, flow type) with building exposure (location, proportion of building types and orientation within blocks) and vulnerability (building type and strength). The almost total simulated building loss for all scenarios indicates that substantial losses can be expected in the event of inundation. Furthermore, lahar depths and pressures obtained from simulations were much greater than most of the buildings in the study area could withstand, even if retrofitting to improve structural strength was undertaken. This suggests that, in this study area at least, exposure and lahar hazard have a larger role in determining building loss than vulnerability.

Source code and data used in this publication are available from the Zenodo repository 10.5281/zenodo.155144 (Mead, 2016).

Calculating ultimate bending moment and shear force

The ultimate bending moment ($M_{\mathrm{u}})$ and ultimate shear force ($V_{\mathrm{u}})$ are calculated using the following equations (Roos, 2003): $$\begin{array}{}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{u}}=(f_{\mathrm{t}}+f_{\mathrm{d}}){\displaystyle \frac{w{b}^{\mathrm{2}}}{\mathrm{6}}}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{u}}=f_{\mathrm{v}}wb,\end{array}$$ where $f_{\mathrm{t}}$ is the tensile strength of the masonry wall, $f_{\mathrm{d}}$ 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. The shear strength of the masonry wall ($f_{\mathrm{v}})$ is related to the tensile and compressive stress through (Roos, 2003) $${\displaystyle }{f}_{\mathrm{v}}=\mathrm{0.5}{f}_{\mathrm{t}}+\mathrm{0.5}{f}_{\mathrm{d}}.$$ The tensile strength is assumed to be 0.2 MPa as, according to AS 3700–2011, the tensile strength should be no greater than this value without testing. The wall thickness, $b$, is between 150 and 250 mm for terracotta bricks (Martelli, 2011) and is assumed to be similar for ignimbrite bricks observed in the study area. The design compressive stress, $f_{\mathrm{d}}$, can be determined by calculating the vertical forces (i.e. building weight) acting on the walls. This can be estimated from building properties such as number and weight of floors, weight of the masonry and building design (e.g. Roos, 2003). However, such detailed building data are lacking here and carries considerable uncertainty for a heterogeneous urban area with varied construction materials, building ages and designs such as Arequipa. Instead we use the design compressive capacity ($f_o)$, specified in AS 3700–2011, to determine the design compressive stress: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{f}_o=\mathit{\varphi }{f}_{\mathrm{c}}{A}_{\mathrm{b}}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{d}}=k{f}_o,\end{array}$$ where $f_{\mathrm{c}}$ is the characteristic compressive strength of the masonry, $\mathit{\varphi }$ is the capacity reduction factor, $A_{\mathrm{b}}$ is the bedded area of the masonry (brick width $\times$ length) and $k$ is a reduction factor based on the wall design. The characteristic compressive strength is determined using the unconfined compressive strength tests of Martelli (2011) on building materials sourced from Arequipa. Presuming the mortar is of relatively low quality (M2), the characteristic compressive strengths (according to AS 3700–2011) are 3.8 MPa for ignimbrite masonry and between 3.5 and 4.54 MPa for terracotta masonry. The slenderness reduction factor, $k$, describes the susceptibility to buckling. Following AS 3700–2011, this factor is calculated as $${\displaystyle }k=\mathrm{0.67}-\mathrm{0.02}\left(S_{\mathrm{rs}},-,\mathrm{14}\right)$$ for buildings with a reinforced concrete roof or floor (i.e. typologies 3–6C), and $${\displaystyle }k=\mathrm{0.67}-\mathrm{0.025}(S_{\mathrm{rs}}-\mathrm{10})$$ for buildings with other roof or floor supports (typologies 1–2). This factor requires a calculation of the slenderness ratio, $S_{\mathrm{rs}}$: $${\displaystyle }{S}_{\mathrm{rs}}={\displaystyle \frac{a_{\mathrm{v}}H}{k_{\mathrm{t}}b}},$$ where $H$ is the height between floors or supports, taken as 2.8 m for reinforced concrete type buildings and 3 m for non-reinforced buildings (Martelli, 2011). The vertical slenderness coefficient, $a_{\mathrm{v}}$, is determined from the lateral support along the top edge of the wall. Walls with roof support (types 3–6) have a coefficient of 1, while unsupported walls (types 1A–2B) act as a cantilever and have a coefficient of 2.5. Considering the thinnest bricks, the slenderness coefficient is negative for building types 1A–2B as the design is out of the range of those considered in AS 3700–2011. Acknowledging the low strength of these frequently makeshift structures, the slenderness coefficient is therefore set to 0.01. The thickness coefficient, $k_{\mathrm{t}}$, takes into account the strength of supporting columns. This coefficient is set to 1 for non-reinforced frame buildings and is dependent on the spacing and thickness of reinforced beams within the masonry for reinforced buildings. Estimates of the spacing and thickness suggest that the coefficient will be between 1.4 and 2 for type 6A–6C buildings and between 1 and 1.2 for type 4 buildings. The large spacing between reinforced columns and their relative width, pictured in Thouret et al. (2014), is responsible for the much lower coefficients assigned to type 4 buildings.

The Supplement related to this article is available online at doi:10.5194/nhess-17-703-2017-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

The authors would like to thank the Civil Defence office in Arequipa (Instituto Nacional de Defensa Civil INDECI), in particular the Regional office (A. Arguedas) and the Provincial/City office (J. Vasquez) for support during field work as well as the students of the department of geology of the University Nacional San Aguistin in Arequipa. J.-C. Thouret's work in Arequipa has been supported by the Labex CLERVOLC (contribution number 235), the PICS CNRS programme and the French Embassy in Lima. Edited by: G. Macedonio Reviewed by: two anonymous referees