Assessing lahars from ice-capped volcanoes using ASTER satellite data, the SRTM DTM and two different flow models: case study on

Abstract. Lahars frequently affect the slopes of ice-capped volcanoes. They can be triggered by volcano-ice interactions during eruptions but also by processes such as intense precipitation or by outbursts of glacial water bodies not directly related to eruptive activity. We use remote sensing, GIS and lahar models in combination with ground observations for an initial lahar hazard assessment on Iztaccihuatl volcano (5230 m a.s.l.), considering also possible future developments of the glaciers on the volcano. Observations of the glacial extent are important for estimations of future hazard scenarios, especially in a rapidly changing tropical glacial environment. In this study, analysis of the glaciers on Iztaccihuatl shows a dramatic retreat during the last 150 years: the glaciated area in 2007 corresponds to only 4% of the one in 1850 AD and the glaciers are expected to survive no later than the year 2020. Most of the glacial retreat is considered to be related to climate change but in-situ observations suggest also that geo- and hydrothermal heat flow at the summit-crater area can not be ruled out, as emphasized by fumarolic activity documented in a former study. However, development of crater lakes and englacial water reservoirs are supposed to be a more realistic scenario for lahar generation than sudden ice melting by rigorous volcano-ice interaction. Model calculations show that possible outburst floods have to be larger than ~5×105 m3 or to achieve an H/L ratio (Height/runout Length) of 0.2 and lower in order to reach the populated lower flanks. This threshold volume equals 2.4% melted ice of Iztaccihuatl's total ice volume in 2007, assuming 40% water and 60% volumetric debris content of a potential lahar. The model sensitivity analysis reveals important effects of the generic type of the Digital Terrain Model (DTM) used on the results. As a consequence, the predicted affected areas can vary significantly. For such hazard zonation, we therefore suggest the use of different types of DTMs and flow models, followed by a careful comparison and interpretation of the results.


Introduction
Lahars are mixtures of water and volcanic debris of varying size and solids fraction, which frequently affect the steep, unconsolidated slopes of volcanoes.Lahars represent one of the most far-reaching threats in volcanic terrains.They may occur not only prior to or during eruptions, but also years later, and can be triggered by torrential storms, lake outbursts and earthquakes (Verstappen, 1992;Pierson, 1998;Vallance, 2000).Stream flows are muddy water streams with a sediment concentration below 20% in volume.Lahars however, can be further separated according to their sediment concentration: into hyperconcentrated flows between 20 and 50-60% in volume, and into debris flows if the volumetric sediment concentration is higher than 50-60% (Pierson and Scott, 1985;Pierson, 1986;Smith and Lowe, 1991;Coussot and Meunier, 1996;Vallance, 2000;Lavigne and Thouret, 2002).The lahars simulated in this study represent debris flow types with sediment concentrations similar to the observed ice-melt-triggered lahars on Popocatépetl in 1997and 2001(Julio Miranda and Delgado Granados, 2003;Capra et al., 2004;Julio Miranda et al., 2005).
Satellite remote sensing has been widely used in volcanic studies (Kerle and Oppenheimer, 2002;Pieri and Abrams, 2004) while GIS tools and computer based models became increasingly important during the past years, for example, for calculating the runout paths of potential lahars (Iverson et Published by Copernicus Publications on behalf of the European Geosciences Union.City, the study area and Puebla is derived from SRTM3 elevation model.The dotted area shown as 'prehispanic lahars' represent the approximate extent of areas flooded by repeated lahar events from Popocatépetl, Iztaccíhuatl and La Malinche after Siebe et al. (1996).al., 1998;Julio Miranda and Delgado Granados, 2003;Sheridan et al., 2005;Walder et al., 2005;Hubbard et al., 2006;Davila et al., 2007;Huggel et al., 2008).By applying such tools to Iztaccíhuatl volcano in Mexico, we show how an initial hazard evaluation can be carried out relatively fast and at low cost.

Iztaccíhuatl volcano: setting and background
Iztaccíhuatl (5230 m a.s.l.) is a major Quaternary volcanic complex within the Trans-Mexican Volcanic Belt (TMVB) consisting of various craters, with the principal volcanic vents aligned along a N-NW to S-SE axis (Figs. 1 and  2).The volcanic activity has been variously characterized as essentially extinct (Siebe and Macías, 2006), dormant (Nixon, 1989) or active (Delgado Granados et al., 2005).According to Nixon (1989), the last activity with dacite flows took place at the northern flank at approximately 80 000 years BP.A major structural failure producing a giant debris avalanche to the E-SE side occurred several thousand years ago, just south of La Panza at the site of Las Rodillas (Siebe et al., 1995).Delgado Granados et al. (2005) reported diffuse gas emissions at La Panza of Iztaccíhuatl.Based on this and the fact that lava flows from Iztaccíhuatl overlie volcanic deposits from Popocatépetl which are younger than 5000 years (Siebe et al., 1995), these authors considered Iztaccíhuatl to be an active volcano, following the definition that volcanoes with eruptive activity within the last 10 000 years are active.
After Pico de Orizaba (also known as Citlaltépetl, 5675 m a.s.l.) and Popocatépetl (5452 m a.s.l.), Iztaccíhuatl is the third highest peak in Mexico and, like the other two volcanoes, is still covered by ice (Delgado Granados, 2007).The presence of glaciers in Mexico is noteworthy because at this latitudinal range there are no other glaciers worldwide.In reference to the silhouette of the volcanic edifice (similar to the shape of a woman lying on her back) the two main glacier systems are called El Pecho (the chest) at the summit and La Panza (the belly) 1 km to S-SE (Fig. 2 and Fig. 4).All the glaciers on Iztaccíhuatl have shown extensive retreat during the past 150 years.
Hazard evaluation for ice-volcano interactions is crucial at Mexico's active ice-covered volcanoes.Since 1994, when the recent eruptive activity of Popocatépetl volcano began, lahars related to ice melting developed in 1997 and 2001.Both events reached the nearby village of Santiago Xalitzintla, which is located 14 km away from the summit (Julio Miranda and Delgado Granados, 2003;Capra et al., 2004;Julio Miranda et al., 2005;Andrés et al., 2007).
In prehistoric time, the area around Popocatépetl and Iztaccíhuatl volcanoes was extensively inundated by lahars which destroyed various human settlements (Fig. 1; Siebe et al., 1996).The material (essentially ash and pumice from Popocatépetl) was released from the flanks of Popocatépetl, Iztaccíhuatl and La Malinche volcano, and flooded over 2600 km 2 (GIS calculation from the map in Siebe et al., 1996).If we assume an average thickness of the deposits of 2-5 m and a volumetric sediment content of 60%, the total lahar volume including water is about 9-22 km 3 during numerous events over many years.These lahars are assumed to have been mainly related to heavy precipitation, but also to volcano-ice interactions (see Sect. 1.3).

Predisposition for lahars on Iztaccíhuatl
For the generation of lahars, the following three conditions are important: (1) availability of potentially erodible soil and debris (including thickness of source deposits and their physical characteristics), (2) critical slope and channel gradient, (3) availability and amount of water.
Loose rock material is abundant at most of the summit areas of Iztaccíhuatl, consisting mainly of eroded dacites that are accumulated in large rock talus and moraines (Nixon, 1989).Pyroclastic material is observed on several parts of the volcano.At the lower and southernmost part (<4000 m a.s.l.), tephras from Popocatépetl volcano deposited over the last 5000 years can be found.Above 4300 m a.s.l., old and loose tephras from ancient eruptive activity from Iztaccíhuatl are abundant.Vegetation is lacking, so that the ground is not stabilized by plant matter.Permafrost can be found at northern expositions above 4900 m a.s.l.(Palacios et al., 2007).
A digital terrain model (DTM) derived from the shuttle radar topography mission (SRTM) data of February 2000, reveals inclinations of >35-45 • (70-100%) around the summit area.Although the angle of repose of loose tephra has not been evaluated, we assume conditions close to the stability limits.Given that criteria (1) and ( 2) are fulfilled at the summit area of Iztaccíhuatl, the availability of water has to be the deciding factor for lahar formation.On Iztaccíhuatl volcano it can be released by precipitation or melting of the available ice and snow deposits.

Rain-triggered and glacier-related lahars
Frequency of lahars and precipitation intensity are often strongly correlated (Lavigne and Thouret, 2002;Van Westen and Daag, 2005).On the slopes of Iztaccíhuatl and Popocatépetl, small lahars occur commonly during the rainy  season.For instance, a rain-triggered lahar occurred on 20 June 1994 (six months prior to the initiation of the current eruptive activity at Popocatépetl volcano).After a heavy rainfall, a relatively small lahar reached, entered, and deposited volcanic ash and debris on the streets of the village Santiago Xalitzintla, situated at the confluence of outlets from Iztaccíhuatl and Popocatépetl (Fig. 3).The material was mainly tephra from past eruptions of Popocatépetl but deposited on the southeastern flanks of Iztaccíhuatl and then remobilized by the heavy precipitation event.
The highest climate station around Iztaccíhuatl which delivered data for many years was at about 4000 m a.s.l. on the western slope, but the data were not available for this study.In 1970 the 'Instituto de Geografía de la UNAM' published a climate map at a scale of 1:500 000, showing that September is the rainiest month.It is notable that at elevations above 5000 m a.s.l.precipitation can be >1000 mm/year with an average of 950 mm/year (Guerra-Acevedo, 2006).
Scenarios for water release from ice and snow leading to a potential lahar generation on Iztaccíhuatl can be: 1.An eruptive activity interaction at Iztaccíhuatl volcano.
3. A non-eruptive volcanic and/or a climatically related interaction.
Scenario 1 may have the highest impact but the probability of occurrence is low as Iztaccíhuatl is considered to be active but dormant.Deposition of ash from Popocatépetl onto the glaciers of Iztaccíhuatl as proposed in scenario 2 is a possibility that could occur at any future eruptive event at Popocatépetl if the wind direction is heading northwards as it occurred in the past (Siebe et al., 1996).This could lead to intensified melting of the glaciers and possibly triggering of lahars.The third scenario includes local melting of snow, ice and permafrost connected to geothermal volcanic influence (Abramov et al., 2008 in press) or seasonal melting.This can lead to subglacial, supraglacial or proglacial water reservoirs and moraine dammed lakes which are prone to outbursts.Such observations have repeatedly been made in volcanic (Russell and Knudsen, 2002) and non-volcanic glacial environments (Haeberli, 1983).
The importance of ice, firn, snow and permafrost in lahar formation at Iztaccíhuatl has not yet been investigated and there are no records of such events.In order to provide a related initial hazard assessment we modelled hypothetical flows originating from the glaciated areas.By establishing scenarios using different volumetric lahar magnitudes and DTMs of 30-90 m resolution, we can approximately define the areas around the entire Iztaccíhuatl edifice that could be affected in case of sufficient water release from Iztaccíhuatl's glaciers.

Data and methods
We used the following data and methods to analyze the glacier dynamics for evaluation of present and future ice conditions for lahar generation at Iztaccíhuatl.The satellite imagery used was acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) onboard the NASA Terra satellite.ASTER scenes consist of 3 visible and near infrared channels (VNIR; bands 1/2/3N; 0.52-0.86µm), a back-looking channel (band 3B; 0.76-0.86µm), 6 short wave infrared channels (SWIR; bands 4-9; 1.6-2.43µm), and 5 thermal infrared channels (TIR; bands 10-14; 8.13-11.65 µm) (Toutin, 2002).In false colour composite images (channels 3N/2/1; 9/4/3N and 10/9/4, corresponding to red/green/blue) of an ASTER scene from 17 March 2001, the large Little Ice Age (LIA) moraines are easily visible (White, 2002).Together with the LIA moraines shown in the 1:50 000 geologic map by Nixon (1989) and the 1:50 000 topographic map published by INEGI (1995), the moraines were mapped and the last glacial maximum extent was reconstructed (Fig. 4).The exact year of the LIA maximum for Mexican glaciers is not known but is assumed to be around 1850/1860 AD as is the case for comparable tropical glaciers in the Andes (Ramírez et al., 2001;Ceballos et al., 2006).
The oldest detailed map of the glaciers on Iztaccíhuatl was made by Lorenzo (1959).However, there remains an uncertainty concerning the exact year of data collection.Lorenzo (1959) used aerial photographs taken in 1945 as a guide but the field work was probably carried out in 1958.Therefore, the glacier extent derived from Lorenzo (1959) is assumed to represent the situation at around 1959.Additionally, there exists a 1:20 000 scale orthophoto map, published by INEGI (1983).The photographs used were taken in 1982 and one of the authors (H.Delgado) manually mapped the glacier extent for the same year.The best available data for glacier mapping consisted of aerial photographs, taken on 27 December 1994.They were orthorectified by the authors of this contribution and the glaciers were mapped with a GIS.Discrimination of snow and ice was simplified by analogue stereo interpretation.
The ASTER data with its medium spatial resolution (VNIR 15 m, SWIR 30 m, TIR 90 m; Toutin, 2002), is often not useful for repeated mapping of the changes in glacier extent at a high temporal resolution.However, within periods of several years, these low-cost images are adequate if the changes in glacial extent are significantly larger than the pixel size (15 m in this case).In our study, the time span between the aerial images from December 1994 to the ASTER-scene of March 2001, and the corresponding glacier retreat turned out to be sufficient to quantify changes.To get a distinct classification of the glaciers, band ratio images were calculated by dividing the VNIR channel 3N through the SWIR channel 4.This method takes advantage of the strong decrease of reflectance of snow and ice between the near-infrared channel 3N and the shortwave-infrared channel 4 as compared to other materials (Paul et al., 2002).The resulting ratio images depict ice and snow in bright contrast to the dark grey values for debris and rocks.
In late November 2004, two field visits were carried out.On 23 November 2004, the glaciers were mapped with the aerial photograph of 1994 as a reference.Due to inaccessibility and invisibility from the top, the glacier tongues at the steepest parts in the NW and E of El Pecho could not be mapped exactly.At these areas, the glacier extent mapped from the ASTER image of 2001 was taken as a maximum extent for 2004.However, the glacier tongues very likely retreated during the intervening 3.5 years, so that the glaciated area given for 2004 represents an upper limit for the real area in 2004.
Finally, the time series has been extended by two more recent ASTER-scenes from 15 March 2006 and 2 March 2007.Surprisingly, significant changes on El Pecho glacier could be detected even in the 1 year period between these two most recent images.
The DTM used in this study was taken from the Shuttle Radar Topography Mission (SRTM), acquired in February 2000 at approximately 30 m and 90 m resolution (Rabus et al., 2003).Here, the SRTM version with 3 arcsec spatial resolution, that is approximately 90 m, was used (SRTM3).In order to have a second DTM with a potentially better spatial resolution and to analyse the influence of the DTM accuracy on the modelled flows, DTMs with a grid size of 30 m and 60 m respectively were processed from the March 2001 Terra/ASTER channels 3N (nadir looking) and 3B (looking 27.6 • back from nadir; Kääb, 2005;Kääb et al., 2005).The spatial resolutions of 90 m, 60 m and 30 m for terrain elevation data is considered to be adequate for an initial, largescale hazard assessment of the entire Iztaccíhuatl edifice.Kaser and Osmaston (2002) made similar observations specifically on tropical and subtropical glaciers around the world.The observations emphasize the higher sensitivity of glaciers in lower latitudes to global atmospheric warming.
A regional influence on climate at Iztaccíhuatl by the significant emissions of the megacity of Mexico City (at a distance of only ∼40 km to the NW) might also play a role.

Glacier altitude and thickness
Along with the glacier area shrinkage, the former convex glacier shape on El Pecho has changed now to a thin ice body with a concave surface.Photographs from the mid-20th century show the crater of El Pecho summit covered by a massive convex ice cap.The altitude at the highest part of the glacier was formerly indicated as 5286 m a.s.l.(Lorenzo, 1959), but the current official altitude is 5,230 m a.s.l.(INEGI, 1995).If one excludes a change of the surveying height basis, this difference in elevation can be indicative of the decrease in glacier thickness.Within the last decades, the height of the original peak, which was in the centre of the glacier, decreased greatly.In 2004, the crater rims on the south and north of the El Pecho glacier emerged from the melting ice and finally became the new peaks of Iztaccíhuatl (A and B in Fig. 6).Handheld GPS measurements in November 2004 yielded a height of 5210 m a.s.l.(vertical accuracy in the order of tens of meters) for both crater rims and 5195 m a.s.l. for the centre of the glacier, where the main peak of the volcano was formerly located.The maximum height of this point on the SRTM3 DTM is 5211 m a.s.l. for the date of the acquisition in February 2000, with an absolute vertical accuracy of ±16 m (Smith and Sandwell, 2003).
The ASTER DTMs reveal maximum values of 5196 and 5176 m a.s.l. in March 2001 for the 30 m and 60 m DTM grids, respectively, with a vertical accuracy in mountainous terrain generally being at ±30 m (Eckert et al., 2005;Kääb et al., 2005).Considering that the summit area of El Pecho is represented as an almost flat region, a fact that often improves the height accuracy, and all the above independent measurements, we arbitrarily choose 5200 m a.s.l.±15 m as an approximate reference for the altitude of the centre of the El Pecho glacier between the years 2000 and 2004.The altitude in 2007 could be several meters lower due to further thinning of the ice on El Pecho.
Analysing the contour lines on Lorenzo's glacier map of 1959, we derive a height of 5270 m a.s.l. for the peak in the centre of the glacier.Compared to the value of ∼5200 m a.s.l.around the year 2004, the difference is 70 m within a ∼45 year period, corresponding to more than 1.5 m vertical decrease of ice per year.The loss of ice thickness on La Panza, situated 200 m lower in elevation, is also remarkable.During the same period, there is a vertical decrease of 30 m or 0.7 m per year.This value is confirmed by photographs of a large rock buried in the ice, which was barely visible in early photographs.Today it extends ∼25 m above the ice surface.

Climatic influence
In order to estimate the climatic influence on their retreat, Iztaccíhuatl's glaciers should be compared to similar Mexican glaciers.Unfortunately, this is difficult because there are only two other ice-clad volcanoes and both are not adequate for comparison: -The glaciers of Iztaccíhuatl's southern neighbour Popocatépetl (5452 m a.s.l.) have been strongly affected by the eruptive activity since 1994, in particular in 1997 and 2001 (Julio Miranda and Delgado Granados, 2003;Julio Miranda et al., 2005).
-The glaciers of Pico de Orizaba are much larger, reach up to 5675 m a.s.l. and therefore still benefit from an intact accumulation area (White, 2002;Delgado Granados, 2007).
Hence, we show in Fig. 5 the retreat of two tropical glaciers on the southern hemisphere with similar climatic conditions, comparable size and altitude range: Yanamarey glacier is situated in the Peruvian Cordillera Blanca at 9 • 40 ′ S at an altitude between 4750 and 5000 m a.s.l.(Kaser and Osmaston, 2002) while Chacaltaya glacier is at 16 • S between 5140 and 5360 m a.s.l. in Bolivia (Ramírez et al., 2001).The area of Iztaccíhuatl's individual glacier systems El Pecho and La Panza lies in between the two other glaciers shown, while the total glaciated area on Iztaccíhuatl is similar to the area of Yanamarey glacier during the past few decades.Although the glaciated area of the latter was a fourth of the one on Iztaccíhuatl in 1850 AD, during recent times Iztaccíhuatl's glaciers became smaller than Yanamarey.Chacaltaya glacier also shows a slower retreat than Iztaccíhuatl's glaciers.However, caution is needed for this comparison, because many basic factors such as high altitude precipitation and solar radiation due to cloud cover, etc. are not the same in these different regions.

Volcanic influence
The fact that the higher El Pecho glacier has a vertical melting rate approximately twice as fast as the lower La Panza glacier (1.5 m/yr.and 0.7 m/yr.respectively; see Sect.3.2) leads to the hypothesis that -beyond climatic effects -at least a part of the very fast vertical downwasting of the glacier on El Pecho could be induced by geo-or hydrothermal volcanic heat flow.
In the southern part of El Pecho glacier a notable feature was observed.On aerial images from 1977 a small dark point appears on the glacier of El Pecho.On the airphotos from 1994 there is a large pit at the same location, which was also observed during fieldwork in 2001.In 2004 the pit had developed into a crevasse-like opening of about 50 m in length (Fig. 6).The feature is situated in the middle of the glacier in a concave area where shear and extensional stresses in the ice are expected to be minor.This atypical location for crevasse development suggests that the form could be related to a locally confined heat source, e.g. a fumarole beneath the glacier, where geothermal heat flux is considerably higher (Welch et al., 2007).No vapour or sulfur smell was observed at this location but the influence of volcanic heat flow on glacier retreat at El Pecho is supposed to contribute to the ablation.The glacier may serve here as a giant calorimeter, as has been reported for example at Mt. Wrangell, Alaska, by Benson and Motyka (1978) and Benson et al. (2007).
Close to La Panza, Delgado Granados et al. ( 2005) installed a volatile trap in 1999 and measured Na + , K + , Ca + , Mg + , SiO 2 , SO 2− 4 , as well as F − , which suggest a magmatic source of the vapours.A strong sulfur odour could also be detected at two points around La Panza during field trips in November 2004.Delgado Granados et al. ( 2005) further report earthquakes that have been measured by seismic stations of the National Seismological Service (NSS) at Altzomoni and Tlamacas (A and T in Fig. 1).The origin of these seismic events was identified to lie beneath Iztaccíhuatl volcano (J.Pacheco, former head of the NSS; personal communication).Delgado Granados et al. (2005) reported temperature measurements between 0-10 • C at the ground and 0.1 m below the surface at La Panza region, and glacier-ice temperatures ranging from −2.0 • C and up to 0.5 • C at depths between −6 m and 0 m beneath the glacier's surface.The measurements indicate an increased volcanic heat flow and temperate glacier bed conditions.The existence of liquid water in englacial water reservoirs has already been reported by Álvarez and Delgado (2002) who carried out GPR (ground penetrating radar) profile measurements through La Panza glacier.

Current and future hazard scenarios
Observations of ablation and accumulation area over the past years indicated that the equilibrium line altitude (ELA) rose above the top of El Pecho, transforming the entire glaciers into an ablation area.The ELA only sinks for a short time period below the summit, and throughout the year ablation dominates on all the parts of the glaciers.If we linearly extrapolate the down-wasting trend of the ice, the glaciers will survive no longer than 2020 (Fig. 5).In a long term view this may reduce the risk of any volcano-ice interaction and subsequent lahar formation.Seasonal snow cover, permafrost, and some relict ice fields may however be preserved beyond this date.Firn and snow can be even more important for lahar formation than ice in case of interaction with eruptive volcanic activity (Pierson et al., 1990).
An eruption of Iztaccíhuatl volcano is currently unlikely.More probable is the development of hazards directly or indirectly related to glacier shrinkage, such as formation of a crater lake during or following glacier recession, release of sediment from basal moraines, slope destabilization by glacial retreat and permafrost thawing.For more detailed assessment of these factors, further data, including detailed GPR profiles of the El Pecho area and precipitation data of the summit region, should be collected.

Lahar model description
Two different flow models were used to simulate the potential lahar flow paths.The first one is the modified single-flowdirection (MSF) model by Huggel et al. (2003), which is further developed from the D8 flow direction method (Marks et al., 1984).Model input parameters are the starting cell location, a DTM and an empirical value for the H/L ratio (see below).The flow direction is calculated from each cell to one of its eight neighbours (D8), either adjacent or diagonal, in the direction of the steepest descent.In order to prevent stream width being represented by only one pixel, Huggel et al. (2003) implemented a function that enables flow diversion from the steepest flow direction.The more the flow diverts from the steepest descent direction, the higher the corresponding flow resistance.The resistance values are transformed into a semi-quantitative likelihood for each cell being affected by the lahar (Fig. 7).The modelled flow stops when it reaches a specified average slope value.This slope value is defined as the ratio between the vertical drop H at the flow initiation area and the horizontal distance L along the curving flow path to the lowest deposit (H/L; Fig. 8; Huggel et al., 2003).Actually, depending on the sediment concentration, lahars can be highly mobile and show rather attenuat-Fig.7. MSF model (Huggel et al., 2003) with an H/L of 0.19 (see Fig. 8) on the colour shaded SRTM3 DTM.Lahars were only modelled from the glaciers on El Pecho and La Panza to the east and west side so that the flows only represent lahars related to the glaciers.ing flow characteristics than abrupt stopping as granular debris flows.In the MSF model this is usually accounted for by choosing a minimum H/L value corresponding to a maximum flow reach.
The second model applied in this study was LAHARZ by Iverson et al. (1998).They examined 27 lahar paths on nine volcanoes and derived two semi-empirical equations predicting the inundated valley cross-sectional areas and the planimetric areas as a function of the lahar volume.The model assumes the simplification that there is only erosion above the limit of a 'high-energy cone' and only deposition below it (notwithstanding, erosion and deposition usually occur in conjunction).Similarly to MSF, the high-energy cone has to be defined by an H'/L' ratio which differs from the H/L used for MSF.The H'/L' ratio is the height and length from the highest point down to the limit of the high-energy cone where deposition starts, but not as far as to the lowest deposit (Fig. 8).As soon as the flow exceeds the limit of the high-energy cone, the model starts to calculate 3 valley cross-sections (at angles of 45 • , 90 • and 135 • to the flow direction) which are filled corresponding to the first semiempirical equation (Iverson et al., 1998).All inundated cells are marked as potentially affected by a lahar and added up until an area defined by the second semi-empirical equation is reached and the flow stops.The input parameters are potential lahar volumes, an appropriate H'/L' ratio for the border of the high energy cone, and a DTM.The ratio of sediment to water content of the simulated lahars corresponds to an average of those 27 examined lahars, from which the most diluted ones were excluded (Iverson et al., 1998).Hence, LAHARZ simulates typical debris flow-type lahars and the water content can not be varied.

Results of the MSF model
The MSF model was first run using the SRTM3 DTM from the glaciers on El Pecho and La Panza to the west and east side.H/L ratios for the lahar modelling between 0.2 and 0.15 were used, corresponding to published case studies (Siebe et al., 1996;Vallance et al., 2001;Capra et al., 2004).The volumes of debris flows correlate with the H/L ra-Nat.Hazards Earth Syst.Sci., 8, 559-571, 2008 www.nat-hazards-earth-syst-sci.net/8/559/2008/tio as demonstrated for lahars at Popocatépetl volcano by Huggel et al. (2008).That means that the greater the volume of a debris flow, the smaller is the corresponding H/L ratio (Legros, 2002).After the relation given by Huggel et al. (2008), the chosen values approximately represent volumes between 1×10 5 m 3 and 1×10 7 m 3 .Above 3700 m a.s.l., the results show a strong flow divergence, with high probabilities that the main drainage paths would be affected, and only low probabilities that broader lateral extents would be inundated.Below 3700 m a.s.l., the flows are channelized in deep ravines heading towards populated areas.On the west side the lahars reach some settlements at an H/L of 0.2 (namely San Rafael, San Antonio Tlaltecahuacán and Santiago Cuautenco), whereas the city of Amecameca would be affected with an H/L lower than 0.19.On the east side only San Agustín Atzompa would be hit partially if the H/L value falls below 0.17 (Fig. 7).

Results of the LAHARZ model
For the LAHARZ model the H'/L' ratio, which defines the border of the high-energy cone, had to be determined (Fig. 8).To this end, we analysed the lahars at Popocatépetl which occurred during its eruptions in 1997 and 2001.In both cases, the village of Santiago Xalitzintla SE of Iz-taccíhuatl (Fig. 9) was hit by the lahars with volumes of about 4×10 5 m 3 (Capra et al., 2004;Julio Miranda et al., 2005).We calculated several high-energy cones to ascertain that the modelled LAHARZ flow of 4×10 5 m 3 just reached the village of Santiago Xalitzintla as it happened in 1997 and 2001.Within these simulations we obtained a best-fit-value of 0.28 for the H'/L' at Popocatépetl.Assuming similar conditions at Iztaccíhuatl, we used the same H'/L' value for the LA-HARZ model runs on Iztaccíhuatl.For comparison, Iverson et al. (1998) used an H'/L' of 0.23 at Mount Rainier in the US, and Davila et al. (2007) used a much lower value of 0.04 for Colima volcano in Mexico.The latter value is probably more appropriate for hyperconcentrated flows than for debris-flow-like lahars.With such a value the high-energy cone would lie completely in the alluvial plane and hence, it seems far too low for Iztaccíhuatl volcano.
We assume an average ice thickness of 30 m (GPR data from Delgado Granados et al., 2005) over the whole glaciated area.Consequently, the total ice volume in 2007 was roughly 8.2×10 6 m 3 .On the basis of ice availability, as well as the volume of the 1997 Popocatépetl lahar (4×10 5 m 3 ; Julio Miranda et al., 2005), we defined hypothetical lahar volumes at Iztaccíhuatl as 1×10 5 m 3 , 5×10 5 m 3 and 1×10 6 m 3 .These volumes correspond to flows generated from melting of 0.5, 2.4 and 4.9% of the total ice volume on Iztaccíhuatl, assum-ing a 40% water content in volume (boundary between debris flows and hyperconcentrated flows; Pierson and Scott, 1985;Pierson, 1986;Vallance, 2000;Lavigne and Thouret, 2002).The mass of ice lost is typically larger than the amount of water contributing to a possible lahar because not all meltwater is incorporated into the flows (Thouret, 1990;Julio Miranda et al., 2005).These volumes are plausible for the possible scenarios without any additional water from rainfall (as described in 1.3.)and lie within the range of modelled lahar volumes at Concepción Volcano in Nicaragua (Vallance et al., 2001) and at Colima in Mexico (Davila et al., 2007).Larger or more diluted lahars, which also may reach farther, are assumed to develop only as a combination of glaciervolcano interaction processes and rainfall events with shallow landsliding and additional sediment entrainment.This case seems to have a low probability.For that reason and the fact, that the sediment content can not be varied in LAHARZ (see Sect. 4.1), we do not present model results for lahars larger than 10 6 m 3 .
On the west side, the LAHARZ-flows reach San Rafael, San Antonio Tlaltecahuacan and Santiago Cuautenco with volumes less than 5×10 5 m 3 while Amecameca would only be reached with more than 5×10 5 m 3 .In the eastern part, the same effect as with the MSF model is observed.The material stays in the deep and narrow channels but as a consequence, the flows travel farther.Lahar volumes greater than 5×10 5 m 3 reach San Agustín Atzompa in the NE.The affected town is situated at the mouth of a valley that drains from Iztaccíhuatl, similarly to Santiago Xalitzintla, which is situated completely in an erosion trench.However, most of the towns in the east are situated safely on the terraces in between the canyons, such that the ravines build a reasonable protection from lahars up to hypothetical volumes of 1×10 7 m 3 .4.4 Influence of the DTM on lahar models Huggel et al. (2008) have shown that the quality and resolution of the DTM can have an important effect on the flow paths and runout distances of simulated lahars.For the MSF model, the changes of the runout distance between the runs on the SRTM-DTM and the ASTER-DTM are minor.However, small differences in the model terrain morphology can determine if a flow will find a way over a moraine or a narrow crest line.If on one DTM the flow takes a different path than on the other DTM, the subsequent flow inundates a formerly unaffected valley providing completely different results.In this study, significant differences resulted from flows modelled on the ASTER DTMs with 30 m and 60 m resolution.They were of an unsatisfying quality, mainly due to large data gaps which had to be interpolated.The modelled lahars flowing through these interpolated areas are mostly of limited practical value.Comparison of individual parts of the ASTER-DTMs to the topographic maps and to the SRTM3 DTM helped however to locate a number of crucial points where the natural sidewalls of the downstream valleys are low or narrow, and larger lahars effectively could break through.
Concerning flow channel deviation, the LAHARZ model generally has fewer problems than the MSF because only the lower parts below the high-energy cone are simulated.There, the ravines are wider and better represented by the SRTM3 and ASTER DTMs than the areas in the summit region.However, analysis of the calculated lahars on the SRTM3 and ASTER DTMs revealed significant variations in the flow path at a problematic zone on the east side of Iztaccíhuatl (Fig. 10).To assess the quality of each DTM, contour lines were derived and compared to the topographic map and the ASTER VNIR image.It was found that the ASTER DTM is defective in this area (mainly due to stereo image matching problems), forcing the modelled lahar on a wrong flow path.This is visible in the subtraction DTM of the ASTER-and SRTM-data (Fig. 10).Two erroneous areas with exaggerated heights up to 100 m (green) fall directly within the flow path, deflecting the lahar stream.Despite the abbreviation of the path, the lahar modelled on the ASTER DTM is still stopping earlier.This effect was observed for each modelled flow and we interpret grid cell size to be an important influence on runout distance calculated by the LAHARZ model.The flows inundate a wider path on the ASTER DTM (60 m) than on the SRTM3 DTM (90 m) and hence on the ASTER DTM, less material is transported downwards.This leads to significantly shorter distal hazard zone boundaries.A lahar modelled on the SRTM3 DTM at the northeastern slope and with a volume of 1×10 6 m 3 affects the central part of San Agustín Atzompa while the same volume on the ASTER DTM stops 2.5 km upstream and represents no hazard for the settlement (Fig. 9 and 10; this effect was also observed by Huggel et al., 2008).Stevens et al. (2002) found similar effects concerning the flow length on different DTMs but not to that extent.Their terrain data were derived from topographic maps and an airborne SAR (TOPSAR by NASA) with 20 m and 10 m horizontal resolution respectively.In contrast to this, the hazard zonation of Davila et al. (2007) at Colima volcano in Mexico reveals higher lateral spreading and a 1.5 km shorter runout distance for the coarser DTM.Therefore we deduce that the behaviour of the LAHARZ model is highly dependent not only on the spatial resolution, but also on the acquisition method of the DTM.

Conclusions
The glacial recession on Iztaccíhuatl has been pronounced over the past decades and is likely to lead to total deglaciation within few decades, probably before ∼2020 AD.This trend is in accordance with observed worldwide glacier retreat in tropical regions that is attributed to atmospheric warming.In the case of Iztaccíhuatl glacier retreat may also be affected by an increase in the volcanic heat flow.The crevasse-like opening on El Pecho glacier may be the chimney of an active geo-or hydrothermal vent and should be investigated further, as should the diffuse degassing fields on La Panza.Gas emission measurements, GPR profiles (which have only been performed for one glacier on La Panza by Delgado Granados et al., 2005), and ground temperature measurements could provide a better understanding of the volcanic activity of Iztaccíhuatl, which has been under dispute during the last decades.The volcano is not necessarily ready to become active again but there appear to be some changes in the volcano's hydrothermal system taking place.This could promote the development of subglacial water reservoirs and their possible outbursts may provide a sufficient amount of water for lahar formation.Within a 15 km radius of the summit, a population of approximately 210 000 lives that could be affected by lahars, depending on the lahar volume.
Glacier-related hazards at Iztaccíhuatl are expected to decline over the next decade due to the rapid loss in glacier area and volume, or even total glacier disappearance.However, seasonal snow will still pose a hazard in the case of interaction with the volcano, and several new hazards are not related to volcanic activity: (1) the development of a crater lake at El Pecho due to climate related glacier retreat, which may be prone for outbursts and (2) destabilization of the over-steepened upper flanks due to glacial recession and permafrost thawing, which both could favour the initiation of lahars or landslides.The large amounts of unconsolidated debris which have been revealed by the retreat of the glaciers could be mobilized by heavy precipitation during the rainy season.Accurate precipitation data from the summit of Iztaccíhuatl do not currently exist, but a meteorological station will be installed soon by one of the authors.However, raintriggered lahars are more common on lower slopes (Scott et al., 2001).The hazard potential of this individual trigger mechanism should be investigated in a separate study, especially concerning a future eruption of Popocatépetl with following ash deposition on the slopes of Iztaccíhuatl that could be remobilized by heavy precipitation.
Concerning the safety of the surrounding villages, there is an important morphological difference between the western and eastern flank.Deep incised valleys serve as natural stream channels, and therefore the location of villages relative to the channel is crucial.On the west side, the ravines directly lead to the populated alluvial plane, while on the east side most of the villages are built on the higher, and hence safer, terraces between the ravines.More detailed studies are necessary to evaluate the need for mitigation measures such as deflection dams, retention basins or a warning system that would increase the safety of the settlements around Iztaccíhuatl.
The SRTM3 DTM represents the main drainage channels in a way that is sufficient for a large scale hazard assessment.In the present case, it provides more realistic results than the ASTER DTMs.However, the coarse resolution of the SRTM3 DTM allows less lateral deposition of material, a fact leading to exaggerated runout distances of the LAHARZ flows.As a consequence, for lahar hazard assessment, we emphasize the necessity of using different flow models in conjunction with mid-and high-resolution DTM data from various sources.Flow modelling on high-resolution DTMs would be helpful for detailed mapping of the hazard zone boundaries on the basis of the presented initial hazard assessment, but each simulation should be critically compared to results based on calculations on other DTMs.

Fig. 1 .
Fig. 1.Upper right: Location of the most active volcanoes and the study area within Mexico.The coloured shaded relief showing Mexico City, the study area and Puebla is derived from SRTM3 elevation model.The dotted area shown as 'prehispanic lahars' represent the approximate extent of areas flooded by repeated lahar events from Popocatépetl, Iztaccíhuatl and La Malinche after Siebe et al. (1996).

Fig. 2 .
Fig. 2. The silhouette of Iztaccíhuatl resembling a lying woman as believed by indigenous Mexican people, seen from the west: La Cabeza (head), El Pecho (chest), La Panza (belly), Las Rodillas (knees) and Los Pies (feet).El Pecho and La Panza are still glacier clad.Picture taken on 21 November 2004 by D. Schneider.

Fig. 3 .
Fig. 3. View of the streets of Santiago Xalitzintla showing the deposit of a rain-triggered lahar generated on Iztaccíhuatl's southeastern flanks with deposits erupted from Popocatépetl.Picture taken on 20 June 1994 by H. Delgado Granados.

Fig. 4 .
Fig. 4. Glacier outlines for eight different dates.Contour line interval is 50 m and the background image is the ASTER VNIR scene of March 2001.

Fig. 5 .
Fig. 5. Decreasing glacier extent at Iztaccíhuatl (Mexico, 19 • N) compared to the glaciers Yanamarey (Peruvian Cordillera Blanca, 16 • S; data from Kaser and Osmaston, 2002) and Chacaltaya (Bolivia, 9 • 40 ′ S; data from Ramírez et al., 2001) through time (years AD).On Iztaccíhuatl volcano, the two glaciated parts 'El Pecho' and 'La Panza' are shown separately and summarized as 'Iztaccíhuatl'.These two parts were connected in 1850 AD and hence there is no separate data for 'El Pecho' and 'La Panza' for that date.Note the break in the y-axis and the possible future development of the glaciated area (dotted line).

Fig. 6 .
Fig. 6.El Pecho Glacier from south to north with the southern crater crest line (A) and the former ice clad northern crest line (B).A large crevasse is situated in the middle of the crater glacier (C).Note the current concave glacier shape in contrast to the former convex reconstructed surface of 1959, indicated by the dashed black line.The height difference at the centre of the glacier is around 70 m.Picture taken on 22 November 2004 by D. Schneider.

Fig. 8 .
Fig. 8. Schematic view of a volcano and the difference between the H/L ratio used in the MSF model and the H'/L' required for LA-HARZ.The division into a higher part where lahars erode material and a lower section where deposition dominates is needed to define the high-energy cone used in LAHARZ.

Fig. 9 .
Fig. 9. LAHARZ(Iverson et al., 1998) model on the colour shaded SRTM3 DTM.The dark blue line around the summit is the high-energy cone with an H'/L' of 0.28 (see Fig.8).The lahars were modelled for the same drainage paths as with the MSF model.Red, orange and yellow represent different lahar volumes with corresponding degree of hazard(high, medium, low).LAHARZ does not model erosion in the upper part above the high-energy cone.

Fig. 10 .
Fig. 10.Height difference of ASTER and SRTM DTM (ASTER -SRTM).Green cells indicate areas where the ASTER DTM is higher than the SRTM DTM and the blue cells vice versa.Contour lines are derived from the SRTM DTM at an interval of 50 m.The black lahar was modelled on the SRTM DTM while the blue line (high energy cone with an H'/L' of 0.28) and the colored lahars were calculated on the ASTER DTM.Problems with the ASTER DTM mainly occur in narrow ravines which lead to incorrect deflection of the modelled lahar.

Table 1 .
Glacier areas in km 2 and percentaged in relation to the 1850 extent.The area listed as 'Iztaccíhuatl' is the sum of the individual 'El Pecho' and 'La Panza' glacier systems.
(Zemp et al., 2007)ollected data reveals strong changes in glacier extent during the past decades on Iztaccíhuatl.Figures 4 and 5 show the glacier retreat, and in Table1the corresponding absolute and relative values are given.As a simplification, the different glaciers were grouped into two systems: El Pecho and La Panza.The areas of the smaller glaciers Atzintli and San Agustín (now disappeared) south of La Panza are counted into the La Panza system.If one takes the 1850 extent (6 369 000 m 2 ) as a reference, only 4.3% (273 000 m 2 ) of the glaciated area remained in March 2007.During the more recent 48-year period between 1959 (1 369 000 m 2 ) and 2007 (273 000 m 2 ), around 80% of the glaciers disappeared.This is in accordance with the general worldwide trend of glacier shrinkage(Zemp et al., 2007)but the retreat rate is extraordinarily high compared to the global average.