Introduction
Volcán Villarrica (39°25′12″ S, 71°56′27″ W; 2847 m a.s.l.; Fig. 1) is a highly active ice-capped volcano, which is characterized in historical times mainly by mild strombolian activity (Reference González-FerránGonzález-Ferrán, 1995; Reference Lara, Lara and ClaveroLara, 2004), permanent degassing and periodic explosions, with the lava lake remaining at a high level (90–180 m below the surface), at least since 1984. It is also very sensitive to magmatic-conduit activity (Reference CalderCalder and others, 2004; Reference Witter and DelmelleWitter and Delmelle, 2004). Concentrations of acid gases measured at the summit of the crater have been recognized as a hazard to climbers on the volcano, who may be exposed to concentrations above the limits defined by the US National Institute of Occupational Safety and Health (Reference Witter and DelmelleWitter and Delmelle, 2004). Its historical eruptive activity (Reference Petit-Breuihl and LobatoPetit-Breuihl and Lobato, 1994; Reference Lara, Lara and ClaveroLara, 2004) indicates a low frequency of large explosive eruptions (volcanic explosivity index (VEI) 3–4). More than 50 eruptive events have been documented since 1558 (Reference Petit-Breuihl and LobatoPetit-Breuihl and Lobato, 1994). The latest most violent eruption took place in 1971/72 when lava flows were generated, as well as 30–40 km h−1 laharic flows (Reference Naranjo, Moreno, Lara and ClaveroNaranjo and Moreno, 2004) descending towards Lagos Villarrica and Calafquén (Fig. 1). Lahars produced by eruptions of Volcán Villarrica in 1948/49, 1963/64 and 1971/72 have resulted in the deaths of more than 75 people (Reference SternStern, 2004), and are considered the main hazard of the volcano (Reference MorenoMoreno, 2000). The volcano is covered by a 30.3 km2 glacier (Reference RiveraRivera and others, 2006), mainly distributed towards the south and east where the main glacier basin (Glaciar Pichillancahue–Turbio (17.3 km2); Fig. 2), composed of partially ash/debris-covered ice, is located, partially infilling a volcanic caldera depression (Reference Clavero, Moreno, Lara and ClaveroClavero and Moreno, 2004). The energy balance of this glacier has been monitored since 2003 (Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others, 2007), and global positioning system (GPS) as well as radio-echo sounding (RES) measurements were carried out in January 2005 (Reference RiveraRivera and others, 2006).
One of the most interesting direct effects of volcanic activity on the overlying glaciers is the ice cracking or crevassing observed before or during eruptive events (Reference KlohnKlohn, 1963; Reference Fuenzalida-Ponce and González-FerránFuenzalida-Ponce, 1976; Reference Fuentealba, Riffo, Moreno and AcevedoFuentealba and others, 1984; Reference González-FerránGonzález-Ferrán, 1995). Ice cracking has been detected by seismicity (Reference Métaxian, Araujo, Mora and LesageMétaxian and others, 2003), but here we propose observing possible ice cracking by obtaining daily photographs of the ice. This process could be important for ice flow and for the hydrological balance of the glacier, as crevasses are the main pathway for meltwater to enter the glacier’s en-and subglacial drainage system (Reference Fountain and WalderFountain and Walder, 1998; Reference Fountain, Jacobel, Schlichting and JanssonFountain and others, 2005).
A second direct effect of volcanic activity is the deposition of pyroclastic materials on top of the glaciers, changing the albedo and affecting the energy balance by providing insulation from atmospheric heat and insolation where particles are large or the cover is continuous (Reference Adhikary, Yamaguchi and OgawaAdhikary and others, 2002; Reference Kirkbride and DugmoreKirkbride and Dugmore, 2003). Studies of tephra thermal properties and glacier melt rates at Volcán Villarrica (Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others, 2007) identified a very low thermal conductivity of 0.35 W m−1 K−1 in the lapilli tephra which blankets most of the ablation zone. Furthermore, a critical thickness of just 5 mm of tephra cover was found to reduce the melt rate of the buried ice compared with a bare surface. Consequently, volcanically produced materials appear to have a large positive impact on the mass balance of glaciers on Volcán Villarrica, due to an extensive mantle of insulating tephra in the lower ablation zones. In this study, quantitative assessment of the total area of tephra cover and tephra-free ice slopes and snow albedo has been made through terrestrial photography.
Aims and Methods
The main aim of this paper is to analyse albedo variation and its impact on melt on Glaciar Pichillancahue–Turbio of Volcán Villarrica by means of terrestrial photography. Previous measurements of the glacier areal variations (Reference RiveraRivera and others, 2006) have been updated until 2007, to allow the assessment of recent ice retreat on the volcano’s glaciers.
Meteorological data
The meteorological data were obtained by an automatic weather station (AWS) located during the summer of 2003/04 on the surface of the glacier at 1933 m a.s.l. (near the equilibrium-line altitude, ∼2000 m a.s.l.) (Reference RiveraRivera and others, 2006). During the winter the AWS was moved to a location on a rock outcrop at 1925 m a.s.l., next to the fixed camera (Figs 2 and 3). The AWS recorded incoming and reflected shortwave radiation, net all-wave radiation, air temperature, air humidity and wind speed at 2 m above the surface with hourly mean values recorded on a Campbell CR10 data logger. Albedo was measured using a Kipp & Zonen CM6B albedometer sensitive to radiation in the range 0.3–2.8 μm, with the sensors mounted in a surface-parallel plane. Details of the collected variables and technical details of the sensors at the AWS are described by Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others (2007).
Oblique photography
For measuring albedo, surface changes and tephra cover, an automatic camera (Canon EOS 300D) was installed at the upper part of the rim surrounding the main volcano edifice (upper part of the caldera, Fig. 2), from where daily photographs of the glacier were obtained. The camera has a 6.3 megapixel complementary metal oxide semiconductor (CMOS) sensor and recorded the images on a 2 GB flash-card memory. It was fitted with a high-quality, fixed 24 mm focal length lens to minimize optical distortions. This camera sensor has some sensitivity in the near-infrared spectrum, at least to 1000 nm, and beyond that point if the infrared filter is removed (experimental tests by the authors). The camera was inserted into a Pelikan sealed box where it was controlled by a Canon timer. The system was powered by 12 V batteries which were fed by a solar panel installed nearby.
Conventional photography is a powerful medium for collecting and storing information. If this information can be located precisely in space, then photography becomes a powerful tool for quantitative analysis. Here we use a tool for georeferencing oblique photography developed by Reference CorripioCorripio (2004), using a single image and a digital elevation model (DEM). The accuracy of the technique depends on the accuracy of the DEM and on the quality of the photographic image, especially the degree of distortion and aberration produced by the lens. The technique does not produce elevation data and requires an existing DEM. What the tool does is to locate the geographical position of every pixel in a photographic image. It is therefore useful in mapping land cover and assessing changes in surface cover. The technique follows standard procedures for perspective views in computer graphics and photogrammetry (Reference Slama, Theurer and HenriksenSlama and others, 1980; Reference FiumeFiume, 1989; Reference Foley, van Dam, Feiner and HughesFoley and others, 1990). It creates a virtual photograph of the DEM that can be scaled to the resolution of the photographic image to establish a mapping function between pixels in the photograph and gridcell points. This allows the exact position of pixels to be located on the oblique image. The georeferencing process consists of a viewing transformation applied to the DEM, in which the coordinates of every gridcell are firstly translated to refer them to a coordinate system with the origin at the camera position. A transformation is then applied according to the viewing direction and focal length of the camera. This results in a three-dimensional set of points corresponding to the cells in the DEM as seen from the point of view of the camera. Finally, the resulting viewing transformation is projected into a two-dimensional viewing ‘window’, corresponding to the area of the film and scaled proportionally. The process is explained in detail by Reference CorripioCorripio (2004); it has been coded in Interactive Data Language (IDL) under the Creative Commons license and is freely available from the authors.
Using this technique the evolution of the snow cover was mapped accurately during the acquisition periods. Once the photograph is georeferenced, the reflectance values are normalized according to the viewing geometry, the angle of incidence of the sun’s rays on the slope, the ratio of direct to diffuse radiation, the atmospheric transmittance between the pixel location and the position of the camera, and the effect of radiation reflected from the surrounding slopes. The final result is a map of normalized reflectance values, or relative albedo. The atmospheric transmittance was calculated using the radiative-transfer model MODTRAN (Reference Berk, Bernstein and RobertsonBerk and others, 1989) and general knowledge of the local atmospheric conditions from the AWS. The photographic image was mapped to a resampled DEM from the Shuttle Radar Topography Mission (SRTM, acquired by JPL (Jet Propulsion Laboratory, Pasadena, CA, USA)/NASA in 2000). It is important to note that the resampling procedure to a 10 m resolution DEM cannot increase the resolution of the original DEM, at 90 m, but allows the extraction of more information from the photography within the known spatial limits of the original DEM gridcell. The DEM was not contemporary with the photographic image acquisition. This may introduce errors into the results, as the slope and aspect may be slightly different from one date to another, and the snow accumulation may change. It is beyond our financial capability to acquire a high-resolution DEM for every campaign, but we are working on the design of some terrain measurements that allow the precise evaluation of the errors incurred by using slightly old DEMs, such as simultaneous measurements of albedo at different points and comparison with those derived from georeferenced photography.
Unfortunately this project was a constant battle against the elements. Despite testing the kit under laboratory conditions, the camera case was twice flooded during severe weather conditions, which also prevented collection of data on one occasion. The camera set-up was finally destroyed by a lightning strike in March 2007.
GPS survey
Several kinematic and static GPS surveys were conducted to map the extent of the glacier and to georeference tie points for photogrammetric purposes (Fig. 3). Javad’s dual-frequency GPS receivers and antennas were used exclusively at a sampling rate of at least 2 s and an elevation cutoff angle of 10°. At all times a nearby reference station mounted on bedrock was occupied to provide geodetic quality measurements for differential positioning at the centimetre level. Thereby, baseline lengths never exceeded a few kilometres. GPS data were post-processed using the Waypoints GrafNav software package, version 7.70, where precise ephemeris and clock information provided by the International GNSS (global navigation satellite systems) Service (IGS, final products) was incorporated. In January 2005 the local reference station was linked to the SIRGAS (Geocentric Reference System for the Americas) which is the regional realization of the International Terrestrial Reference Frame.
Satellite imagery
In order to update the glacier variations at the volcano, several satellite images and aerial photographs were used (Table 1). All the satellite images were geolocated and orthorectified using the regular IGM (Instituto Geográfico Militar) cartography and available DEMs (e.g. SRTM and airborne synthetic aperture radar (AirSAR)) (Reference RiveraRivera and others, 2006). Once the satellite images were orthorectified, classification procedures based on spectral band ratios were applied to account for the glacier extent and snow/ice/debris classification (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002). Aerial photographs were stereoscopically analysed and the resulting information was transferred with a zoom transfer scope (ZTS) to the regular cartography, with ice fronts being digitally compared to the satellite (Reference Benson and FollettBenson and Follet, 1986). All glacier limits were analysed using GIS (geographical information system) commercial software, such as IDRISI 32 for Windows, ArcInfo version 8.0.1 and PCI Geomatica, allowing an accurate estimation of areas and frontal changes (Fig. 2).
Energy-balance and melt modelling
We ran a modelling experiment to assess the effect of tephra deposition on the glaciated surface of the volcano and its influence on mass balance and runoff. Using meteorological data for 12 days in January 2007 and the corresponding albedo variations derived from the photographic images, a distributed energy-balance model was applied to the eastern slopes of the Volcán Villarrica. Thus, we could assess the effect of diminished albedo due to volcanic ash deposition on the mass balance and meltwater production from the volcano glacier. In January 2007 there were no available wind-speed data from the AWS. These data were derived from the NOAA (US National Oceanic and Atmospheric Administration) archived 1° resolution aviation (AVN) model outputs for the gridcell corresponding to Volcán Villarrica. The data were chosen at a 750 hPa pressure level, which corresponds to the middle height of the volcano’s cone. The model applied (SnowDEM: Snow Distributed Energy-balance Model) is explained in detail by Reference CorripioCorripio (2003a). It is a distributed, multilayered energy-balance model, which takes into consideration radiative fluxes, heat interchange with the atmosphere, evaporation or sublimation and heat flux due to precipitation. Shortwave radiation is evaluated according to a detailed parametric radiative transfer model of the atmosphere plus terrain effects (Reference CorripioCorripio, 2003b, Reference Corripio2004; Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others, 2004). Longwave radiation is estimated from atmospheric humidity and temperature plus terrain parameters and ground/snow emissivity according to geological characteristics if data are available. Tests on complicated surfaces in the Andes, such as penitentes, testify to the ability of the model to reproduce snow surface temperature (Reference Corripio, Purves, De Jong, Collins and RanziCorripio and Purves, 2005). Application to a watershed in the Alps for the estimation of meltwater runoff gave values within 6% of measured runoff for a runoff gauge with a 10% accuracy (unpublished data).
Results and Discussion
The frontal and areal glacier variations of Glaciar Pichillancahue–Turbio have been updated up to 2007 (Table 2). Both the main arms of the glacier have continued retreating at similar rates, as indicated by frontal length changes measured in recent decades (Reference Casassa, Acuña, Zamora, Schliermann, Rivera, Lara and ClaveroCasassa and others, 2004). These glacier tongues are totally debris-covered and only at steep flanks is bare ice visible, due to the backwasting ablation process (Fig. 2). In spite of the insulation provided by the ash and debris covering the ice, the glacier has lost significant area in recent decades, much more than other debris-free glaciers located on top of active volcanoes (Reference RiveraRivera and others, 2006). The present extent of Villarrica’s glaciers is shown in Figure 2, and recent variations are summarized in Table 2.
In order to monitor the glacier at a daily resolution, we used oblique terrestrial photographs that were georeferenced to a DEM with the help of accurate ground-control points (GCPs) measured on the glacier surface, as shown in Figure 4. This tool can be applied to precisely locate snow features on the surface of the glaciers in areas that are very difficult to access. The eastern upper side of the volcano is constantly swept by ash falls and toxic gases from the crater fumaroles. This made direct surveying impossible without appropriate protective clothing and breathing masks. Thus remote sensing is a more convenient alternative, as shown in Figure 5. The figure shows the changing positions of crevasses on the upper section of the volcano. The left image is from 25 March 2006 and the right image is from 14 January 2007. The derivation of flow rates and surface variations using this approach is currently being studied. This preliminary demonstration is given to illustrate the potential of the technique for high-temporal-resolution surface monitoring in hazardous or inaccessible environments.
The pattern of tephra dispersion is clearly visible on the georeferenced image of the volcano of 25 December 2005 (Fig. 6). The conical shape of the georeferenced image is due to the field of view of the camera. There is a thick tephra layer around the crater rim and a band of darker snow to the southeast following the prevailing west and northwest winds. The amount and area of tephra deposition depends on the intensity of fumarolic activity and the concurrent winds. It would be possible to model wind flow across the volcano, but it is far more difficult to predict the spatial and temporal variability of fumarolic activity. It is therefore very difficult to anticipate the ash dispersion pattern. The only solution is to observe it at relatively high temporal resolution and incorporate the results in any model of the glacier surface.
Here we present the results for a set of 11 images during a clear-sky period from 5 to 16 January 2007. The albedo of the visible area of the volcano was derived from the photographic images and incorporated into the energy-balance model. The albedo derived from the images reveals increasing values toward the fringes of the visible area (the northern and southern slopes). We believe this is a realistic result, as those slopes are away from the prevailing winds and suffer less tephra pollution by volcanic fallout. In fact, a visual inspection of the northern slopes while skiing Villarrica a week earlier showed thick ash layers on the uppermost section of the volcano near the crater and clean, metamorphosed granular snow below that point to about 1800 m a.s.l. On the upper section snow was fine-grained, interspersed with small wind deposits of highly broken precipitation particles (Reference ColbeckColbeck and others, 1990); on the lower section it consisted of larger snow grains with high water content.
Previous work suggests that the critical ablation threshold is passed as soon as the tephra forms a continuous layer at the surface, with insulation and melt reduction overriding the influence of albedo lowering (Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others, 2007). In our analysis we assume the critical debris thickness to correspond to an albedo of 11% (based on the broadband albedo of andesitic basaltic tephra, ASTER spectral library, http://speclib.jpl.nasa.gov/), below which snowmelt is reduced relative to bare-snow surfaces. Spurious high values may be due to errors in the precise boundary between tephracovered and bare-snow areas. These errors are of the order of two pixels of the original DEM resolution at the borders of the image, or ±180 m. Areas where the viewing angle is very shallow have been masked, but sub-pixel variation in slope may add to the error. It is unlikely that the precision can be increased without the use of photogrammetric cameras and a very up-to-date DEM. The resulting ablation map (Fig. 7) shows high spatial variability. This variability may be enhanced by post-depositional reworking of the tephra layers, and facilitated by a positive feedback, as surface particles will tend to aggregate while melting on concave surfaces (Reference DrewryDrewry, 1972), which are initially caused by differential melt. Of particular note are areas of locally enhanced melting downwind of exposed ash banks in the lower third of the image, and accelerated melting over large areas downwind of the crater, which are exposed to almost continuous airfall tephra deposition. These irregularities are also discernible on any transect along and across the eastern slopes of the volcano, as shown in Figure 8. This figure shows the modelled differences in accumulated ablation for a transect along the 2500 m elevation contour and along a transect from the summit crater towards the camera. Along the elevation isopleth the pattern is approximately symmetrical, with minimum values towards the slopes that are less subject to ash deposition and a maximum on the southeast slopes. The horizontal gradient is rather large, with differences in ablation of several centimetres over a few tens of metres. On the vertical transect we can observe a clear anomaly, as ablation increases with elevation and reaches a local maximum near the crater (distance 0–500 m), where ash depositions are more intense. It then follows an irregular but decreasing trend downslope, as ash gets more dispersed with distance from the source. There are two peaks ∼5 km from the crater, associated with older ashes that are resurfacing after the overlying snow has disappeared. In these areas the tephra is probably thick enough to insulate the underlying ice and reduce melting. Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others (2007) measured a mean daily melt rate of 46 mm w.e. at three stakes set on snow with mean albedo 0.51 in the vicinity of the AWS, under similar meteorological conditions in the second half of January 2004. The measured melt rate (12 day cumulative melt = 506 mm) corresponds well with modelled values using photographic-derived albedo in January 2007 over large areas of snow with relatively light tephra cover (Figs 7 and 8).
The effects of variable deposition, and redistribution of fallen tephra, on spatial patterns of melt are illustrated by comparison with modelled melt rates when these spatial variations are neglected. Figure 9 shows the differences in modelled cumulative melt rates for the same period, replacing photographic-derived albedo with the Reference Brock, Willis and SharpBrock and others (2000) ageing-curve albedo parameterization. This parameterization calculates albedo as a function of temperature since the most recent snowfall, and assumes albedo decay is caused by snow metamorphism and the build-up of impurities over time, but does not account for spatial variation in snow impurity content. While the model with parameterized albedo is able to account for some of the along-glacier (vertical) variation in melt rates associated with slower snow metamorphism at higher elevations, cross-glacier variation in melt rates, other than that due to aspect and shading, is not incorporated and, in particular, the high spatial variability in melt rates immediately below the crater and in the vicinity of exposed ash banks is missed (Fig. 9). This demonstrates that using photographic-derived albedo provides an improved depiction of the actual spatial variation of surface albedo.
The differences in cumulative ablation between models accounting or not for tephra deposition reach peak values of >40 cm in some spots, but these peak values are probably an overestimation. The errors in this computation are due to (1) sub-pixel slope and albedo variation; (2) local areas of very shallow viewing angle where the average DEM slope is steeper and (3) projection of cones and mounts into adjacent grids due to almost parallel viewing angles. Therefore only the areas where the viewing angle of the camera is >60° are reliable. This results in differences in ablation of ∼20– 30 cm between model runs considering ash cover or clean snow. While at the end of the accumulation season the surface of the glacier is relatively smooth, at the end of the ablation season it presents many concavities a few metres in depth and a few tens of metres in diameter (Fig. 10). The negative values are areas where the photographic-derived albedo is higher than the parameterized albedo. This is likely to happen in small concave areas which get more shading from the sun. These small concavities are not registered in the DEM; a more precise estimation of their albedo would require a higher-resolution DEM. The irregular snow surface is likely to be due to a combination of effects from glacier flow and crevasse formation together with differential ablation.
The increase in ablation of snow due to tephra impurities was estimated by Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others (2007) to be 8% at the location of the AWS. The same authors recorded rates of melting similar to those produced here for the areas with higher ash depositions. Ablation certainly seems higher further up the slope, as the albedo is lower. This can also be seen from the longitudinal transect (Fig. 8). The average modelled increase in ablation for the pixels visible from the camera viewpoint is 13.6% in this study. It is questionable whether this increased ablation and related meltwater production might affect the glacier dynamics, given the high porosity of the ground. It is, however, likely that it has an impact on the groundwater recharge and on runoff further down the valley, with additional implications for the ecological system (e.g. Reference HauerHauer and others, 1997). A network of piezometers together with runoff gauges would be the ideal instrumentation to detect this effect on water levels, and hopefully will be incorporated in future research programs in the region.
Conclusions
A number of photographs were obtained from a fixed camera installed at Volcán Villarrica. These photographs have been used to incorporate daily spatial albedo variations into a glacier-melt model based on AWS data, SRTM elevation data and GPS measurements on the ground. The modelled ablation shows significantly higher values than previously measured on the volcano (Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others, 2007) in areas directly downwind from ash sources, i.e. the crater and exposed tephra banks. These ashes reduce albedo and as a result increase the surface melt rate. The maximum values were found near ash/debris hummocks where winds blow material onto the snow. However, large areas of the glacier experience high ablation (low albedo) in the areas located downstream to the east and southeast of the crater, where airfall tephras are frequently deposited. Another interesting result obtained here is the possibility of detection of crevassing at the flanks of the volcanic cone. Changing patterns of crevasses could indicate subglacial volcanic activity, as well as glacial ice dynamics. Further research will be needed to extract surface velocities using, for example, feature tracking in the daily photographs. The frontal variations of the glacier have also been updated, confirming previous receding trends, in spite of the thick insulating debris which mantles most of the glaciers’ ablation zones. Continued glacier retreat in spite of the insulating tephra is likely to be due to a combination of reduced accumulation due to decreasing precipitation trends in the region (Reference Bown and RiveraBown and Rivera, 2007), rapid melt at exposed ice cliffs along crevasses and possibly enhanced basal melting in areas of high geothermal heating (Reference RiveraRivera and others, 2006). The role of ash deposition on top of the snow surfaces has been discussed here, but more detailed analysis of the mass balance of the glacier needs to be done in order to understand the specific driving factors for the retreating trend, especially how the decrease in precipitation observed in the Chilean lake district (Reference Bown and RiveraBown and Rivera, 2007) is affecting accumulation on the glacier.
Acknowledgements
This work was sponsored by Fondo Nacional de Ciencia y Tecnología, Chile, (FONDECYT 1040515 and 7050177) and Centro de Estudios Científicos (CECS), Chile. CECS is funded in part by the Millennium Science Initiative and the Centers of Excellence Base Financing Program of Comisión Nacional de Investigación Científica y Tecnológica de Chile (CONICYT). Centro de Ingenierrá de la Innovación del CECS is funded by CONICYT and the Gobierno Regional de Los Rios. M. Rodriguez, F. Bown, A. Uribe, G. Gazitúa, Corporación Nacional Forestal, Chile, (CONAF), Servicio Nacional de Geología y Minería, Chile, (SERNAGEOMIN), C. Acuña and R. Koschitzki helped with data collection, preparation of some of the figures, field campaigns and the installation of the equipment. The paper was written while one of the authors (J.G.C.) was funded by a grant from the Austrian Science Foundation (FWF project No. M952-N10). We are grateful to the scientific editor, M.T. Gudmundson, and to L. Nicholson and M. Nolan for helpful comments.