The use of FLO 2 D numerical code in lahar hazard evaluation at Popocatépetl volcano : a 2001 lahar scenario

Lahar modeling represents an excellent tool for designing hazard maps. It allows the definition of potential inundation zones for different lahar magnitude scenarios and sediment concentrations. Here, we present the results obtained for the 2001 syneruptive lahar at Popocatépetl volcano, based on simulations performed with FLO2D software. An accurate delineation of this event is needed, since it is one of the possible scenarios considered if magmatic activity increases its magnitude. One of the main issues for lahar simulation using FLO2D is the calibration of the input hydrograph and rheological flow properties. Here, we verified that geophone data can be properly calibrated by means of peak discharge calculations obtained by the superelevation method. Digital elevation model resolution also resulted as an important factor in defining the reliability of the simulated flows. Simulation results clearly show the influence of sediment concentrations and rheological properties on lahar depth and distribution. Modifying rheological properties during lahar simulation strongly affects lahar distribution. More viscous lahars have a more restricted aerial distribution and thicker depths, and resulting velocities are noticeably smaller. FLO2D proved to be a very successful tool for delimitating lahar inundation zones as well as generating different lahar scenarios not only related to lahar volume or magnitude, but also taking into account different sediment concentrations and rheologies widely documented as influencing laharprone areas.


Introduction
Lahar phenomena represent one of the major threats at volcanoes.They are mixtures of water and sediments that flow on volcano slopes (Smith and Fritz, 1989), remobilizing volcanic products.Their triggering mechanisms include glacial melting during volcanic activity, intense rainfall or dam break.These triggering mechanisms allowed them to occur without being related to volcanic eruptions, which makes them very hazardous.Lahars have caused several catastrophes worldwide vastly reported in the literature (Voight, 1990;Siebe et al., 1996;Lavigne et al., 2000;Scott et al., 2005).Additionally, small eruptions can trigger very large events (e.g., the 1985 Nevado del Ruiz eruption; Williams, 1987).Based on these facts, several tools need to be used in order to delimitate lahar inundation zones precisely around volcanoes and to minimize our vulnerability to these phenomena.
Areas prone to lahar inundation can be adequately delineated by different methods.Field data are useful in our understanding of the magnitude and distal reach of past events in a specific volcano, but they have the disadvantage that, in very active gullies, lahar deposits are easily eroded by subsequent events (i.e., Graettinger et al., 2010).Misinterpretation of lahar deposits, including their magnitude and recurrence, could lead to exclusion of some areas that could be affected during major events.In recent years, one of the major instruments used to replicate old events and accurately delineate lahar active paths is lahar modeling (i.e., Iverson et al., 1998;Davila et al., 2007;Muñoz-Salinas et al., 2009a;Williams et al., 2008;Worni et al., 2012).Different simulations codes include LaharZ (Schilling, 1998), Titan2D two-phase flow (Williams et al., 2008), andFLO2D (O'Brien et al., 1993).
Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Caballero and L. Capra: The use of FLO2D in lahar hazard evaluation at Popocatépetl volcano All of them have been documented as reproducing, with good to high degrees of accuracy, lahar events.
In this work, the 2001 lahar that occurred at Popocatépetl volcano just after the emplacement of a pumice flow from an 8 km eruptive column is modeled using FLO2D numerical code (O'Brien et al., 1993), a tool recently successfully used to simulate debris flows in volcanic environments (Worni et al., 2012).FLO2D allows reproduction of debris and hyperconcentrated flows and offers the advantage of modifying flow sediment concentration and debris flow rheology, factors that influence lahar distribution and depth.Geophone data are used here to reconstruct initial hydrograph and sediment concentrations.Results obtained were compared with field data and demonstrated a good agreement in thickness and flow distribution.A comparison with previously published data related to lahar velocity (Muñoz-Salinas et al., 2007) is also made.Additionally, lahars with fluctuating sediment concentrations but with similar volumes are simulated to observe the influence of sediment concentration and rheological behavior on lahar distribution and its hazard.Finally, a sensitivity analysis of digital elevation model (DEM) resolution is also performed here to better constrain all the factors influencing numerical simulation with FLO2D code.

Background
Popocatépetl volcano is one of the most active volcanoes in Mexico, whose activity has increased since 1994, representing a threat to more than a million people in its surroundings.It belongs to Sierra Nevada, a N-S 20-volcano chain, located in the central sector of the Trans-Mexican Volcanic Belt (Fig. 1).Lahars at Popocatépetl volcano have been triggered mainly due to volcanic activity (Siebe et al., 1996;González, 2000;Capra et al., 2004;Muñoz-Salinas et al., 2010).The magnitudes of these events in Popocatépetl's history have varied from small lahars generated by intense rainfall (Capra et al., 2004) to very huge events that travelled more than 50 km away from the volcano associated with Plinian activity (González, 2000;Siebe et al., 1996Siebe et al., , 1999)).Large lahars triggered by major eruptions affected mexican pre-Hispanic civilizations, and their deposits are distributed where major cities, like Puebla and Cuatla, are now settled (Siebe et al., 1996).
On 7 May 2013, Popocatépetl volcano incremented the intensity of its volcanic activity with intra-crater lava dome extrusions rapidly destroyed by moderated explosions that formed columns of up to 3-4 km in height.This behavior put the scientific community on alert, causing the raising of volcanic alert to yellow phase three by 12 May 2013 (CE-NAPRED, internal reports).One of the possible scenarios considered during this crisis was an event similar to the 2001 explosive activity, which was characterized by an 8 km eruptive column and the subsequent formation of pumice flows up to 6 km from the crater (Fig. 2; Sheridan et al., 2001).Partial melting of the glacier remobilized the new deposits, forming a lahar a few hours after on the northeastern flank of the volcano, along Huiloac Gorge, almost reaching Santiago Xalitzintla (Capra et al., 2004).A similar event occurred in 1997, when intense ash fall on top of the glacier promoted a rapid water release that resulted in a lahar that also inundated Santiago Xalitzintla (with more than 2196 inhabitants; INEGI, 2010) along the Huiloac-Tenenepanco ravine (Capra et al., 2004).Nowadays, a glacier is no longer present on the Popocatépetl cone (Julio- Miranda et al., 2008) but, as at other Mexican active volcanoes such as the Volcán de Colima, lahars are very frequent during the rainy season (Davila et al., 2007;Capra et al., 2010), especially when pyroclastic material from an ongoing activity accumulates on the upper slope of the volcano.The 2001 explosive episode probably represents the biggest eruption since 1994, when the volcano reactivated after more than 70 years of inactivity (since 1927) (Martin-Del Pozzo, 2012).The occurrence of a similar scenario makes the reproduction of this event of fundamental importance, in order to delimitate lahar hazard zones accurately.

2001 lahar event
On 22 January 2001, an explosive episode (Fig. 2), characterized by an 8 km eruptive column, gave place to pumice flows that emplaced mostly in the northern sector of the cone up to 4-6 km in distances from the crater (Espinasa Pereña, 2012).This episode of activity had an estimated VEI of 3-4 (Espinasa Pereña, 2012), and has been one of the most violent phases since 1994.
As the flows descended over the glacier, they eroded and melted part of it, and after approximately 5 hours, a lahar descended along the Huiloac Gorge for 15 km, almost reaching Santiago Xalitzintla (Capra et al., 2004;Tanarro et al., 2010).Based on geophone data, the event started at 22:32 and ended at 23:20 GMT.The load of the produced lahar was calculated as 1.6 × 10 5 m 3 and reached speeds from 1.3 to 13.8 m s −1 (Muñoz-Salinas et al., 2007, 2009b).
The deposit has been described as a massive debris flow, with an average thickness of 70 cm up to a maximum of 150 cm in the middle reach.It consists of a matrix-supported unit, with rounded to sub-rounded clasts embedded in a matrix that consist of up to 70 % of sand and silt, and up to 1 % of clay (Fig. 3).Components consist of pumices and andesitic fragments (Capra et al., 2004).
Six lahars have occurred after the 2001 lahar events; four of them occurred during the dry season, and were caused by snowmelt.The others were generated in the middle to final periods of the rainfall season, and were associated with torrential rains (Muñoz-Salinas et al., 2010).Based on this, lahar hazard evaluation should take into account different triggering scenarios and the sediment available for remobilization inside volcano gorges.Special attention should be paid along Huilac Gorge, the most active basin in regards to lahars.

Methodology
The FLO2D numerical model (O'Brien et al., 1993)  It ranges between 6 up to 17 wt %.This fact makes viscous forces very important in comparison to dispersive forces.Based on that, FLO2D is an appropriate tool for simulating this event.
FLO2D uses an established grid over the topographic area (computational domain) and calculates flow depth and velocity between the grid elements.Flow velocity in eight di-rections is calculated independently, with an algorithm that includes flow geometry, roughness (Manning n value), the slope between two grid elements, and the wetted perimeter.Once the velocity is solved, it is multiplied by the crosssectional area, and discharge is calculated.Flood-wave advance and direction are determined by topography and resistance to flow.Debris and hyperconcentrated flows can be simulated by a quadratic rheological model that involves sediment concentration, yield stress and viscosity terms.FLO2D predicts viscous fluid motion as a function of sediment concentration.
Input data required for mud/debris flow simulations comprise inflow (hydrograph) and outflow cells, Manning n values, rheological parameters (α and β), laminar flow resistance (K), and a limiting Froude number.Rheological parameters (α and β) are empirical parameters that relate yield stress and viscosity to sediment concentration by the empirical relationships (O'Brien and Julien, 1988) In a FLO2D mud/debris flow simulation, the water and sediment volumes as well as sediment concentrations for every grid element and time step are computed.
Analysis of the different input parameters was performed by Boniello et al. (2010), Hsu et al. (2010), and Worni et al. (2012).They include data from the FLO2D user's manual, as well as back analysis to fit field and debris flow data from the literature in order to obtain confident rheological coefficients.All of their results point to DEM resolution, rheological parameters and Manning coefficients like the most sensitive parameters for a successful simulation in FLO2D.
In this work, debris flow simulation was performed using a 10 m digital elevation model (DEM).The 10 m DEM was created from contour lines at 10 m provided by INEGI (Instituto Nacional de Estadística y Geografía).The grid imposed over topography for input data and FLO2D depth and velocity calculations was also at 10 m.A 10 m DEM and grid allowed flow simulation results to be obtained on the same scale.

Inflow data
For hydrograph reconstruction, data recorded from geophone stations installed by CENAPRED and USGS were used (Fig. 4).Data from geophone PFM3 were chosen because of their proximity to the point where all tributary rivers converge to Huiloac Gorge.Besides, Huiloac Gorge has no tributary gullies; hence, its budget is determined by water availability behind this point (Tanarro et al., 2010).
Data from the full-frequency band (10-250 Hz) were correlated with peak discharge calculations from Muñoz-Salinas et al. (2007) obtained by the superelevation method.To avoid surging during flow simulations, extreme peaks observed in raw data were suppressed, but the overall shape of the hydrograph was maintained.The conversion from AFM units (cm s −1 × 10 −6 ) to discharge units (m 3 s −1 ) was made assuming that the maximum peak discharge (Muñoz-Salinas et al., 2007) corresponded to the maximum value of the hydrograph (Suwa et al., 2000).Flow concentration was established assuming a maximum flow concentration of 0.5 during the peak flow and the main body of the flow, which is in agreement with field observation and textural characteristics of the deposit (Capra et al., 2004).The recessional part of the flow was assumed to be a hyperconcentrated flow (Fig. 4).

Manning n values
FLO2D requires Manning n values to calculate flow resistance of the turbulent and dispersive shear stress by the following equation: where n td is the flow resistance of the turbulent and dispersive shear stresses, n is the Manning n value, and C v is the sediment concentration.
Manning values were established using USGS empirical tables.Base values of Manning's n for natural channels composed of cobbles and boulders range from 0.030 to 0.070 (Phillips and Tadayon, 2006), so an intermediate value of 0.05 was chosen.Huiloac Gorge has different properties related to channel irregularities, channel alignment, obstructions, and vegetation, so some adjustments to base n values  were needed in order to characterize channel roughness better.Three different Manning n values were used, based on the channel characteristics mentioned above (Fig. 5).A value of 0.064 was used for the segment were the channel is very narrow, the base is composed of cobble and boulders, and a moderate degree of irregularities is observed.Variations in channel cross sections alternates occasionally, and the effects of obstructions and vegetation are negligible.The floodplain on both sides of Huiloac Gorge has severe vegetation obstructions, giving a Manning value of 0.118, which is also in accordance with Chow (1959) values for floodplain analysis.Last adjustment in n values was done where Huiloac Gorge opens to different branches, where channel segments are more sinuous, shallower and wider, with medium vegetation degree giving a value of 0.081.

Rheological coefficients α and β
Rheological coefficients α and β were chosen based on the O' Brien and Julien (1988) data.Four different values of α and β for the 2001 lahar were simulated in order to observe the influence of rheological characteristics on lahar distribution (Table 1).The first scenario was based on the repro-duction of the 2001 lahar.Values α and β were chosen using the sample, with the more similar granulometric composition of this lahar corresponding to the Glenwood 2 sample, also in accordance with values used in previous works (Worni et al., 2012).The second scenario was a more diluted lahar; α and β values were maintained, and only a hydrograph with minor sediment concentration was used.Third and fourth scenarios were done for more viscous lahars, and α and β values correspond to samples with more clay content in order to raise the yield strength and viscosity of the flow.Clay contents of 4.8 (Glenwood 3 sample) and 6.8 (Glenwood 1 sample) were used for simulations 3 and 4, respectively.

Laminar flow resistance (K) and limiting Froude number (F )
The laminar flow resistance value was set in 2000 according to the values suggested in the user's guide based on substrate characteristics.A limiting Froude number of 0.9 was set to maintain a subcritical flow regime.Field data point to this behavior, since flows in the supercritical regime develop sedimentary structures like cross-stratification and dunes.Be-sides, supercritical flows are suppressed by high rates of sediment transport or high fine sediment content like in the 2001 lahar (up to 16 wt % of silt and clay).Based on the previous statements, a limiting Froude number of 0.9 reflects an accurate flow regime of the real lahar.

2001 event
Results of the 2001 lahar simulation obtained by FLO2D are shown in Fig. 6.Four aspects were compared between simulation and field and indirectly estimated data (Muñoz-Salinas et al., 2007) to observe simulation accuracy: spatial distribution, flow depth, flow velocity, and flow volume.
Flow simulation reaches approximately 9.5 km from the hydrograph point, and stops almost 1.7 km before reaching Santiago Xalitzintla.Although the last outcrops of the 2001 lahar deposit are found approximately 3.5 km before Santiago Xalitzintla (Fig. 1), distal reaches obtained by FLO2D simulation would represent the watery recessional parts of the lahar.Flow depth comparison was made using the river transversal profiles and stratigraphic columns described by Capra et al. (2004).River profiles (Fig. 6) show flow depths obtained during simulations and deposit thickness at the same point.Flow depths consider sediment plus water, so the inundation line must always be above the deposit.At proximal reach, the flow is 4.2 m, while deposit thickness is reported to be from 0.5 to 0.8 m.These variations are because the simulated maximum flow depths are located in the central channel, while deposit thickness was measured near the walls where deposition takes place during lahar emplacement.In fact, direct observation of lahars clearly shows that rapid and frequent flow discharge fluctuations occur as the deposit is progressively piled up on the channel sides (Vázquez et al., 2014), behavior that can clearly explain differences in flow depth and deposit thickness.At medial reach, simulation gives a maximum flow depth of 3.2 m, and deposit thickness varies from 0.15 to 1.90 m.The distal part of the simulation has a maximum flow depth of 1.6 m, and deposit thickness varies widely, from 0.65 to 0.12 m.
Lahar velocity distributions observed during flow simulations were compared with velocity data reported by Muñoz-Salinas et al. (2007) by the superelevation method (Fig. 7).Maximum lahar velocities calculated by FLO2D are circa 9.5 m s −1 , and are achieved in proximal facies.Muñoz-Salinas et al. (2007) estimated velocities between 13.8 and 1.5 m s −1 for the 2001 lahar.Considering the wide range of estimated velocities and the error estimated for the superelevation method of 15 % (Pierson, 1985), velocities estimated during the lahar simulation are a good approximation of estimated flow velocities.
Finally, the flow volume calculated by FLO2D is circa 3.6 × 10 5 m 3 (water bulked with sediments).The mean amount of water is 2.1 × 10 5 m 3 , and the sediment volume is 1.5 × 10 5 m 3 , considering the longitudinal profile in sediment concentration from the hydrograph input data (Fig. 4).Capra et al. (2004), Julio Miranda et al. (2005), and Muñoz-Salinas et al. (2009b) estimated a lahar volume within 1.6-2.4× 10 5 m 3 that is on the same order of magnitude as that calculated by FLO2D.
Based on previous statements, results of the FLO2D 2001 lahar simulation are in good agreement with previous works.

Different lahar scenarios at Popocatéptl volcano
Flow simulations allow the generation of different lahar scenarios.These scenarios should be created based on the magnitude of past events and the present conditions in lahar-active paths.Most lahar simulations are magnitudebased scenarios or reproduction of some events (Aguilera et al., 2004;Williams et al., 2008;Worni et al., 2012).Strong variations in lahar sediment concentration or transformation from debris flow to hyperconcentrated flow or vice versa can occur (Pierson and Scott, 1985).Based on this, we believe that different scenarios should consider different sediment concentrations for lahars of the same magnitude, since this characteristic will influence travel distance and therefore will strongly influence hazard maps.
Since its reactivation in 1994, the largest lahars generated at Popocatépetl volcano along Huiloac Gorge were triggered by peaks in volcanic activity that melted part of the glacier.The most hazardous scenarios are similar to the 1997 and 2001 lahars, with flow volumes estimated between 2 and 3.3 × 10 5 m 3 (Capra et al., 2004;Julio Miranda et al., 2005;Muñoz-Salinas et al., 2009b).In addition, based on lahars that occurred at Popocatepetl volcano in 1997and 2001, Capra et al. (2004) point to the evidence that, depending on the nature of the material contributing to the initiation of a lahar, the flow behavior can change greatly, mostly depending on fine (silt and clay) content.The 1997 lahar was fine depleted and rapidly diluted to a hyperconcentrated flow at the same distance, whereas the 2001 lahar still maintained a high sediment concentration in the range of a debris flow.Therefore, fine sediment content influences dynamic behavior, transformations and dynamic rheological behavior.
FLO2D allows modification of concentration and rheological properties of the lahars modifying initial input hydrograph and α and β coefficients.Based on that, α and β coefficients and sediment concentrations were adapted to consider these possible different behaviors observed in the 1997 and 2001 lahars (Table 1).Three alternative scenarios (AS) were proposed for events of the same magnitude as the 2011 lahar.AS1 has a lower sediment concentration but the same rheological coefficients of the 2001 lahar.AS2 has rheological properties equivalent to a fine content of 4.8 %.Finally,  AS3 has rheological properties equivalent to a fine content of 6.8 %.The sediment concentrations of AS2 and AS3 were the same as the 2001 lahar.
Simulation results, which include inundation areas, flow depths and velocities, are presented in Figs. 8 and 9.The influence of concentration and rheological properties on lahar distribution, depth and velocity is clear to see.It is worth mentioning that no variation in flow width is seen to be related to rheological properties of the lahar.AS1, the more diluted lahar (Figs.8a and 9a), gives a spatial distribution, flow depths, and velocities very similar to that of 2001.Based on that, hydrograph sediment concentration did not deeply affect lahar distribution and flow properties.On the other hand, AS2 and AS3, where rheological parameters are modified, display different flow behaviors and distribution.AS2, a more viscous lahar, has a more restricted aerial distribution, and it stops approximately 3.2 km before reaching Santiago Xalitzintla (Figs. 8b and 9a).In the proximal parts, it achieves almost 5.2 m in flow depth, with a clear steepest front.There is no recessional watery part of the flow like in the more diluted 2001 lahar.Resulting velocities are smaller, reaching a maximum value of 8.4 m s −1 .AS3, the more viscous lahar scenario (Table 1, Figs. 8c and 9a), reaches a maximum flow depth of 6.2 m, almost 2 m thicker than that of the  2001 lahar.It stops at 3.5 km before arriving at Santiago Xalitzintla, and is the slowest lahar, with a maximum velocity of 6.8 m s −1 .
A comparison of the 2001 lahar simulation and the alternative scenarios shows how important it is to consider sediment concentration and rheological properties of lahars for hazard evaluation (Fig. 8).Differences in distal reaches are almost 2 km between more diluted and more viscous lahars (Fig. 9a).A point analysis shows that rheological properties can cause strong differences in flow depth.Differences in flow depths  are up to 2.5 m, especially in the proximal zone and in the central portion of the channel (Fig. 9b).
Based on the above, it is important to consider these strong differences when lahar hazard maps are made.A final map integrating the results of the 2001 lahar and the alternative scenarios simulated is shown in Fig. 10.The map uses the maximum flow depths at each cell based on the four different scenarios simulated here and the maximum aerial distribution.

FLO2D sensitivity analysis to DEM resolution
The influence on DEM resolution has been widely documented for granular flows (Stevens et al., 2002;Huggel et al., 2008;Capra et al., 2011).To analyze FLO2D sensitivity to DEM resolution, flow simulations were performed by resampling the 10 m DEM to a 30 m DEM, and the two extreme lahar scenarios (2001 lahar and AS3) were used.The grid imposed over DEM for computation was also 30 m in resolution.Input parameters (inflow hydrograph, Manning n value, and Froude number) were kept constant in these simulations.The results are displayed in Fig. 11.Lahar distributions obtained from the 30 m DEM resolution as well as flow depths are lower for both scenarios (Fig. 11).Lahar distribution is reduced by almost 15 %, distal reaches diminish for up to 1.5 km, and maximum flow depths are reduced by 2-2.5 m, so a low-resolution DEM led to underestimation of maximum reaches.This effect was also observed using the laharZ code (Schilling, 1998), since lowresolution DEM smooths valley channels, promoting lateral inundation in spite of longer runout (Davila et al., 2007).

Conclusions
Lahar modeling represents an excellent tool for designing hazard maps.It allows, if they are properly calibrated, the definition of potential inundation zones for different lahar magnitude scenarios and sediment concentrations.
Results presented here proved FLO2D to be a useful tool in lahar modeling at Popocatépetl volcano.Very good agreement between field data (Capra et al., 2004), flow behavior (Muñoz-Salinas et al., 2007) and the 2001 lahar simulation was observed.One of the main issues for lahar simulation using FLO2D is the calibration of the input hydrograph.Here, we verified that geophone data could be properly calibrated by means of peak discharge calculations obtained by the superelevation method.
Lahars in active volcanoes represent a major threat.They can be triggered by volcanic activity, i.e., glacier melting, by intense rainfall, by a crater dam break or by the transformation of debris avalanches into these phenomena.The consequences of such diverse scenarios are lahars with important variations in sediment concentration and behaviors, varying between debris and hyperconcentrated flows.Both types of flows have different dynamic and rheological behaviors, resulting in different flow velocities and maximum runouts; hence, they represent different threats.One of the shortcomings of lahar modeling is the creation of different scenarios based only on lahar magnitude.Results of lahar simulation presented here for the 2001 lahar involve lahars of the same magnitude, but with different sediment concentrations and fine contents.Differences in travel distances are up to 1.8 km for events of the same magnitude, and maximum flow depths can have differences of more than 2.5 m.Besides, modifications in lahar rheology not only affect its distribution, but they also modify flow thicknesses and velocities.More viscous lahars have more restricted distributions, but are deeper than more fluid ones, as is observed by different scenario simulations.
Results of lahar modeling can be made with an uncertainty of 2-4 km (Huggel et al., 2008).This effect could be enhanced by not taking into account lahar rheology.Modeling lahar events based only on volumes can enhance the uncertainty in the results, but this can be avoided by recreating different rheological scenarios.DEM resolution as well as its quality is also a key point in flow simulation, as also observed for other numerical models for granular flow (i.e., Stevens et al., 2002;Davila et al., 2007;Capra et al., 2010;Worni et al., 2012).
FLO2D proved to be a very successful tool for delineating lahar inundation zones as well as generating different lahar scenarios not only related to lahar volume or magnitude, but also taking into account different sediment concentrations and rheologies widely documented as influencing lahar-prone areas.

Figure 1 .
Figure 1.(a) Localization of Popocatépetl volcano in the Trans-Mexican Volcanic Belt.(b) Popocatépetl volcano.The orange line represents the 2001 lahar deposit distribution.Red circles are CENAPRED geophone locations.Green circles indicate the location of stratigraphic sections collected by Capra et al. (2004).
routes flood hydrographs over unconfined channels solving the full dynamic wave equation.FLO2D uses a quadratic shear stress model that involves five stress terms, including the cohesive yield stress, Mohr-Coulomb shear, viscous shear stress, the turbulent shear stress, and the dispersive shear stress.The usage of such terms describes the continuum of flow regimes from viscous to turbulent/dispersive flow.The incorporation of such terms allows a good reproduction of lahars where interaction of non-cohesive sediments develops dispersive stresses to the flow.It also favors the simulation of lahars with a high proportion of fine sediments where viscous forces dominate flow behavior.The 2001 lahar of Popocatépetl volcano has a high proportion of fine sediments (silt and clay).

Figure 3 .
Figure 3. Photographs of the 2001 lahar deposit.(a) Proximal facies of the deposit.(b) Textures of the 2001 lahar deposit where the massive structure and clasts dispersed in a sand-silty matrix are observable.

Figure 4 .
Figure 4. Geophone data used to reconstruct the hydrograph of the 2001 lahar.Q max used to calibrate geophone data was obtained by Muñoz-Salinas et al. (2007).Left scale refers to the flow discharge (m 3 s −1 ), and the right scale AFM units.In the upper part, sediment concentration for each portion of the hydrograph is indicated.Geophone data are from CENAPRED.

Figure 5 .
Figure 5. Distribution of Manning coefficients used during FLO2D simulations.

Figure 6 .
Figure 6.Simulation results of the 2001 lahar for flow depth.Green circles correspond to field data published by Capra et al. (2004).Red circle is geophone PFM3 from CENAPRED.Flow sections compare flow depths calculated by FLO2D with field data collected by Capra et al. (2004).Shaded areas represent the 2001 lahar deposit.Unit A refers to the 1997 lahar deposit and unit B to the 2001 lahar deposit.

Figure 7 .
Figure 7. Velocity distribution obtained during flow simulation (m s −1 ).Yellow triangles refer to the locations of flow velocity calculated by Muñoz-Salinas et al. (2007) used here for comparison.

Figure 8 .
Figure 8. Results of alternative scenarios for lahars with an equivalent volume to the 2001 lahar, but with different sediment concentrations or rheological properties.(a) AS1, more diluted lahar.(b) AS2, lahar with 4.8 wt % of clay.(c) AS3, lahar with 6.8 wt % of clay.AS2 and AS3 represent high yield strength and more viscous lahars.The first column represents flow depth, the second column, flow velocity.Green circles correspond to field data published by Capra et al. (2004).Red circles are geophone PFM3 from CENAPRED.

Figure 9 .
Figure 9.Comparison of the alternative scenarios with the 2001 lahar simulation.(a) Area of inundation.(b) Difference in lahar flow depth (m) between the different alternative scenarios.

Figure 10 .
Figure 10.Final map constructed based on maximum flow depths (m) from the different lahar scenarios simulated here.

Figure 11 .
Figure 11.Lahar simulation with the 30 m DEM.(a) 2001 lahar simulation.(b) Simulation of AS3 (more viscous lahar).Green (a) and red (b) lines indicate simulation results with the 10 m DEM.

Table 1 .
Brien and Julien (1988)oefficients used for the different flow scenarios.Yield strength and viscosity were obtained by the empirical formulas used by FLO2D and from O'Brien and Julien (1988).