2.1. Theory: turbulent flux measurements over glacier surfaces

Calculations of turbulent sensible and latent heat fluxes between a glacier surface and the air above it are commonly made using the ‘bulk’ aerodynamic method, which assumes airflow in the SABL is turbulent and fully adjusted to the underlying terrain (e.g. Reference MunroMunro, 1990; Reference Ishikawa, Owens and SturmanIshikawa and others, 1992; Van de Reference Van de Wal, Oerlemans and van der HageWal and others, 1992; Reference Hock and HolmgrenHock and Holmgren, 1996; Reference Hock and NoetzliHock and Noetzli, 1997; Reference Brock, Willis, Sharp and ArnoldBrock and others, 2000). Provided the influence of atmospheric stability is accounted for, this method is the most appropriate on sloping glacier surfaces where the wind-speed maxima are within a few metres of the surface (Reference Denby and GreuellDenby and Greuell, 2000). Its principal advantage is that measurements of horizontal wind speed, temperature and humidity need only be made at one height (usually 1 or 2 m) above the surface, as long as the z ₀ of the glacier surface in question is known. The value of z ₀ can be combined with a surface-renewal model (Brutseart, 1975; Reference AndreasAndreas, 1987) to determine the roughness lengths of temperature and humidity, also required in the flux calculations (Reference Denby and GreuellDenby and Greuell, 2000; Reference Denby and SnellenDenby and Snellen, 2002).

The accuracy of the bulk method is dependent on the accuracy with which z ₀ can be specified. An order-of- magnitude increase in z ₀ will more than double the value of the turbulent fluxes (Reference Brock, Willis, Sharp and ArnoldBrock and others, 2000) and an error in z ₀ of this magnitude is more significant to the turbulent flux calculation than neglect of atmospheric stability (Braith- waite, 1995). Aerodynamic roughness values recorded over melting glacier surfaces vary over three orders of magnitude, in the 0.1–10mm range (Table 1). At high latitudes, the recorded z ₀ range is five orders of magnitude from 0.001 to 10 mm (Table 2). Often z ₀ is not measured in glacier energy- balance studies and, due to the lack of a suitable parameterization scheme, a published value from another study, which may not necessarily be appropriate, must be used instead.

Table 1. Published aerodynamic roughness lengths recorded over mid- and low-latitude glaciers. The measurement method is indicated by letter as follows: e – eddy covariance; m – microtopographic; p – wind profile; r – residual in closed energy balance. Where available, the 1 standard deviation range is given in parentheses after the mean z ₀ value

2.2. Controls on z ₀ variation on glaciers

Under the normally turbulent flow conditions over melting snow and ice surfaces (Reference AndreasAndreas, 1987), z ₀ depends solely on the dimensions, form and density distribution of surface roughness elements (Reference OkeOke, 1987; Reference StullStull, 1988). The value of z ₀ increases with increasing height, surface area and density of surface roughness elements, until the ratio of the silhouette area (upwind face of elements) to unit ground area covered by each element reaches 0.4, when a transition to ‘skimming’ flow occurs (Reference OkeOke, 1987; Reference GarrattGarratt, 1992) and z ₀ begins to decrease.

Over mid-latitude glaciers, z ₀ values recorded over smooth fresh snow surfaces are at the 0.1 mm scale, but lower values at the 0.01, or even 0.001 mm, scale have been recorded on snow over polar glaciers and ice sheets (Table 2). Values reported for melting snow surfaces are in the 1–10mm range, due to the development of ablation hollows and other microtopographical features in the snow surface, with extremely high z ₀ values inferred for snow penitentes (Table 1).

Table 2. Published aerodynamic roughness lengths recorded over high-latitude glaciers and ice sheets. The measurement method is indicated by letter as follows: e – eddy covariance; m – microtopographic; p – wind profile; r – residual in closed energy balance. Where available, the 1 standard deviation range is given in parentheses after the mean z ₀ value

On melting ice surfaces, dirt cones and boulder tables caused by surface insulation, and cryoconites and other features resulting from local melt differentials, create small- scale morphological features. Ice strain also creates surface morphology (e.g. crevasses and longitudinal foliae) which may enlarge through enhanced ablation along darker albedo bands. Correspondingly, while z ₀ values recorded over smooth ice are at the 0.1 mm scale, the majority of z ₀ values recorded on melting glacier ice are in the 1–10mm scale range (Tables 1 and 2). The very large z ₀ values recorded by Reference Smeets, Duynkerke and VugtsSmeets and others (1999) and Reference ObleitnerObleitner (2000) relate to an area with very large roughness elements of 1–2 m height in the ablation zone of Breidamerkurjökull, Iceland (Table 1), while extremely low z ₀ values are reported over Antarctic blue ice (Table 2).

Reference Smeets, Duynkerke and VugtsSmeets and others (1999) observed z ₀ to increase from a few millimetres to several tens of millimetres over the ablation season at Breidamerkurjökull in response to the growth of ice roughness elements from the 0.1 to 1 m scale. Similarly, Reference Arnold and ReesArnold and Rees (2003) recorded an increase in snow z ₀ from 0.04–0.05 mm to 0.2–0.3 mm between spring and mid-summer at midre Lovénbreen on Svalbard, with development of ablation hollows in the snow surface. In the ablation zone of the Greenland ice sheet, Reference Grainger and ListerGrainger and Lister (1966) observed z ₀ to decrease from 11 to 6.8 mm, then to 5.8 mm, with changes in surface material from coarse snow sastrugi to melting snow and finally to rough ice. In contrast, Reference Denby and GreuellDenby and Smeets (2000) and Reference Greuell and SmeetsGreuell and Smeets (2001) recorded no variation in ice z ₀ over several months of measurements at Breidamerkurjökull and at Pasterzenkees, Austria, respectively, which corresponded with no visible changes to the roughness of the surface. Similarly, Reference Grainger and ListerGrainger and Lister (1966) reported no significant change in ice z ₀ in the lower ablation zone of the Greenland ice sheet over an ablation season. Overall, therefore, it is unclear whether there is a typical pattern of z ₀ evolution over glaciers during the ablation season.

The value of z ₀ varies with wind direction over irregularly shaped obstacles (e.g. sastrugi) (Reference Jackson and CarrollJackson and Carroll, 1978; Reference InoueInoue, 1989; Reference King and AndersonKing and Anderson, 1994). For many glaciers wind direction is dominated by katabatic flows and topographic control, and is fairly constant (Reference Greuell, Knap and SmeetsGreuell and others, 1997; Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others, 2004). Therefore, dependence of z ₀ on wind direction may not be of great significance to turbulent flux calculations in most cases.

Only a few of the studies quoted in Table 1 indicate the range of uncertainty in z ₀ measurements; in most cases a mean value is quoted. The most accurate values are likely to be obtained from Antarctic studies, where long homogeneous fetch and a relatively deep SABL are favourable to z0 measurement using sonic anemometers. The vertical profile measurements carry greater uncertainty due to problems of sloping surfaces, atmospheric stability in the surface layer and difficulty in defining the base level for instruments on an uneven glacier surface (Reference MorrisMorris, 1989; Reference Smeets, Duynkerke and VugtsSmeets and others, 1999).

2.3. Measurement of z₀

Direct measurement of z ₀ is possible, using eddy covariance instruments, such as sonic anemometers, which respond to vertical wind-velocity fluctuations on an instantaneous basis. However, these instruments are difficult to deploy on valley glaciers where the surface layer is thin and the instruments are prone to failure and damage (e.g. Reference MunroMunro, 1989; Reference Plüss and MazzoniPlüss and Mazzoni, 1994; Reference Smeets, Duynkerke and VugtsSmeets and others, 1999). Furthermore, the need for careful setting and calibration of the instruments, and long measurement periods mean this method is unsuitable for measuring z0 at a large number of sites across a glacier.

The standard method is to derive z ₀ from the vertical profiles of horizontal wind speed and air temperature, using measurements at two or more heights in the SABL. The logarithmic wind-speed profile can be adjusted for surface layer stability, using Monin–Obukhov similarity theory, enabling z ₀ to be found (e.g. Reference MunroMunro, 1989; Reference Bintanja and van den BroekeBintanja and Van den Broeke, 1994; Reference Plüss and MazzoniPlüss and Mazzoni, 1994; Reference Denby and GreuellDenby and Smeets, 2000; Reference ObleitnerObleitner, 2000; Reference Greuell and SmeetsGreuell and Smeets, 2001). The instrumentation is more robust than for the eddy covariance method, and hence more suited to measurement over glaciers, but the calculation of z ₀ is very sensitive to errors in the instrument heights. A height error of just 0.1 m may change z ₀ by an order of magnitude, but defining the zero-reference plane for the instruments can prove difficult on a rough glacier surface (Reference MunroMunro, 1989; Reference Smeets, Duynkerke and VugtsSmeets and others, 1999). A further problem is the shallow and variable nature of the SABL over mid-latitude glaciers, which may be shallower than the measurement heights (Reference Grainger and ListerGrainger and Lister, 1966; Reference Munro and DaviesMunro and Davies, 1977; Reference MorrisMorris, 1989; Reference Denby and GreuellDenby and Greuell, 2000; Reference Arck and SchererArck and Scherer, 2002). Therefore, long measurement periods are needed to obtain the mean and standard error of the z ₀ value (Reference WieringaWieringa, 1993), and hence the vertical profile method is also unsuitable for recording z ₀ at a large number of sites across a glacier.

Several workers have sought to overcome the z ₀ measurement problem through microtopographic measurements of the glacier surface (e.g. Reference FöhnFöhn, 1973; Reference PricePrice, 1977; Reference MunroMunro, 1989; Reference Arnold and ReesArnold and Rees, 2003). Several relationships between surface roughness element geometry and z0 have been proposed, but it is the empirical relationship of Reference LettauLettau (1969) which has gained widest acceptance in glacial and other studies:

(1)

in which h* is the average vertical extent, or effective obstacle height, of the roughness elements, s is the silhouette area (area of upwind face of an average element) and S is the unit ground area occupied by each element.

The two main challenges for the microtopographic approach are to describe the surface roughness elements and to obtain a representative sample for accurate modelling of the glacier surface element dimensions (Reference MunroMunro, 1989; Reference Smeets, Duynkerke and VugtsSmeets and others, 1999). A sampling method for Equation (1) suitable for glaciers was developed by Reference MunroMunro (1989, 1990). The only measurements required are of the variation of surface elevation made at regular intervals relative to a horizontal reference, in a plane perpendicular to the prevailing wind; h ^* is calculated as twice the standard deviation of the elevations, with the mean elevation set to zero. The number of continuous groups of positive height deviations above the mean elevation defines the frequency, f of roughness elements, and the width of a typical element is defined as the length of the traverse, X, divided by 2f. Equation (1) is solved by substituting s = h*X/2f and S = (X/f)². Despite the simplification of surface element form, the estimation of the silhouette area of the elements using this approach is only ~12% different from its true value (Reference MunroMunro, 1989).

Microtopographic z ₀ measurements have been independently verified using eddy correlation instruments at Peyto Glacier, Canada, by Reference MunroMunro (1989) and with independent profile and eddy covariance measurements at Breidamerkurjökull (Reference Smeets, Duynkerke and VugtsSmeets and others, 1999). However, these comparative measurements were limited to rough ice surfaces. Snow surfaces generally have a fairly isotropic distribution of roughness elements, but the applicability of the microtopographic method to anisotropic ice surfaces has been questioned (Reference FöhnFöhn, 1973; Reference StullStull, 1988; Reference Smeets, Duynkerke and VugtsSmeets and others, 1999), although Reference WieringaWieringa (1993) claims it is applicable in moderately anisotropic situations.

3.1. Field site

Fieldwork was undertaken at Haut Glacier d’Arolla, Valais, Switzerland, a ~6.3 km² valley glacier with an elevation range of ~2550–3500m above sea level (a.s.l.), consisting of an upper basin with northwesterly aspect feeding a glacier tongue flowing to the north (Fig. 1). The main field-data collection periods were May–September 1993 and July and August 1994. Preliminary fieldwork was also conducted in September 1992 and under winter conditions in November 1992 and January and March 1993. The glacier has been the site of several research projects into glacier hydrology, dynamics, meteorology and melt in recent years (Reference RichardsRichards and others, 1996, Reference Brock, Willis, Sharp and ArnoldBrock and others, 2000; Reference Mair, Nienow, Sharp, Wohlleben and WillisMair and others, 2002; Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others, 2004). Below about 3000 m a.s.l. the surface gradient is shallow (generally <10º), but the upper accumulation area contains steep icefalls, particularly on the north face of Mont Brulé. Most of the glacier’s surface can be accessed with relative ease and safety, enabling changing surface conditions to be monitored over large areas.

Fig. 1. Site map of Haut Glacier d’Arolla. The rectangle encloses the area of the glacier displayed in Figure 3.

3.2. Monitoring glacier-wide and seasonal z ₀ variations

To determine glacier-wide variations in z ₀ and related changes in surface conditions, 68 sample points, ranging in elevation from 2572 to 3002 m a.s.l., were established (Fig. 1). Measurements could not be made safely above 3000m a.s.l. due to steep slopes and crevasses. The western margin of the glacier tongue was not sampled as it is completely moraine-covered. Preliminary fieldwork in 1992 revealed that the spatial variability of z ₀ was greatest at low elevations. Accordingly, the spacing of sample points was increased from ~-50m on the snout to ~-200m in the upper basin. The entire network could be sampled in 2–3 days, producing an almost instantaneous picture of spatial z0 patterns. The sample point locations were surveyed onto the Swiss Grid using a Geodimeter 400 total station.

The network of points was sampled at 2–3 week intervals throughout the 1993 ablation season (Table 3). The proportion of the sample points monitored increased during the ablation season in response to the increasing variability in surface conditions. Two glacier-wide surveys were also conducted during the 1994 ablation season, at a smaller number of sample points (Table 3), to enable the broad patterns of z ₀ variation during 1993 to be compared with those during a second ablation season. Measurements were also made at higher spatial resolution over areas between sample points in 1993, prior to glacier surveys 2 and 4 to assess small-scale z ₀ variation. To study the impact of new snowfall and its subsequent melting on z0, additional point measurements were made on the days following summer snowfalls on: 21 May, 3 and 13 June and 28 August 1993, and 3 September 1992.

Table 3. Dates and number of points sampled in 1993 and 1994 glacier surveys

3.3. Microtopographic measurement of z₀

Surface microtopography was measured manually using a 3 m horizontal reference pole and a metal tape measure. The pole was made out of a hollow plastic tube, rectangular in cross-section, and marked at 100 mm intervals along its length. At each sample point, surface microtopography was measured by placing the reference pole horizontally on the glacier surface, perpendicular to the prevailing wind direction, and measuring the distance from the base of the reference pole to the glacier surface to the nearest millimetre at horizontal intervals of 100 mm to generate a 3 m profile. The 30 vertical distance measurements generated were substituted into Equation (1) following Reference MunroMunro (1989, 1990) to calculate z ₀ for the sample point. This method is quite insensitive to measurement errors. An error of ±5mm at one, several or all of the 30 vertical distance measurements made along the reference pole varies z ₀ by at most ±3%, which results in an uncertainty in ln(z₀) of <1%. Nevertheless, great care was taken to stop the pole from sinking into soft snow surfaces when the pole was supported by only a few points of contact with the surface beneath.

At Haut Glacier d’Arolla, as on most valley glaciers, it was assumed that the prevailing wind was constrained by large-scale topography and flowed either straight up- or down-glacier. Analysis of wind direction data recorded at an automatic weather station (AWS) located just in front of the glacier snout supports this assumption, with 90% of recorded hourly wind directions in the ablation season within ±30° of two principal modes (up- and down-glacier).

To test whether a 3 m horizontal profile is long enough to generate a representative z ₀ value, and whether microtopographic z ₀ is independent of the length of the profile used, pairs of z ₀ measurements were made at 20 sites using both 3 and 9m profiles. This comparative sample included fresh snow, 3 day old snow, rough snow and ice and debris surfaces. No significant difference was found between z ₀ values calculated from 3 and 9 m profiles (t test, pH₀ < 0.05). Furthermore, at seven sites 3m horizontal profile measurements of z ₀ were compared with measurements of z ₀ using 5, 6, 12 and 15 m profiles. Although some small differences occurred, there was no systematic variation between z ₀ generated from 3 m and longer profiles. Based on these data, microtopographic measurement of z ₀ is independent of the length of profile, for lengths between 3 and 15 m.

The roughness pole technique was also able to record small-scale microtopography, since vertical height deviations were measured to the nearest millimetre. Thus, z ₀ was also measured over very smooth snow surfaces, during the winter and following fresh snowfall.

3.4. Explaining z₀ variation through surface properties and meteorological variables

An AWS was installed ~200 m in front of the glacier snout at 2547ma.s.l., and operated continuously throughout the fieldwork periods, to enable assessment of the effects of meteorological conditions on z0 and the development of z0 parameterizations (LMS in Fig. 1). The AWS recorded halfhourly averages of 1 s samples of incoming shortwave radiation (Wm^–2), air temperature (°C), relative humidity (%), wind speed (ms^–1) and direction (°) at 2 m height. An identical meteorological station (UMS in Fig. 1) was located on the glacier at 2884 m a.s.l. from 4 July to 25 August 1993 and from 5 July to 23 August 1994. Data from this station were used to determine the local temperature and incoming shortwave radiation lapse rates to extrapolate air temperature and incoming shortwave radiation to all glacier sample points. To determine the relationship between z ₀ and accumulated melt, regular measurements of surface lowering and snow density were made at 16 sample points along the glacier centre line, using ablation stakes. To investigate the relationship of z ₀ to snow depth, snow depth was also recorded using a 3 m avalanche probe (error = ±10 mm).

3.5. Comparison of microtopographic z₀ with wind profile z₀ measurements

To determine whether the microtopographic measurements generated reliable z ₀ values, comparison was made with wind profile z ₀ measurements derived from horizontal wind- speed and temperature profiles recorded over melting snow and slush, between 9 and 29 July 1994, and over melting ice, between 11 and 24 August 1994 (snow and ice profile sites in Fig. 1). The specific objectives of the comparison were to: (i) determine whether microtopographic z0 values agreed with profile z ₀ values recorded over the same surface and (ii) determine whether microtopographic measurements made in a plane perpendicular, or parallel, to the prevailing wind corresponded with the profile z0 over an anisotropic surface. Microtopographic measurements were made over the area upwind from the snow profile site (SPS) on 30 occasions between 9 and 29 July 1994 and on 26 occasions over the area upwind from the ice profile site (IPS) between 11 and 24 August 1994. The upwind area location was determined by the dominant wind direction in the previous 24 hours recorded in the wind profile measurements.

At the SPS and IPS, horizontal wind speed and air temperature were measured at 0.5 and 2.0 m above the surface. Samples at 1 Hz were recorded and averaged at 10 min intervals on a data logger (Campbell Scientific Inc., model CR10, USA). A wind vane (Vector Instruments, model W200P, UK; precision 6°, threshold wind speed 0.6m s^–1) recorded wind direction at 1.0 m height. Wind speed was measured using pulse output type anemometers (Vector Instruments, model A100M, UK; threshold wind speed 0. 15 m s^–1, precision 0.1 m s^–1) and air temperature was measured using resistance temperature-curve matched thermistors (precision 0.4°C) mounted in naturally ventilated radiation shields (Environmental Measurements Ltd, UK). All instruments were mounted on thin aluminium arms (25 mm diameter) supported by a 3 m steel mast (25 mm diameter). A plastic sleeve was drilled into the surface at each site and the base of the mast rested inside the sleeve, supported by a steel screw. Holes were drilled in both the mast and the sleeve at 100 mm intervals, which enabled the mast to be lowered regularly. This, together with adjustments to the heights of the aluminium arms, ensured that the instruments were kept at an approximately constant distance from the glacier surface as it melted.

The 10 min averaged profile data were assembled into half-hour mean datasets, and Monin-Obukhov similarity theory (e.g. Högstrom, 1988; Reference GarrattGarratt, 1992) was used to solve iteratively for friction velocity, u*, and temperature scale, T*, initially setting bulk stability corrections for momentum, α_M, and heat, α_H, to zero:

(2)

(3)

where k is von Kármán’s constant (0.40) and Pr is the Prandtl number (0.95). The subscripts 1 and 2, respectively, refer to the lower and upper level measurements of wind speed, u, and temperature, T, at height, z. Then, the first set of u* and T* values was used to make an initial estimate of the Monin- Obukhov length, L:

(4)

in which is the mean absolute temperature in the surface layer and g is acceleration due to gravity. The sequence of calculations was repeated, with α_M = α_H = 5, using each new value from Equation (4), until there was no further change in u* and T*. A range of stability correction functions exists, but the use of Equation (4) with α_M = α_H = 5 is consistent with the experience that various models give virtually identical results in near-neutral conditions (Reference AndreasAndreas, 2002). Near-neutral conditions are here defined for 0 > z/L < 0.03, taking z to be the height of the upper measurement level. Rearrangement of the log–linear wind profile for use in this range of z/L yields:

(5)

to generate a z ₀ value for each dataset.

In addition to near-neutrality, profile z ₀ measurements were only used in the subsequent analyses if the following criteria were met:

1. 1. Height of maximum wind speed >2 m, assumed when u ₂ > u ₁.

2. 2. Non-obstructed airflow over fetches of at least 500 m.

3. 3. Natural ventilation of radiation shields with u ₁ > 3.5 m s^–1(>4.5ms^–1 for reflected shortwave radiation >50Wm^–2) to minimize radiative heating effects (Reference Georges and KaserGeorges and Kasser, 2002). The manufacturer of the radiation shields specifies an error of 0.4°C for a shortwave radiation flux of 1000 Wm^–2 at a wind speed of 3 ms^–1, hence the temperature measurement error is estimated to be ≪0.4°C.

4. 4. Wind direction range <20° to allow close identification of the upwind surface cover for comparison with microtopographic measurements.

4.1. Surface conditions at the snow and ice profile sites (SPS and IPS)

Initially, at the SPS there was a snowpack of ~1 m depth marked with surface ablation hollows with vertical relief of ~0.1 m and horizontal spacing of ~0.5–1.0 m. The hollows formed a fairly regular pattern with no obvious alignment along or across glacier. After 20 July the hollows collapsed as the remaining snowpack turned rapidly to slush, presenting a much smoother surface with vertical dimensions of roughness elements ~0.01 m. No measurements were recorded between 20 and 24 July due to a power supply failure. Given the obvious difference in microtopography between the rough snow (before 20 July) and smoother slush surfaces (after 24 July), data for these periods are analyzed separately below. Microtopographic measurements were taken perpendicular and parallel to wind directions on snow, but only in the perpendicular-to-wind plane on slush, due to the uniform nature of this surface. The ice surface at the IPS was characterized by foliation bands which formed a series of parallel hummocks and troughs, aligned along glacier, of height ~0.2m and horizontal spacing ~1.0m. Hence, the long axes of the surface roughness elements had strongly preferred orientation aligned with the dominant up- and down-glacier winds. The hummocks continued for over 500m up- and down- glacier from the IPS, but were interrupted every 10–20 m by narrow troughs cutting transverse to the ridges, which probably marked the locations of former crevasses. In contrast to the SPS there was no visible change to the surface microtopography over the measurement period. Apart from a few days of cyclonic weather, conditions were predominantly fine throughout.

4.2. Profile ln(z₀)) values

The ln(z _0p) values generated from the profile measurements are plotted against z/L in Figure 2. On snow and ice the ln(z _0p) values are fairly tightly scattered, between 0.79 and 1.81 mm, and 0.63 and 2.73 mm, respectively, but on slush the ln(z _0p) values display larger scatter between –1.45 and 1.38 mm. None of the ln(z _0p) results show any trend across the stability range.

Fig. 2. Wind-profile derived ln(z ₀) values plotted against z/L for snow, slush and ice surfaces. The ranges of wind speed (u) and temperature (T) corresponding to the ln(z _0p) values are: snow, u = 4.1-8.1 ms^–1 and T = 0.1-5.1 °C; slush, u = 3.5-5.6m s^–1 and T = –0.1 to 1.8°C; ice, u = 5.1-12.1 ms^–1 and T = –2.1 to 3.6°C. The ranges of wind speed and temperature differences between upper and lower measurement levels for each set of ln(z _0p) values are as follows: snow, u ₂–u ₁ = 1.6-2.2ms^–1, T ₂ –T ₁ = 0.3–1.2°C; slush, u ₂–u ₁ = 0.9–1.6ms^–1, T ₂ –T _i = 0.1–0.3°C; ice, u ₂–u ₁ = 1.5–2.9 ms^–1, T ₂–T ₁ =0.1–1.7°C.

4.3. Comparison of microtopographic and wind-profile ln(z₀)) on snow

The mean ln(z _0m) values from perpendicular (ln(z _0m) = 0. 84mm) and parallel (ln(z _0m) = 0.48 mm) microtopographic profiles are slightly lower than the mean ln(z _0p) value of 1.27 mm, but there is a large overlap in the ranges (mean ± 1 standard deviation of the mean) of ln(z _0m) and ln(z _0p) (Table 4). The ln(z _0p) and ln(z _0m) ranges are similar at the upper end, whereas the bottom end of the ln(z _0m) range is much smaller than the lowest ln(z _0p) value. Statistically, there is no significant difference between ln(z _0m) in perpendicular and parallel profiles (t test, pH₀ > 0.05), as expected from the homogeneous nature of the snow surface. Visually, the mean perpendicular ln(z _0m) value corresponds most closely with the mean ln(z _0p) value.

Table 4. Comparison of wind-profile and microtopographic ln(z ₀) over rough snow, slush and ice surfaces at the snow and ice profile sites; σ – standard deviation of the sample

4.4. Comparison of microtopographic and wind profile ln(z₀) on slush

The mean ln(z _0p) value of –0.13 mm on slush is significantly lower than the corresponding value on the rough snow surface (t test, pH₀ < 0.001;Table 4). The mean ln(z _0m) value of –0.42 mm corresponds closely with the mean ln(z _0p) and the ln(z _0m) range is completely within the ln(z _0p) range (Table 4).

4.5. Comparison of microtopographic and aerodynamic ln(z₀) on ice

The mean ln(z _0m) value from perpendicular microtopographic profiles of 1.94mm is almost exactly equal to the mean ln(z _0p) value of 1.93 mm, and the ln(z _0m) range for perpendicular profiles is entirely within the ln(z _0p) range (Table 4). In contrast, the mean ln(z _0m) value of –0.13 mm from parallel microtopographic profiles is significantly lower than the equivalent ln(z _0p) and perpendicular microtopographic profile ln(z _0m) values (t test, pH₀ < 0.0001;Table 4).

4.6. Comparison of microtopographic and aerodynamic ln(z₀): discussion

The results support the application of the microtopographic method to measurement of ln(z ₀) over melting glacier surfaces. While the range of mean ln(z ₀) values in the study was not very large, it spans the 0.1–1.0mm z ₀ scale, and surface types, typical of glacier surfaces during the ablation season (Table 1). The ln(z _0p) values (and perpendicular ln(z _0m) values) are significantly different between the rough snow, slush and ice surfaces (t test, pH₀ > 0.5), indicating that these are distinct surface types with their own characteristic aerodynamic roughness length values. The ln(z _0m) values generated from profiles made perpendicular to the prevailing wind are statistically the same as the ln(z _0p) values recorded over the same surface type (t tests, pH₀ < 0.01;Table 4). However, ln(z _0m) values recorded from profiles parallel to the prevailing wind were significantly lower than ln(z _0p) on the anisotropic ice surface (t test, pH₀ < 0.0001;Table 4), but similar to ln(z _0p) on the more isotropic snow surface. This implies that, where roughness elements have a strong orientation, microtopographical measurements made in a wind-parallel plane do not effectively record the upwind face areas of the surface roughness elements; typically the areas are underestimated. No glacier surfaces were encountered where the long axis of roughness elements was aligned across glacier. It cannot therefore be determined whether parallel or perpendicular microtopographic measurements would correspond to ln(z _0p) where microtopography is rougher in the parallel- to-wind direction than in the perpendicular-to-wind direction. Such a configuration of roughness elements is not likely to be common on mountain glaciers, however, since ice dynamics and the action of meltwater tend to generate ridges and troughs aligned along glacier (Reference Goodsell, Hambrey, Glasser, Nienow and MairGoodsell and others, 2003). Hence, microtopographic ln(z0) measurements should be made using roughness pole profiles aligned perpendicular to the prevailing wind, as defined in section 3, particularly where roughness elements do not form a homogeneous pattern.

The results suggest that accurate (by comparison to vertical wind profile measurements) ln(z ₀) values can be obtained from samples of about six microtopographic measurements. On both slush and ice surfaces, the ln(z _0m) range was smaller than that of ln(z _0p) (Table 4). On both of these surfaces there was little spatial variation in the vertical dimensions of surface roughness elements. In contrast, on rough snow the range of ln(z _0m) was larger than that for ln(z _0p) (Table 4). On this surface there was some spatial variation in the vertical extent of roughness elements, i.e. ablation hollows developed to varying sizes over different areas upwind from the SPS. It appears that the larger element sizes controlled ln(z ₀) on this surface given the close correspondence between the upper range of ln(z _0m) and the mean ln(z _0p) (Table 4).

Incorrect identification of the base level for the instrument heights is a possible error source in the ln(z _0p) measurements. Reference MunroMunro (1989) added 0.17 m (the typical vertical extent of the surface roughness elements) to instrument heights in the calculation of ln(z _0p) at Peyto Glacier, while Reference AndreasAndreas (2002), in a re-analysis of the same dataset, deemed such a height adjustment to be unnecessary. Doubt over the instrument heights adds uncertainty to ln(z _0p) values. An adjustment to instrument heights of ±50 mm leads to a large change in the mean ln(z _0p) value (Table 5). The height uncertainty is not great enough, however, to explain the large decrease in ln(z0p) between rough snow and slush at the SPS. Reduction of the instrument heights by 0.1m (the typical vertical extent of roughness elements) over rough snow reduces the mean ln(z _0p) to 0.41 mm, which is still much larger than the mean ln(z _0p) recorded over slush.

Table 5. Variation in mean ln(z _0p) (mm) and mean z _0p (mm) with adjustment to instrument base height level for snow, slush and ice surface types

6.1. Snow ln(z₀): ln(Z_0S)

The following independent variables were used to explain ln(z _0S) variations: accumulated melt (M _a ; millimetre water equivalent (mm w.e.)), accumulated daily maximum temperature (T_a; °C), accumulated daily mean incoming shortwave radiation (R _a ; W m^–2) and accumulated days (D _a), each of which increases from a value of zero at the time of the most recent snowfall, and snow depth (d; m). These variables may account for the increasing roughness of snow surfaces with time through the formation of ablation hollows associated with local melt rate variations (Reference HuntHunt, 1993). The variables M _a, T _a and d might also explain spatial patterns of ln(z _0S) through their correlations with the up-glacier decrease in the surface melt rate. There was no significant variation in shortwave radiation receipts between the LMS and UMS and it was therefore assumed that incoming shortwave radiation was uniform across the glacier. Based on the mean difference in temperatures between the UMS and LMS, a uniform lapse rate of 0.9°C per 100 m rise in elevation was applied to the temperature data.

All independent variables are correlated significantly with ln(z _0S), and the strongest correlations are those for the four accumulated independent variables (Table 6). However, the relationships between ln(z _0S) and accumulated variables are non-linear, being characterized by three distinct phases in each case (Fig. 6a–d). At both low and high values of the accumulated independent variables, ln(z _0S) varies little (for ln(z _0S) values of about –5 to –2.5 mm and about –0.5 to 2 mm, respectively), but these phases are separated by a period of rapid increase in ln(z _0S) at medium values of each accumulated independent variable. These graphs suggest that the formation of ablation hollows in a melting snow surface begins slowly, but, once hollows have initiated, their growth proceeds rapidly, until some self-limiting condition is reached. This might occur when shading of the bottom of the hollows, or concentration of impurities in the snow acts to halt their growth (Reference HuntHunt, 1993). It can also be seen that ln(z _0S) increases with decreasing snow depth, although there is large scatter in this relationship (Fig. 6e). At snow depths less than ~0.7 m, there are no ln(z _0S) values lower than –2.0 mm (z _0S = 0.14mm), probably due to a combination of well- developed ablation hollows on old snow surfaces and the influence of the underlying ice microtopography on fresh shallow snow covers.

Table 6. Correlations of dependent variables with independent variables. Dependent variables are: snow ln(z ₀) (ln(z _0S)), ice ln(z ₀) (ln(z0ı)), and ln(z ₀) following snowfall on an ice surface (ln(z _0SI)). Independent variables are: accumulated melt (M _a); accumulated daily maximum temperatures (T_a), accumulated daily mean incoming shortwave radiation (R_a), accumulated days (D_a), snow depth (d), and underlying ice ln(z ₀) (ln(z _0])). See text for full definitions. Correlations significant at the 0.05 level are shown in bold. The degrees of freedom for each correlation are given in parentheses. A dash indicates insufficient data to attempt a correlation

Fig. 6. Relationships between snow ln(z ₀s) and (a) accumulated melt, (b) accumulated daily maximum temperature, (c) accumulated daily mean incoming shortwave radiation, (d) accumulated days and (e) snow depth. (f) Relationship between ice ln(z _0I) and accumulated daily maximum temperature.

Two forms of parameterization were applied. First, a linear equation of the form:

(6)

was applied, in which a and b _x are coefficients and V _x are independent variables. Second, a non-linear equation was applied to explain the stepped form of the relationship between ln(z _0S) and the accumulated independent variables:

(7)

in which b _x are coefficients and V is an independent variable. A stepwise regression procedure was used to identify the relationships that explain the largest amount of ln(z _0S) variation using any combination of the independent variables (Table 7).

Table 7. Parameterizations of ln(z ₀): coefficient values and summary statistics; R² is the coefficient of determination. The standard error is given in parentheses after each coefficient value

The stepped form of the relationships of ln(z _0S) to M _a, T _a, R _a and D _a is better represented by the non-linear regressions (Equation (7)) than by the linear regressions (Equation (6)). The non-linear parameterization based on M _a explains the largest amount of ln(z _0S) variation, as indicated by its R ² value but, since the parameterizations based on T _a, R _a and D _a are calibrated with much larger datasets, these relationships are probably better parameterizations of ln(z _0S) variation (Table 7). Furthermore, a ln(z0S) parameterization based on M _a may introduce a circularity problem in a melt model, since the melt rate, i.e. the output from the melt model, must be known a priori. Errors in the initial ln(z _0S) value will generate errors in the melt rate which, in turn, might lead to greater error in ln(z _0S), and thus amplify over time. Overall, the most successful parameterization for ln(z _0S) is the non-linear equation based on T _a (Fig. 7):

Fig. 7. Variation of the non-linear ln(z ₀s) parameterization (Equation (8)) and measured ln(z _0S) values, with accumulated daily maximum temperatures since snowfall.

(8)

The non-linear parameterizations using R _a and D _a offer good alternatives depending on data availability and the modelling approach used. For parameterizations using T _a (Equation (8)) under the condition T _a < 0, the minimum asymptotic value of ln(z _0S) of –3.5 mm (z _0S = 0.03 mm) should be applied, as negative number logarithms cannot be found.

6.2. Ice ln(z₀): ln(z_0|)

The independent variables M _a, T _a, D _a, each of which increases from a value of zero at the time the ice surface is first exposed following melting of the overlying snow cover, and elevation, E, were used to explain ln(z _0l) variations. The accumulated variables are included as surface roughness may increase over time due to small- scale melt rate differentials. All measurements made over ice surfaces were initially included in the analyses. Subsequently, the dataset was divided between: (i) initial ln(z _0l) values recorded at each point immediately following melting of the overlying snow cover, in order to examine ln(z _0l) as a function of E alone and (ii) the change in ln(z _0l) from its initial ln(z _0l) value over time, Δln(z _0l), in order to examine temporal ln(z _0l) trends alone.

The relationships between ln(z _0l) and the independent variables are weak; only the positive correlation with T _a is significant (Table 6; Fig. 6f). The initial ln(z _0l) values are not significantly correlated with elevation and ∆ln(z _0l) is not significantly correlated with any independent variable at the 0.05 significance level. Figure 6f shows little evidence for the pattern of increasing ln(z _0l) between June and August, followed by decreasing ln(z _0l) towards the end of the ablation season, which was suggested by Figure 3c–e. Instead, ln(z ₀|) both increased and decreased over time on different parts of the glacier following melting of the overlying snow cover, with no general trends apparent.

An attempt to parameterize ln(z _0l) as a function of T_a was unsuccessful as this variable explained an insignificant amount of ln(z _0l) variation. One approach to ln(z _0l) parameterization in a numerical model is to use the mean ln(z _0l). The mean ln(z _0l) at Haut Glacier d’Arolla is 0.81 mm (z _0l = 2.24 mm). The standard deviation of ln(z _0l) of 0.89mm gives a range of ln(z0I) of 0.08–1.7mm (z0I range: 0.92–5.47mm). Although errors in calculating spatial and temporal variations in turbulent fluxes will arise from using a constant mean ln(z _0l) value, these errors will tend to cancel when making calculations over the ablation season for an entire glacier. An alternative is to sample ln(z _0l) randomly from a frequency distribution defined by the mean and standard deviation of ln(z _0l), in order to better simulate the likely range of turbulent flux values over ice (Reference Brock, Willis, Sharp and ArnoldBrock and others, 2000).

6.3. The ln(z₀) following fresh snowfall on an ice surface: ln(z_0SI)

When fresh snow falls on a rough ice surface, ln(z _0SI) is strongly influenced by the underlying roughness elements if the fresh snow is too shallow to blanket the underlying ice roughness elements. Following snowfalls on ice surfaces ln(z _0SI) was most strongly correlated with the underlying ice ln(z ₀) (Table 6), demonstrating the important influence of the underlying ice topography on the surface roughness of shallow snow covers. There is a tendency for ln(z _0SI) to increase following snowfall, as the fresh snow melts and more of the underlying roughness elements are exposed, as demonstrated by the negative correlation with d and the positive correlation with T _a (Table 6). The tendency for summer snowfalls to be followed by clear days with subzero temperatures and little or no surface melt could explain why ln(z _0SI) is not correlated with R _a or D _a (Table 6).

It is possible to parameterize ln(z _0SI) through a multiple regression relationship (Equation (6)) on the independent variables underlying ln(z _0l), T _a and d, which explains >45% of the variation in ln(z _0SI) (Table 7). If the underlying ln(z _0SI) is not known, ln(z _0SI) can be parameterized using T _a and d alone (Table 7). Hence, the ln(z _0SI) of fresh snowfalls on rough underlying ice surfaces can be parameterized separately from ‘deep’ snowpacks (section 6.1) in a numerical melt model.

7.1. Validity of the microtopographic z ₀ measurement technique

The close agreement of microtopographic z ₀ measurements with independent wind profile z ₀ measurements over snow, slush and ice surfaces (section 4), provides strong support for the use of microtopographic z ₀ measurements over melting glacier surfaces. Indeed, the microtopographic measurements had lower scatter than the profile measurements over slush and ice, despite the careful selection criteria applied to obtain reliable wind profile z ₀ estimates. The microtopographic technique may well be the better of the two techniques over melting glacier surfaces, if the uncertainty of the base level for instrument heights is considered in wind profile z ₀ calculations. Further validation of the microtopographic technique is needed, however, in particular over surfaces where roughness element size is spatially variable. Our results suggest that z ₀ is controlled by the larger-sized elements over such surfaces.

The wind profile z ₀ measurements relied on only two measurement levels and naturally ventilated temperature shields, which represents the minimum instrumentation necessary for this technique, creating difficulty in determining the instruments’ base height. The strict data selection criteria applied, in particular, the requirement of high wind speeds and near-neutral atmospheric conditions, together with careful monitoring of instrument heights throughout the experiments ensure that the resulting profile z ₀ values are reliable, a conclusion which is corroborated by the close similarity of our profile z ₀ values with z ₀ values reported over similar surface types in other studies using more detailed profiles or eddy covariance measurements. Recent improvements to micrometeorological instrumentation deployable on glaciers should facilitate more reliable comparisons between the wind-profile and microtopographic techniques in the future.

The validation of the microtopographic technique was limited to fairly rough surfaces with roughness elements at the centimetre to decimetre scale (z ₀ at the 0.1–1 mm scale), but we demonstrate the microtopographic technique is also applicable to smoother surfaces in the z ₀ = 0.01–0.1 mm range. In order to obtain reliable microtopographic measurements, a 3 m pole is sufficient for surfaces where z ₀ ≤ 10 mm. If a shorter pole is used it is unlikely that a sufficient sample of surface roughness elements will be recorded, while for surfaces where the vertical dimensions of the elements are >1 m, a longer pole should be used. Microtopographic measurements should be made with the pole aligned perpendicular to the prevailing wind direction, particularly where roughness elements’ long axes have a preferred orientation. If the pole is aligned parallel to the wind, the upwind face area of elements may be underestimated, leading to an underestimate of z ₀. The conversion of height deviations recorded in a microtopographic profile to z ₀, following the method of Reference MunroMunro (1989), is based on a simplification of element forms into regular cube shapes. This is a necessary generalization given that the original formula of Reference LettauLettau (1969) was developed by placing bushel baskets on a frozen lake, and due to the difficulty of measuring and converting irregularly shaped elements into a z ₀ value.

The microtopographic method has a distinct advantage over other z ₀ measurement methods in that it enables repeated measurements to be made at many points across a glacier surface. Given the large spatial and temporal variations in z ₀ recorded across Haut Glacier d’Arolla during two ablation seasons, such a sampling strategy is essential to generate representative z0 values for the modelling of turbulent fluxes and surface melt rate variations.

7.2. Seasonal patterns of z ₀ variation: description and parameterization for numerical melt models

Values of z ₀ in the range 0.01–0.10mm, recorded at Haut Glacier d’Arolla over fresh snow and during the winter, are lower than those previously reported for valley glaciers, but similar to those recorded over Antarctic snow surfaces (Table 2). On older melting snow our z ₀ values of 0.1–5 mm are similar to those reported for other mountain glaciers. For ice, the z ₀ values at Haut Glacier d’Arolla are similar to old melting snow surfaces, but slightly larger (mean = 2.24 mm) with a higher upper limit of ~10mm. The ice z ₀ values at Haut Glacier d’Arolla are similar to those reported for other glaciers, although well below the extremes reported by Reference Duynkerke and Van den BroekeDuynkerke and Van den Broeke (1994), Reference Smeets, Duynkerke and VugtsSmeets and others (1999) and Reference ObleitnerObleitner (2000) (Fig 4; Table 1).

Snow z ₀ exhibited the same clear trend during both the 1993 and 1994 ablation seasons at Haut Glacier d’Arolla. Initially, z0 increased very gradually from sub-millimetre values early in the ablation season or following fresh snowfall, then underwent a period of rapid increase over a few weeks before stabilizing at values of a few millimetres in the mid- to late ablation season. This pattern appears to be controlled by the initiation and growth of ablation hollows in snow, until they reach a self-limiting condition, which in turn may be controlled by solar elevation angle, snow dust content or the magnitude of the turbulent fluxes.

Temporal snow z ₀ variation may be successfully explained by independent variables which accumulate from the time of the last snowfall: melt, daily maximum temperature and incoming shortwave radiation, and days. Parameterizations based on accumulated melt and accumulated daily maximum temperature can also account for along-glacier spatial z0 variations. Non-linear relationships between these two variables, which model the variable rate of snow z ₀ increase over the ablation season, explain >80% of snow z ₀ variation (Fig. 7; Table 7).

Patterns of ice z0 variation at Haut Glacier d’Arolla were less systematic than those for snow. Locally, some marked spatial patterns appeared, related to areas of rough or smooth ice (possibly a function of ice flow and foliation bands) and debris cover, which persisted from one season to the next. Temporally, z ₀ increased over some areas of the glacier during the first half of the ablation season and decreased slightly towards the season’s end, but in other areas different temporal trends occurred. Consequently, it was not possible to explain ice z ₀ variation in terms of any independent variables.

The z ₀ of fresh, shallow snowfalls is strongly controlled by the underlying roughness elements. Variation of snow z ₀ following fresh snowfalls on ice, or rough snow, surfaces can be parameterized separately from ‘deep’ snow in a numerical melt model.

7.3. Implications for numerical melt models and future studies

The lack of a suitable parameterization scheme for glacier z ₀ variations has been widely acknowledged as a problem in the physically based modelling of glacier surface melt rates (Reference BraithwaiteBraithwaite, 1995; Reference Hock and HolmgrenHock and Holmgren, 1996; Reference Hock and NoetzliHock and Noetzli, 1997; Reference Samuelsson, Bringfelt and GrahamSamuelsson and others, 2003). The parameterizations developed here should improve the accuracy of turbulent flux calculations in energy-balance models. The snow z ₀ parameterizations calculate an increase in snow z ₀ of three orders of magnitude between the early and midablation season (or following a mid-summer snowfall), which results in more than a doubling of the turbulent fluxes. This form of snow z ₀ variation (Fig. 7) implies glacier snowmelt models must accommodate a step change in the rate of turbulent heat transfer to melting snow if they are to accurately calculate the time of underlying ice exposure. The temporal pattern of z ₀ variation on melting snow, and its causes, demands further investigation. The errors in turbulent flux calculations resulting from the use of a constant mean ice z ₀ are relatively small, due to the smaller range of variation of ice z ₀ than snow z ₀. Based on a mean z ₀ of 2.24 mm, variation of z ₀ in the one standard deviation range (of 0.92–5.46 mm) alters the turbulent fluxes by at most ±20%.

The transferability of the snow z ₀ parameterization should be tested at other sites. The size to which snow ablation hollows grow is controlled by local environmental factors such as insolation and snow dust content (Reference HuntHunt, 1993), so parameter values in the snow z ₀ parameterization may differ between regions with factors such as latitude and proximity to sources of dust and soot. Further study of ice z ₀ on different glaciers is warranted, particularly if such work aims to identify characteristic ice z ₀ types, which may be usefully incorporated into numerical melt models. Given that areas of debris, and areas of rough ice controlled by dynamics and foliation bands, tend to persist, a single microtopographic survey of a glacier could be used to generate a map of the spatial variation of z ₀ for use in a distributed melt model. There is evidence from this study that areas of relatively rough, or smooth, ice are preserved from one year to the next. Mapping areas of ice z ₀ is time-consuming, even using microtopographic methods, so the application of radar from satellite, airborne or ground-based platforms to glacier z ₀ measurements should be investigated.