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 (Reference Georges and KaserGeorges 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.

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. Reference Greuell, Genthon, Bamber and PayneGreuell and Genthon, 2004). Consequently, only important details will be described here. In this study, the energy balance model of Reference Brock and ArnoldBrock 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:

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 (Reference Brock, Willis and SharpBrock and others, 2006). The Reference Brock, Willis and SharpBrock 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 (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti 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 (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005) was used. The snow melt rate, in units of mmw.e., is calculated from hourly data as follows:

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 Reference Brock and ArnoldBrock 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:

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:

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 (Reference Nakawo and YoungNakawo and Young, 1982), the sub-tephra ice melt rate (mmw.e. s^–1) can be calculated as:

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:

where D is a bulk exchange/transfer coefficient ( Jm^–3 K^–1). Reference Nakawo and YoungNakawo and Young (1982), Reference Mattson and GardnerMattson and Gardner (1989) and Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others (2000) used a value for D of 4.9, obtained from studies of sensible heat transfer over snow surfaces (Reference Naruse, Oura and KojimaNaruse and others, 1970), while Reference Takeuchi, Kayastha and NakawoTakeuchi 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 (Reference OkeOke, 1987). Following Reference DeardorffDeardorff (1968), the bulk form of the Richardson number R _b is defined by:

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:

where k is von Karman’s constant (0.40, Reference OkeOke, 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, Reference OkeOke, 1987). This approach has been applied successfully to snow and ice surfaces (e.g. Reference AndersonAnderson, 1976; Reference Moore and OwensMoore and Owens, 1984; Reference Marcus, Moore and OwensMarcus 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 Reference DrewryDrewry (1972), assuming a purely conductive system, thermal conductivity K can be found from:

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 Reference Mattson and GardnerMattson and Gardner, 1989:

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

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

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 (Reference Brock, Willis and SharpBrock 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’ (Reference Müller and KeelerMüller and Keeler, 1969; Reference MunroMunro, 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 (Reference PatersonPaterson, 1994; Reference Brock and ArnoldBrock 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 (Reference Brock, Willis and SharpBrock 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:

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 (Reference Brock and ArnoldBrock 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 (Reference Clark and ClarkClark, 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; Reference Nakawo and YoungNakawo and Young, 1982; Reference Conway and RasmussenConway 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.

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 (Reference DriedgerDriedger, 1981; Reference Adhikary, Yamaguchi and OgawaAdhikary 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. Reference DriedgerDriedger (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, Reference ListerLister (1953) and Reference Kirkbride and DugmoreKirkbride 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 Reference Adhikary, Nakawo, Seko and ShakyaAdhikary 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 (Reference Mattson and GardnerMattson and Gardner, 1989; Reference Takeuchi, Kayastha and NakawoTakeuchi 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.

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. Reference Mattson and GardnerMattson 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. Reference Clark and ClarkClark, 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 Reference DriedgerDriedger (1981) and Reference Kirkbride and DugmoreKirkbride 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. Reference Kirkbride and DugmoreKirkbride 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 Reference RiveraRivera 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 (Reference RiveraRivera 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 (Reference Bown and RiveraBown and Rivera, 2007).