The surface energy balance of an active ice-covered volcano: Villarrica Volcano, southern Chile

Abstract

The energy balance of bare snow and tephra-covered ice near the glacier equilibrium line elevation on Villarrica Volcano, southern Chile, was investigated during 2004 and 2005, combining meteorological, surface temperature and ablation measurements with energy balance modelling. A tephra thermal conductivity of 0.35 Wm^–1 K^–1, and a critical tephra thickness of <5mm at which ablation is reduced compared to bare snow, were obtained from field data. These low values are attributable to the highly porous lapilli particles which make up most of the surface material. Modelled melt totals in the January to March period were 4.95 m and 3.96 m water equivalent (w.e.) in 2004 and 2005, respectively, compared with ∽0.5mw.e. melt for ice buried by >0.1m tephra. Windblown tephra impurities lowered snow albedo, but increased snowmelt by only an estimated 0.28mw.e. over the same period. The net mass balance impact of supraglacial tephra at Villarrica Volcano is therefore positive, as thick ash and lapilli mantle most of the glacier ablation zones, probably reducing annual ablation by several metres w.e. In the accumulation seasons, frequent melting events were recorded with modelled daily snowmelt rates of up to 50 mmw.e.

Introduction

While several workers have examined the effects of (mainly thick) supraglacial rock debris on mountain glacier melt rates (e.g. Nakawo and Young, 1982; Mattson and Gardner, 1989; Conway and Rasmussen, 2000; Takeuchi and others, 2000) the impact of supraglacial tephra deposits on glacier surface energy and mass balance remains poorly researched and is limited to a few short term studies (Lister, 1953; Driedger, 1981; Kirkbride and Dugmore, 2003). The impact of supraglacial tephra on ablation is a balance between melt enhancement due to increased absorption of shortwave radiation, which dominates under thin covers, and melt reduction through insulation from atmospheric heat and insolation, which dominates under thick covers. The thickness at which the sub-debris ablation rate equals the bare snow or ice ablation rate is termed the critical thickness, while the thickness at which sub-debris ablation is maximized is termed the effective thickness. Driedger (1981) identified effective and critical thicknesses of 3 and 24 mm, respectively, for Mt St Helens (1980) tephra overlying snow at South Cascade Glacier, Washington, while Kirkbride and Dugmore (2003) reported respective values of 2 and 5.5 mm for Hekla (1947) tephra overlying ice at Gigjökull, Iceland. The effective and critical thicknesses for tephra are lower than those for rock debris, typically ∽10mm and ∽15–50 mm, respectively (Kirkbride and Dugmore, 2003), due to its relatively low thermal conductivity. Tephra is typically darker than rock debris and small quantities can dramatically reduce snow or ice albedo, increasing ablation even at low concentrations (Driedger, 1981). Consequently, widely diverging surface melt rate responses to relatively small variations in tephra cover can occur if the critical thickness threshold is crossed, resulting in complex spatial and temporal patterns of ablation.

Better understanding of the surface energy balance of glaciers affected by supraglacial tephra deposition is needed for a number of reasons. Melting of snow and ice cover on volcanoes in the summer/dry season makes an important contribution to exploitable freshwater resources in many populated regions, including the Andes, while acceleration of ablation on darkened snow or ice surfaces may contribute to flood hazards during eruptions or extreme climatic events. The extent to which thick supraglacial debris covers can offset glacier retreat due to climatic warming is not known, while climatic interpretation of past glacial advances in volcanic regions from morainic evidence (e.g. Evans and others, 1999) is dependent on an appreciation of the impact of tephra fallout on surface melt.

This paper presents the results from a 2 year investigation of the surface energy balance on the glacier-covered area of Villarrica Volcano in southern Chile, which has been affected by frequent eruptions and tephra deposition throughout the past 500 years (Petit-Breuih and Lobato, 1994; González-Ferrán, 1995; Lara, 2004). The main aims are to quantify the sources of surface melt energy over bare and tephra-covered snow and ice areas and to establish the relationship of the surface melt rate to supraglacial tephra cover thickness. Monitoring of meteorological variables, surface temperature and ablation close to the glacier equilibrium line altitude (ELA) were conducted to enable an assessment of the net impact of the tephra cover on the mass balance of the glacier and its likely response to future climatic warming.

The Study Site

Villarrica Volcano (3982501200 S, 7185602700W, 2847 m a.s.l.) has a glacier-covered area of 30.3 km² (value for 2005; Rivera and others, 2006) mainly distributed toward the south and east where the main glacier basin, Glacier Pichillancahue-Turbio (17.3 km² in 2005), is located (Fig. 1). The ablation zone is covered by an almost continuous ash and lapilli mantle with thickness ranging from a few millimetres to >1 m, except for occasional ice cliffs which are too steep to support significant cover. These surface tephras are composed mainly of andesitic-basaltic material with a silicate content in the range 53–56%, which probably originate as pyroclastic deposits and are advected to the surface by ice flow (Rivera and others, 2006). The volcano is classified as highly active, with over 50 eruptive events documented since 1558 (Petit-Breuih and Lobato, 1994) and permanent degassing, however high magnitude explosive eruptions are infrequent (Lara, 2004). Lahars produced by eruptions of Villarrica Volcano resulted in the deaths of more than 75 people in the 20th century (Stern, 2004).

Fig. 1. Study site: (a) location map; (b) March 2004 Aster image of Villarrica Volcano including an outline of Pichillancahue-Turbio Glacier. Note the extensive tephra cover in the ablation area. Black star = Glacier AWS Site; white star = Ridge AWS Site.

The climate of the Chilean Lakes District (38–418 S) is temperate maritime, characterized by mild winters with high levels of precipitation (April–September precipitation >2m at Valdivia on the coast) and relatively dry, warm and stable summers. The ELA is approximately 2000m (value for 2005; Rivera and others, 2006) and glaciers, many of which are located over active volcanic edifices, have high mass turnover. Most glaciers in the region have been retreating in recent decades in response to atmospheric warming and decreasing precipitation trends since the mid-20th century (Rivera and others, 2002). Rivera and others (2006) calculated an ice thinning rate of 0.81 ±0.45ma^–1 on the tephra-covered area of Villarrica Volcano between 1961 and 2004 and a decrease in glacier-covered area of 0.090±0.034km² a^–1 between 1976 and 2005.

Methods

Field data collection

Meteorological, ablation and surface temperature measurements were made during the 2004 and 2005 ablation seasons at the Glacier AWS (automatic weather station) Site on snow at 1900 m a.s.l. and on nearby tephra banks (Fig. 1). To avoid burial by winter snow accumulation (>5m at this site) the AWS was moved to the exposed Ridge AWS Site at 1890ma.s.l. between April and December 2004 and 2005 (Fig. 1). The AWS recorded fluxes of incoming and reflected shortwave radiation, SWR↓ and SWR↑ respectively (Wm^–2); net allwave radiation AWR* (Wm^–2); air temperature T _a (8C); relative humidity RH (%); vapour pressure e (hPa) and wind speed u (m s^–1) at 2 m (Table 1). The net flux of longwave radiation, LWR ^* (Wm^–2) was calculated by subtracting net shortwave radiation SWR ^* from AWR ^*. The sensors were mounted on cross-arms attached to a tripod which stood freely on the snow or tephra surface at the Glacier AWS Site, and was bolted onto the underlying rock platform at the Ridge AWS Site. Both sites were level and radiation sensors were mounted parallel to the surface. Instrument failure and the difficulty of accessing the site during the winter affected the collection of data from the AWS, but substantial records were obtained for most sensors in both ablation and accumulation seasons (Table 2). A sample interval of 15 s was used and half hour means stored on a Campbell CR10 data logger. The temperature sensor was housed in a naturally ventilated shield (Table 1). Wind speeds were generally sufficient (Table 3) to ensure adequate ventilation for temperature measurement (Georges and Kaser, 2002). Subsequent analysis of meteorological data showed no inverse dependence of T _a on u on any days. Hence, overheating of the temperature sensor due to radiation absorption is likely to be only an occasional problem. Additional point albedo measurements, using a portable albedometer consisting of 2 CM6B pyranometers (Table 1), were taken on snow, tephra and mixed snow/tephra close to the AWS to characterize the albedos of these surfaces.

Table 1. Meteorological instrument specifications. Accuracy indicates the likely error range under field temperature and humidity conditions according to the manufacturer’s specifications (RMS = root mean square error for SWR↓ = 1000Wm^–2)

Table 2. Details of meteorological data obtained

Table 3. Summary of meteorological conditions during 2004 and 2005 at the Glacier AWS Site (ablation season) and the Ridge AWS Site (accumulation season). Daily albedo is the mean of 0900 to 1700 readings

Tephra surface temperature T _s was sampled at a 30 s interval and half hourly means stored on a Campbell CR500 datalogger at 4 sites of 0.025, 0.11, 0.19 and 0.32m thickness on a tephra bank about 100 m east of the Glacier AWS Site in 2004. The 0.32m site was in an exposed position at the top of a ridge, while the other three sites had a more sheltered south-westerly aspect. The probes were placed in the tephra so that the tips were level with the surface, and the cables buried to provide a stable anchor. The main problem was exposure of the thermistors, due to movement or blowing of tephra, which would lead to an underestimate of daytime T _s. Probe positions were maintained during field visits, but checks were not possible between 4 February and 18 April. T _s measurements were repeated at two level sites of 0.09 and 0.13 m thickness on an ash bank approximately 50 m south of the AWS in 2005. Temperatures at the base of the tephra were measured simultaneously with identical probes at each site, and these data demonstrated that the ice surface in contact with the tephra was continuously at melting point throughout the measurement periods.

Snow ablation was recorded at 3 stakes of white, low conductivity material, set in a triangle around the AWS between 17 January and 3 February 2004, and sub-tephra ice ablation recorded at the thermistor sites on the same dates and again on 17 January 2005. Units of surface lowering were converted to mm water equivalent (w.e.) melt using snow density measured in pit (mean density = 585 kgm^–3, range 531–648 kgm^–3) and an ice density of 890 kgm^–3.

The relationship of the snow ablation rate to supraglacial tephra cover thickness was investigated at 10 plots close to the Glacier AWS Site in January 2005. Seven plots were artificially covered with tephra of 5, 10, 20, 30, 40, 50 and 90mm thickness obtained from a nearby tephra bank. An eighth plot was sprinkled with a discontinuous tephra cover, with an estimated loading of 0.1 kgm^–2. These artificial tephra plots were created as similar as possible to the loosely compacted and well-sorted ‘natural’ tephra cover. Daily surface lowering rates were recorded at stakes set in the centre of each plot and compared with daily melt rates recorded on ‘natural’ snow, which incorporated trace amounts of tephra, dust and other airborne impurities, and ‘clean’ snow which had all visible impurities artificially removed.

Energy balance modelling: snow

Energy balance modelling of ice and snow surfaces using meteorological data is well established, and the principles are explained elsewhere (e.g. Greuell and Genthon, 2004). Consequently, only important details will be described here. In this study, the energy balance model of Brock and Arnold (2000) was used. Measured inputs of SWR↓, SWR↑, AWR ^*, T _a, e and u were used to calculate the surface fluxes (in Wm^–2) of net shortwave radiation SWR ^*, net longwave radiation LWR ^*, the turbulent fluxes of sensible and latent heat, SHF and LHF respectively, and the surface melt rate MELT assuming the surface was at melting point:

(1)

The latter assumption is acceptable during summer melt periods when T _a is generally >08C, as the conductive heat flux into snow is very small. During the ablation season, heat advected by precipitation is small and was consequently ignored. Fluxes are positive when directed towards the surface. The model is run using hourly-averaged values of meteorological variables. The turbulent fluxes are calculated from the ‘bulk’ aerodynamic method using measurements of u, T _a and e at 2m, and assuming the respective surface values are 0ms^–1, 08C and saturation vapour pressure at 08C. Atmospheric stability in the surface layer is accounted for using the Monin-Obukhov similarity equations. An aerodynamic roughness length z ₀ of 1 mm was used based on published values for melting snow on mid and low latitude glaciers (Brock and others, 2006). The Brock and Arnold (2000) model has been shown to calculate surface melt rates of comparable accuracy to stake and ultrasonic ablation measurements when local meteorological data are used (Pellicciotti and others, 2005).

Due to the absence of wind speed measurements during the accumulation seasons (Table 2) an enhanced temperature index melt model which requires only T _a and SWR ^* as input data (Pellicciotti and others, 2005) was used. The snow melt rate, in units of mmw.e., is calculated from hourly data as follows:

(2)

In a comparison at 5 automatic weather station sites on Haut Glacier d’Arolla, Switzerland, this model explained ∽90% of the melt rate variations calculated by the Brock and Arnold (2000) model. The model parameters were recali-brated for Villarrica Volcano using melt output from the energy balance model for periods when all meteorological variables were available. The SWR ^* multiplier is almost identical to the original value (0.0094), but the temperature multiplier is higher than in the original study (0.05), probably due to higher wind speed and/or higher humidity and cloud cover at Villarrica Volcano.

Energy balance modelling: tephra

The balance of energy inputs and outputs at a tephra surface can be expressed as:

(3)

where COND is the conductive heat flux in the tephra. All fluxes are positive when directed towards the surface and calculated in units of Wm^–2. Incoming SWR and LWR (after subtraction of 316Wm^–2, the assumed outgoing flux for a snow surface at 08C, from LWR ^*) were measured directly at the AWS. SWR↑ was calculated using an albedo of 0.05, the mean measured tephra albedo, while outgoing LWR was calculated from T _s using the Stefan-Boltzmann relationship, assuming the tephra surface emits as a black body.

Assuming a uniform layer with a linear temperature profile with depth, heat conduction in tephra was calculated from measured T _s values as:

(4)

where K is thermal conductivity (Wm^–K^–1) and d is tephra thickness (m). Since heat conduction into the ice can be ignored, and the net change in heat stored in a debris layer over long time periods (>1 week) is ∽0 (Nakawo and Young, 1982), the sub-tephra ice melt rate (mmw.e. s^–1) can be calculated as:

(5)

where L ^f is the latent heat of fusion of water (334 000 J kg^–1).

The tephra surface was completely dry except during infrequent rainfall events, so LHF was assumed to be zero at the surface. Hence, SHF can be found as a residual since all other terms in Equation (3) are calculated. Previously, SHF has been calculated above supraglacial debris using the bulk aerodynamic approach:

(6)

where D is a bulk exchange/transfer coefficient ( Jm^–3 K^–1). Nakawo and Young (1982), Mattson and Gardner (1989) and Kayastha and others (2000) used a value for D of 4.9, obtained from studies of sensible heat transfer over snow surfaces (Naruse and others, 1970), while Takeuchi and others (2000) calculated D assuming a logarithmic (neutral atmosphere) wind speed profile and measured z ₀. These approaches neglect variation in atmospheric stability in the surface layer, which will become moderately to strongly unstable as T _s rises well above T _a during the daytime, leading to an underestimation of SHF as instability enhances convection.

The bulk aerodynamic approach can be extended to non-neutral conditions using the Richardson number (Oke, 1987). Following Deardorff (1968), the bulk form of the Richardson number R _b is defined by:

(7)

where g is acceleration due to gravity (9.8m^–1 s ^–2) and z _a is measurement height (m). D is evaluated as follows, based upon the classification of the surface layer:

Stable case, R _b positive:

Unstable case, R _b negative:

(8)

where k is von Karman’s constant (0.40, Oke, 1987), ρ _a is air density (assumed to be 1.26 kgm^–3 for low positive air temperatures) and c _a is the specific heat of air (1010 J kg^–1 K^–1, Oke, 1987). This approach has been applied successfully to snow and ice surfaces (e.g. Anderson, 1976; Moore and Owens, 1984; Marcus and others, 1985). Despite its ability to account for both stable and unstable atmospheric regimes, the R _b approach to calculating turbulent fluxes has not previously been applied to debris-covered glaciers. To test the validity of the different approaches, SHF values were determined using the simple bulk aerodynamic calculation without stability correction (Equation (6)), with D = 4.9 Jm^–3 K^–1, and using stability functions based on R _b (Equations (7) and (8)) and compared with the ‘residual’ SHF values obtained from Equation (3). Initially, a z ₀ of 1mm was used in Equation (8). Subsequently, a sensitivity test was performed using z ₀ values of 0.5, 2, 3 and 5 mm.

Calculation of tephra thermal conductivity

Following Drewry (1972), assuming a purely conductive system, thermal conductivity K can be found from:

(9)

where COND is heat flux per unit area averaged over a measurement period, assumed equal to the sub-tephra MELT flux averaged over time, and is the average thermal gradient through the medium, normal to the surface. Temperature and sub-tephra ablation (melt energy) were obtained at the four-instrumented tephra sites between 17 January and 3 February 2004. The net change in heat energy stored in tephra, әSTOR (Wm^–2) was also calculated, following Mattson and Gardner, 1989:

(10)

where ρ _t and c _t are the tephra density (250 kgm^–3) and

specific heat (804 J kg^–1 K^–1; Clark, 1966), respectively, and is the average rate of temperature change (K s^–1) where T̄ _t is defined as T _s2.

Results

Snow energy balance model validation and sensitivity analysis

Modelled snow ablation closely matched measured ablation between 17 January and 24 January 2004, but exceeded measured ablation on the final stake measurement on 3 February (Fig. 2). The largest potential error source in the model is LWR ^* (Table 1) which varies the modelled cumulative melt by ±14%, large enough to explain the difference between modelled and measured melt on 3 February (Fig. 2). Varying z ₀ from 0.1–10.0 mm, a larger z ₀ range than previously recorded on melting snow on mid-latitude glaciers (Brock and others, 2006), varies the modelled cumulative melt by ±8%, but is not large enough to explain the modelled melt overestimation. Possible overheating of the temperature sensor at times of low wind speed and high insolation is a small melt error source (1–2%). Some melt overestimation is likely due to the model’s assumption that T _s never falls below <08C. Ablation stake measurements are sensitive to the snow density value used, and modelled cumulative melt matches measured melt for the maximum snow density value of 0.648 kgm^–3 (Fig. 2). Furthermore, it is well established that surface lowering recorded at ablation stakes lags ablation under conditions of strong solar radiation and surface longwave radiative losses, as experienced on the days prior to 3 February due to the formation of a surface ‘weathering crust’ (Müller and Keeler, 1969; Munro, 1990). In contrast, the days prior to the 24 January ablation stake reading, which closely matched modelled melt, were dominated by high wind speeds in excess of 10ms^–1, which probably removed the ‘weathering crust’. Potential error sources in the ablation stake readings therefore make further evaluation of the performance of the melt model difficult.

Fig. 2. Comparison of modelled snow ablation (solid line) and ablation measured at 3 ablation stakes for minimum, mean and maximum snow density values (black squares with bars) at the Glacier AWS Site. Modelled melt sensitivity to measurement errors and specified z ₀ value is shown as: outer dashed lines = variation in LWR* by ±30Wm^–2; inner dashed lines = variation in z ₀ from 0.1–10 mm; broken line below solid line = reduction in air temperature by 18C when SWR↓>200Wm^–2.

Meteorological conditions, snow energy balance and melt: ablation season

Meteorological conditions were characterized by almost continuously positive T _a, high SWR↓ and moderate to high u (Table 3). Although the periods for which data are available are not identical in each year (Table 2), SWR↓, u and e were very similar, but 2004 was >18C warmer than 2005. Correspondingly, the modelled mean daily ablation season snowmelt rate was higher in 2004 (55mmw.e. d^–1) than in 2005 (44mmw.e. d^–1), due to smaller surface longwave radiative losses and larger gains from the turbulent fluxes (Fig. 3). Modelled melt totals over the 3 month January– March period were 4.95 and 3.96 mw.e. in 2004 and 2005, respectively. SWR ^* dominated surface energy receipts, contributing 67% of total snowmelt energy in 2004 and 73% in 2005.

Fig. 3. Modelled mean daily snow surface energy fluxes at the Glacier AWS Site during the January to March periods in 2004 and 2005. SWR* = net shortwave radiation flux; LWR* = net longwave radiation flux; SHF = sensible heat flux; LHF = latent heat flux. All fluxes were converted fromWm^–2 to units of mmw.e. melt using L _f.

The Pearson’s correlation coefficient between daily albedo and daily snowmelt was –0.43 in 2004 and –0.47 in 2005, with snowfall events such as those on 5 March and 3 April 2005 having a dramatic impact on melt rates (Fig. 4). Snow albedo decreased gradually from >0.6 in January to <0.2 in mid March, lower than minimum snow albedo of ∽0.4, typical for alpine glaciers (Paterson, 1994; Brock and others, 2000), due to the increasing concentration of tephra in the snowpack. In order to estimate the increase in snowmelt caused by snow darkening by tephra impurities, SWR ^* was recalculated using an albedo parameterization for an alpine glacier (Brock and others, 2000) in place of the measured albedo. This parameterization, which was developed using over 300 measurements of albedo over three ablation seasons on an alpine glacier, calculates daily snow albedo as a function of accumulated daily maximum temperatures (ADT) since the last significant snowfall:

(11)

Fig. 4. Overview of modelled daily melt totals and mean daily albedo at the Glacier AWS Site.

Although the parameterization is not truly physical, ADT correlates strongly with the physical causes of snow albedo variation on an alpine glacier, namely snow metamorphism, melting and incorporation of impurities such as organic matter and small particles of rock and soot, but not tephra (Brock and others, 2000). Based on Equation (11), tephra impurities cause a modest increase in total ablation season snowmelt of 0.28mw.e. relative to an alpine glacier under the same meteorological conditions.

Meteorological conditions and snow ablation: accumulation season

The variable but often mild mid to late winter conditions in 2004 contrasted with the less variable and colder early to mid winter conditions in 2005 (Fig. 5, Tables 2 and 3). August and early September 2004 were characterized by four melt events with mean daily temperature >58C and daily melt totals >20mmw.e. The melt event of 8–11 August 2004 is notable as positive air temperatures and high melt rates were sustained for 4 days, with a maximum air temperature of 14.68C (night time minimum 10.18C) and minimum relative humidity of 3% on 9 August, producing an estimated 160 mmw.e. melt. This event was triggered by an eastward shift of the southern Pacific high pressure zone into Argentina, which introduced a north-easterly airflow and the incursion of a tropical air mass into southern Chile. Several melt events also occurred throughout the 2005 winter season. Although these events were of smaller magnitude and shorter duration than in 2004, they demonstrate that winter ablation is an important characteristic of the glacier surface energy balance on Villarrica Volcano (Fig. 5).

Fig. 5. Daily melt totals (columns) and mean daily temperature (lines) during the (a) 2004 and (b) 2005 accumulation seasons at the Ridge AWS Site.

Tephra conductivity, K

The mean K for the four measurement sites was 0.35Wm^–1 K^–1, with a range of 0.23 to 0.46 (Table 4), comparable to the K of pumice (Clark, 1966). At Sites 3 and 4, however, there were large changes in tephra thickness over the measurement period, probably due to melting out of englacial tephra and the addition of material from upslope. This leads to uncertainty in the calculated K value at these sites, given the sensitivity of K to errors in thickness measurements (Table 4). The K of tephra obtained is between one half and one sixth of published K values for rock debris overlying glaciers (0.85–2.6Wm^–1 K^–1; Nakawo and Young, 1982; Conway and Rasmussen, 2000). Tephra therefore provides much more effective insulation of ice from solar radiation and atmospheric heat than rock debris, for an equivalent thickness.

Table 4. Tephra conductivityK. Results from measurements at 4 sites between 17 January and 3 February 2004

Tephra surface temperature showed large diurnal fluctuations, responding to the cycle of SWR ^* during the daytime, while night time minima were controlled by T ^a, falling to 08C when T ^a was <58C (Fig. 6). Hence, the build up of stored heat energy in the tephra during the daytime was quickly released to the atmosphere during the night. Over the period 17 January to 3 February, ∂STOR was between –0.011 and 0.113 MJm^–2 for the four monitored sites, equivalent to the energy required to melt just 0.03–0.3mm of ice, confirming the assumption that ∂STOR is ∽0 over periods >1 week.

Fig. 6. Air temperature (solid line, upper trace) and tephra surface temperatures at four sites (broken lines, lower traces) over the tephra conductivity calculation period.

Relationship of buried snow and ice ablation to supraglacial tephra thickness

Snow ablation rates increased by 97% relative to a ‘natural’ bare snow surface, when the surface was artificially sprinkled with a discontinuous tephra cover (estimated loading of 0.1 kgm^–2; Fig. 7), which is consistent with a reduction in snow albedo from ∽0.55 to ∽0.1. The percentage increase in snowmelt rate for light tephra or dust loading corresponds with previous results (Driedger, 1981; Adhikary and others, 2002). Thereafter, for increasing tephra thickness, ablation is reduced relative to bare snow. The relative ablation reduction is –33% for a 5 mm tephra cover and –90% for a 50mm tephra cover, with ablation rates just 5% of ‘natural’ snow estimated for d > 100 mm. As discontinuous tephra does not have a measurable ‘thickness’ as such, no critical thickness value for Villarrica tephra was identified. Rather, as soon as tephra forms a continuous cover, ablation is reduced relative to bare snow. This result is in contrast to previous studies. Driedger (1981) identified a critical thickness of 24mm for Mt St Helens tephra, a thickness at which Villarrica tephra reduces ablation by ∽75% (Fig. 7). Similarly, Lister (1953) and Kirkbride and Dugmore (2003) identified critical thicknesses of ∽5mm for Icelandic tephras. Only the critical thickness of 1.33 mm identified from artificial loading of glacier ice and seasonal snow patches with dust (largely organic matter) of Adhikary and others (2000) is similar.

Fig. 7. Relationship between tephra thickness and ablation relative to a tephra-free snow surface.

Measured total annual buried ice ablation for the four tephra sites with mean depths 0.11–0.32m were between 0.42–0.56mw.e., just 8–14% of the modelled bare snowmelt in the January–March period. Modelled mean daily ablation season (January to March) melt rates for these four sites, of 3–8mmw.e., are an order of magnitude lower than for comparable thicknesses of supraglacial rock debris (Mattson and Gardner, 1989; Takeuchi and others, 2000). A daily melt rate of ∽1mmw.e. is extrapolated beneath tephra 1m thick. Tephra thickness increased significantly over time at all sites, suggesting a long-term decrease in buried ice ablation rates.

Surface energy balance of tephra and buried ice melt rates

The ablation season energy balance at tephra surfaces was dominated by SWR ^*, which was the only net energy input (Fig. 8). Most of the incoming solar energy is converted to outputs of LWR and SHF, due to high daytime T _s, which was regularly >308C and occasionally >408C (Fig. 6). LWR ^* and SHF were of similar magnitude, but slightly smaller in 2004, due to lower T _s at the 3 sites with south-westerly aspect, and higher T _a (Table 3). COND is relatively small compared to the other fluxes. The mean daily fluxes of 4.4 mm (mean d = 0.20 m) and 6.8 mm (mean d = 0.11 m) w.e. melt in the 2004 and 2005 ablation seasons, respectively, would melt totals of 0.40 and 0.61 mw.e. ice over the January to March period. The fact that measured annual ice melt beneath 0.11–0.32m tephra was 0.42–0.56mw.e. implies that COND drops to negligible values outside of this 3 month ablation period, due to snow cover and lower T _a and SWR↓.

Fig. 8. Modelled mean daily surface energy fluxes at tephra measurement sites, January–March periods in 2004 and 2005. SWR* = net shortwave radiation flux; COND = conductive heat flux in tephra; LWR* = net longwave radiation flux; SHF = ‘residual’ sensible heat flux. Mean tephra thickness: 0.2m in 2004 (4 sites); 0.11 m in 2005 (2 sites). All fluxes were converted from Wm^–2 to units of mmw.e. melt using L _f.

Up to 2 ^–2 of heat energy can be stored in the tephra layer during the daytime if T̄ _t rises to >208C. This energy is returned to the atmosphere as longwave radiation and sensible heat during the late afternoon and night and ∂STOR is <±0.1^–2 on a daily time step. Consequently, ∂STOR is negligible in melt calculations over periods >1 week.

The simple bulk aerodynamic formula, assuming a neutral atmosphere (Equation (6)), severely underestimated the ‘residual’ SHF at the tephra surface by 55% in 2004 and 37% in 2005 (Fig. 9a). SHF calculations which accounted for atmospheric instability using R _b (Equations (7) and (8)) equalled almost exactly the ‘residual’ SHF in 2005, but underestimated it by 25% in 2004. SHF calculations using R _b are dependent on the z ₀ value specified (Fig 9b), matching the ‘residual’ SHF for a z ₀ value of 3mm in 2004 and 1 mm in 2005. The larger z ₀ value inferred in 2004 corresponds with the visibly larger size of roughness elements on the tephra bank compared with 2005. The ‘residual’ SHF cannot be assumed to be the true reference value, as the calculation is prone to errors in the measurement of the other surface energy fluxes. In particular, the assumption of a tephra emissivity of 1.0 probably led to an overestimate of the outgoing LWR flux and an underestimate of the ‘residual’ SHF. Overall, however, these results support the application of instability functions based on R _b for the calculation of turbulent fluxes over debris-covered glacier surfaces, as long as a representative z ₀ value is used. SHF calculations made assuming a neutral atmosphere severely underestimate the true flux, due to strong atmospheric instability when daytime T _s is high.

Fig. 9. (a) Comparison of calculated and residual SHF values in the January to March periods in 2004 and 2005. ‘Neutral’ = simple bulk aerodynamic formula assuming a neutral atmosphere (Equation (6)), ‘Richardson’ = bulk aerodynamic formula using R _b to account for atmospheric instability (Equations (7) and (8)). (b) Dependence of the R _b SHF calculation on z ₀ in 2004 (dashed line) and 2005 (solid line). The diamond and square mark the residual SHF values in 2004 and 2005, respectively.

Discussion

Impact of supraglacial tephra cover on glacier mass balance

The low albedo and conductivity K of tephra create diverging melt responses, dependent upon whether the supraglacial tephra cover is continuous or not. Aeolian redistribution of ash from tephra banks and fallout of fine particles from fumarolic activity in the volcano darkens snow surfaces, particularly in the latter stages of snowmelt. The resulting increase in snowmelt (decimetres w.e.) is, however, an order of magnitude lower than the reduction in ice melt under a continuous tephra cover a few centimetres thick (metres w.e.). Possible differences in the supply of non-volcanic impurities to the snow between the alpine site where the albedo parameterization (Equation (11)) was developed and Villarrica Volcano adds some uncertainty to the estimation of snowmelt enhancement due to tephra impurities, but this is unlikely to affect total modelled snowmelt by more than a few centimetres w.e. Considering that continuous tephra mantles most of the ablation area (Fig. 1) the net mass balance impact of supraglacial tephra on Villarrica Volcano is positive, probably reducing negative specific mass balances in the ablation zone by a few metres w.e. This interpretation, however, ignores the roles of ice cliffs in the ablation zone, which are too steep to support continuous tephra cover. Ablation at ice cliffs is extremely rapid due to the low albedo of ash-darkened slopes and it is possible these features account for more ablation than stable tephra-covered areas, despite being of limited spatial extent. Further research is needed to assess the contribution of these features to total ice mass loss.

A linear vertical temperature gradient through tephra was assumed in both the calculation of ice ablation beneath a tephra layer (Equations (4) and (5)) and in the calculation of tephra K (Equation (9)) as temperature was measured only at the surface and base of the tephra layer. Previous work, e.g. Mattson and Gardner (1989), has demonstrated that temperature profiles within supraglacial debris deviate from the linear, in particular the gradient becomes steeper towards the surface during the daytime. In the absence of internal tephra temperature measurements it is impossible to assess the impact of the linear temperature profile assumption on melt and K calculations, but it is likely that over- or underestimation of the temperature in different parts of the profile will tend to compensate over the periods of >1 week for which calculations were made.

The mean K value of 0.35 Wm^–1 K^–1 for Villarrica tephra is corroborated by the low buried ice melt rates recorded and the published K for pumice (e.g. Clark, 1966). The greater insulating properties of Villarrica tephra compared with tephras reported in earlier studies may be due to particle size, which is a function of airborne transport distances. In the studies of Driedger (1981) and Kirkbride and Dugmore (2003), supraglacial tephra had travelled several 10s of kilometres by air from the eruption sites and consisted almost entirely of ash <1.0mm in diameter. In contrast, supraglacial tephra on Villarrica Volcano consists of local fallout <2 km from source, and is composed mainly of lapilli >2.0mm in diameter. The lapilli particles are highly porous and it is probable that, as a surface layer, they trap a greater volume of air which has a lower thermal conductivity than an equivalent layer of ash particles. The implication is that tephra covers originating from local fallout have a thinner critical thickness than tephra covers resulting from distal fallout (cf. Kirkbride and Dugmore, 2003).

The annual ablation at four tephra-covered sites accounts for only 52–69% of the mean 1961–2004 ablation zone thinning rate of 0.81 ma^–1 calculated by Rivera and others (2006) from comparison of airborne synthetic aperture radar and GPS datasets with earlier topographic maps. The remaining thinning is probably caused by subglacial melting over areas of high geothermal flux (Rivera and others, 2006).

Implications for melt modelling of debris-covered glaciers

In common with most previous studies, we used measured T _s and d to calculate the conductive heat flux into supraglacial debris (tephra). For distributed melt modelling, or melt calculations for remote debris-covered glaciers, however, such measurements are impractical. Thermal satellite measurements are unable to provide the necessary data as their temporal resolution is insufficient considering the large daily variation in T _s. It is unlikely that the spatial distribution of debris thickness across a glacier, its thermal properties and moisture content, nor their variation in time, will be known with sufficient accuracy to enable a physical-based calculation of T _s based on meteorological observations alone. Alternatively, T _s could be estimated from an empirical, or semi-empirical, relationship to meteorological data.

Comparisons with ‘residual’ SHF values have shown that from knowledge of T _s, and T _a and u at one height, the SHF at a tephra surface can be modelled successfully by accounting for changing stability in the atmospheric surface layer using R _b. For this method to be applied more widely, measurements of z ₀ over tephra and other supraglacial debris are needed, given the dependence of SHF on surface roughness. Modelling of latent heat fluxes is more problematic since the surface humidity value is not known. While the latent heat flux was ignored as the tephra surface was almost continuously dry, evaporation will be an important flux following rainfall and also within tephra layers, due to the presence of meltwater. In order to calculate latent heat fluxes on tephra, or other types of supraglacial debris, measurement of air vapour pressure within debris layers is needed.

The tephra conductivity K is assumed constant with depth. This assumption is unlikely to be true as K probably increases with depth as low conductivity air spaces are filled with higher conductivity (melt) water. The method applied, however, provides an in situ value of K under field conditions, which can be used in energy and mass balance models.

Response to climatic change

By temporarily storing energy from solar radiation and returning it to the atmosphere as longwave radiation and sensible heat, tephra cover provides effective insulation for ice and snow. Consequently, tephra-covered glaciers are less sensitive to atmospheric temperature changes than ‘clean’ glaciers. In the ablation season, the main effect of future climate warming will be to accelerate the melting of the winter snow cover, and increase melt on any bare ice surfaces exposed above or below the zone of continuous tephra cover. Climate warming is likely to have a greater impact on the winter mass balance. An amount of 0.37mw.e. melt was calculated in July and August 2004 at 1890m a.s.l., and in several of the winter melt events it is probable that melting extended to the top of the accumulation zone. While much of this surface melt would refreeze within the snowpack, continued warming is likely to increase the frequency and magnitude of winter melt events. The net effect will be to reduce snowfall accumulation, reinforcing the lower winter balances caused by decreasing precipitation trends over the last 50 years (Bown and Rivera, 2007).

Conclusions

This study has investigated the energy balance of bare snow and tephra-covered ice on an active volcano in southern Chile over a two year period, combining meteorological, surface temperature and ablation measurements with energy balance modelling.

The mean K of in situ tephra calculated was 0.35Wm^–1 K^–1, much lower than published K values for rock debris and the inferred K for ash tephra at other sites. The critical thickness of Villarrica tephra is <5mm and melt rates of ice buried beneath >80mm tephra are <10mmw.e. d^–1, an order of magnitude lower than for an equivalent thickness of rock debris. It is likely that the highly porous lapilli traps a larger volume of air compared with equivalent thicknesses of ash tephra, reducing K. One consequence is that larger sized tephra resulting from local fallout could have a more positive impact on mass balance than fine airborne tephra resulting from distal fallout.

At the ELA, shortwave radiation provided 67–73% of total ablation season snowmelt energy and summer air temperatures were almost continuously positive. Mean daily snowmelt rates were 55 and 44 mmw.e. in the 2004 and 2005 ablation seasons, respectively. Total modelled melt in the January to March period was 4.95 and 3.96mw.e. in 2004 and 2005, respectively, compared with ∽0.5mw.e. for ice buried by >0.1m tephra. Small quantities of tephra impurities increase annual snowmelt by about 8%, compared with a melt reduction of >90% beneath tephra >0.1m thick. Considering the extensive cover of thick tephra in the ablation zone, the net effect on glacier mass balance is positive, probably reducing mass loss by several metres w.e. in the ablation zone. Winter melt events occurred frequently in both accumulation seasons, with snowmelt rates of up to 50 mmw.e. d^–1. Consequently, winter accumulation and the existence of glaciers in this region are vulnerable to rising temperatures, despite the stabilizing influence of tephra in ablation zones.

Future work should seek to measure tephra thermal properties at other sites and develop methods for calculating tephra surface temperature from meteorological data for melt modelling of ice-covered volcanoes. The contribution of back-wasting of bare ice cliffs to total mass loss in the ablation zone should also be assessed. Further monitoring of meteorological conditions on glaciers in the region is needed given the potential significance of winter melting events to glacier mass balance. While turbulent sensible heat fluxes can be successfully modelled using a correction for unstable atmospheric conditions, measurement of humidity at the surface of, and within, tephra layers is needed to calculate latent heat fluxes. Measurement of z ₀ over tephra and other supraglacial debris is also required for the calculation of turbulent fluxes.