1. Introduction

The mass and energy balance of glaciers and ice caps are to varying degrees determined by the turbulent sensible- and latent-heat fluxes. Measurements made on glacier surfaces show that turbulent fluxes contribute 20–40% of the total energy balance (see, e.g., Reference KuhnKuhn, 1979; Reference Ohata, Zhongyuan and LiangfuOhata and others, 1989; Reference Van den BroekeVan den Broeke, 1996; Reference OerlemansOerlemans and others, 1999; Reference Wagnon, Ribstein, Francou and PouyaudWagnon and others, 1999). They are thus an important component especially in relation to the climate sensitivity of the mass balance since it is these energy fluxes, along with longwave radiation, which are directly affected by temperature changes. It is therefore important in energy-balance studies to correctly measure and model these fluxes.

Turbulent fluxes can be determined in the field in several ways. The first is by direct eddy-correlation methods using sonic anemometers. Measurements using these instruments are the most direct, but also the most difficult to carry out, chiefly due to the fragility of the instruments, the continuous maintenance and the problem of interpretation. Attempts to directly measure turbulent fluxes on glacier surfaces (e.g. Reference MunroMunro, 1989; Reference Smeets, Duynkerke and VugtsSmeets and others, 1998; Reference Van der Avoird and DuynkerkeVan der Avoird and Duynkerke, 1999) are generally only for short periods of time and do not always agree with other methods. The second, and most widely used, method is to make use of mean wind, temperature and humidity measurements (made with robust instruments) and convert these mean values to surface fluxes by way of Monin–Obukhov (M-O) similarity theory. The third, and most inaccurate, method is the residual method whereby all the other components of the energy balance are measured, along with the mass balance, and the remaining melt explained by way of the turbulent heat fluxes. This is a useful check, but errors in radiation measurements are generally so large that it cannot be used to determine these fluxes to sufficient accuracy.

The second method described above, using mean values and M-O similarity theory, is the subject of this paper. Recent work with second-order models (Reference DenbyDenby, 1999), which describe in more detail the turbulent structure of the atmospheric boundary layer (ABL) above sloping terrain, has shown that assumptions made in M-O theory are often invalid on sloping glacier surfaces. This is the result of katabatic forcing in the ABL which produces the katabatic, or glacier, wind. Katabatic winds are gravity flows caused by the turbulent cooling of air close to the surface. This cooler, denser air sinks down-slope, producing the well-known glacier wind which is characterized by a low-level wind-speed maximum. The wind-speed maxima found on glaciers with slopes of around 5° can be as low as 2 m above the surface and are thus close to or even below standard measuring heights. As is pointed out in section 2, measurements made on Pasterzenkees show wind maxima below 13 m more than 75% of the time.

The presence of the wind-speed maximum alters the turbulent scaling laws, used in M-O similarity theory for two reasons. The first is the non-negligible turbulent transport term in the turbulent kinetic energy (TKE) budget. Since mechanical production of turbulence by shear is zero at the height of the wind-speed maximum (∂U/∂z = 0), turbulent transport, a second-order term in the TKE budget, will dominate in this region, so similarity arguments which assume this term to be negligible cannot be used. The second reason is that M-O similarity theory assumes a constant turbulent-flux layer, the so-called surface layer where fluxes change in the vertical by <10%. This layer can be very thin when a katabatic wind-speed maximum is present. As mentioned above, wind shear reduces to zero at the height of the wind-speed maximum, as will the vertical flux of horizontal momentum .. This indicates a strong divergence in and thus a very thin surface layer. For a wind-speed maximum of 5 m, M-O theory would then be valid only in the lowest 0.5 m.

In spite of this, M-O theory is still used to calculate turbulent fluxes under these conditions. The question posed here is what sort of influence the presence of a wind-speed maximum has on the bulk and profile methods commonly used to determine these fluxes. By making use of observations and second-order modeling results it will be shown that profile methods will always severely underestimate turbulent fluxes when a wind-speed maximum is present but that the bulk method still gives quite good estimates of these fluxes up to at least the height of the wind-speed maximum.

2. Observations

In this paper, use is made of experimental data from the PASTEX meteo-glaciological experiment carried out on Pasterzenkees, Austria, during summer 1994 (Reference Greuell, Knap and SmeetsGreuell and others, 1997). Six weather stations for determining the energy balance of the glacier were stationed along the central flowline. Most of these had just two measurement heights, 0.5 and 2 m, but at one site, known as A1, a 13 m profile mast with eight measurement heights (0.4, 0.7, 1, 2, 4, 6, 8 and 13 m) and a balloon sounder were placed at a point roughly 1 km from the end of the glacier tongue at 2200 m a.s.l. Reference Smeets, Duynkerke and VugtsSmeets and others (1998) provide more information concerning the measurements made using this mast. The local slope at this site was approximately 3.5°.

During the observational period 16 June–11 August, >2600 half-hourly average observations were made. Of these, 77% showed the existence of a wind-speed maximum below the highest profile level when wind direction was downslope, indicative of katabatic flows. Air temperatures were well above 0°C during the observational period, and the surface was almost always melting. Figure 1b shows the average downslope wind-speed and temperature profiles during a 2 day fair weather period at this site. During this period a wind-speed maximum was present >90% of the time.

Fig. 1. (a) Schematic diagram of the glacier wind showing terms described in the text, (b) Average observed downslope wind speed (squares) and temperature (circles) at site A1 on Pasterzenkees during a 2 day fair weather period. Mast profile measurements are indicated by large shapes, and balloon soundings, made approximately every 3 hours, by smaller shapes.

3. Mean and Turbulent Structure of Katabatic Flows

Katabatic flows on glaciers are driven from beneath by buoyancy, or katabatic, forcing. Warm air overlaying the glacier is cooled by turbulent exchange with the ice or snow surface, becoming denser and resulting in the downslope flow of cool air. This katabatic flow, often referred to as glacier wind, is a dominant feature of most temperate glaciers. The low-level forcing combined with the turbulent exchange of momentum will intrinsically lead to the development of a wind-speed maximum at some level. It is the nature of the turbulence that determines which form the mean wind and temperature profiles will adopt.

A one-dimensional second-order turbulence closure model which simulates the vertical profile of the ABL is used here to help describe the turbulent structure of katabatic flows. The one-dimensional model is discussed in detail in Reference DenbyDenby (1999) where it has been shown to describe quite well the mean and turbulent structure of katabatic flows. It makes use of second-order closures, including closures for the turbulent transport of TKE, so M-O similarity is not assumed in order to calculate the katabatic profiles. Only at the lowest level of the model, z = 0.2 m, is the assumption of similarity used to calculate lowest-level fluxes. The roughness length for momentum in the model is predefined and the roughness length for temperature is based on the surface renewal model from Reference AndreasAndreas (1987). The flux-profile relationships, as defined in section 4, are used for the stability correction at this lowest level.

Equations (1) and (2) show the simplified one-dimensional equation for downslope momentum U and temperature perturbation Θ.

Θ is defined as the potential temperature difference from a predefined ambient temperature profile given as Θ_o(x3) = Θ_o +γ_θx ₃, where x ₃ is the true vertical component, i.e. parallel with the gravitational vector, and where the coordinate system (x, z) is defined as a Cartesian coordinate system orthogonal to the sloping surface with slope angle α. V _g is the geostrophic wind component, and the overbar indicates the turbulent-flux terms, as defined by Reynolds decomposition. These are the vertical flux of horizontal momentum and the vertical flux of potential temperature which is directly converted to the turbulent flux of sensible heat where ρ is the density and C _p is the specific heat of air. A schematic representation of the glacier wind is given in Figure 1a indicating the parameters mentioned above.

In Figure 2a the simulated downslope wind and temperature profiles are shown for a typical case of katabatic forcing on a sloping glacier with surface temperature perturbation Θ_s = −10 K, surface roughness length z _o = 2 mm, temperature lapse rate γ_θ = 3 K km⁻¹ and slope sin(α) = −0.1 (6°). The major features of the mean flow are of course the wind maximum, in this case at a height of around 5 m, and the strong temperature inversion beneath the wind-speed maximum, of approximately 0.5 Km⁻¹.

Fig. 2. (a) Simulated downslope wind speed (U) and temperature-deficit (Θ) profiles for the one-dimensional katabatic case described in the text. The wind-speed maximum is at 5 m height. (b) Simulated and profiles, (c) Simulated TKE budget profile.

On typical temperate glaciers katabatic forcing is the dominant driving force in the momentum budget (Reference Denby and SmeetsDenby and Smeets, in press) close to the surface. This is balanced by the flux divergence of horizontal momentum in Equation (1) resulting in the strong gradient in seen in Figure 2b. will pass through zero at, or at least close to, the height of the wind maximum. Though there is a steep gradient in , the vertical profile for sensible-heat flux, , is almost constant in the region below the wind-speed maximum, varying in this case by just 15% (Fig. 2b). Though this is a one-dimensional simulation and other terms in the temperature budget are also important in a two- or three-dimensional case, such as horizontal advection of temperature perturbation, the basic turbulence structure remains the same.

In order to indicate the role played by turbulent transport in katabatic flow, the TKE budget is also shown in Figure 2c. As the height of the wind-speed maximum is approached, the mechanical production through shear approaches zero, and the turbulent transport term transports TKE into this region. This is where M-O theory, as well as local scaling arguments, break down since they do not take into account the transport terms. Above the wind-speed maximum, where shear dominates the TKE budget, local scaling arguments such as those described by Reference NieuwstadtNieuwstadt (1984) become applicable again.

4. Suitability of Bulk and Profile Methods

Before comparing the bulk and profile methods with observations and model results, it is useful to give a brief description of these two methods. Both the bulk and profile methods use M-O similarity theory which is based on scaling arguments using the turbulent velocity and temperature scales, u _* and θ _*, and the length scales z (height), L (Obukhov length) and z _o (surface roughness length). These turbulent scales define completely, under idealized horizontally homogeneous and quasi-steady-state conditions, the turbulent structure of the surface layer, the lowest region of the atmospheric boundary layer (ABL).

The flux-profile relationships, Equations (3) and (4), are definitions relating the vertical gradients of mean wind and temperature to these turbulent scales.

where and .

Under neutral conditions ϕ _m is defined as being equal to unity, and the von Kármán constant (κ = 0.4) is determined from observations to fulfill this definition. Integration of these relationships, under stable conditions and the assumption of a constant flux layer, leads to the well-known log-linear wind and temperature profiles if a linear relationship is assumed for ϕ _m,h of the following form (Reference GarrattGarratt, 1992):

where .

The constants in Equations (5) and (6) are determined solely from observations and vary to some extent in the literature. We use here a value of 5 for both α _m and α _h and a near-neutral turbulent Prandtl number Pr = 1. Integrating Equations (3) and (4) from level z ₁ to z ₂ leads to the following log-linear profile equations:

and integrating these equations from the roughness-length heights of z _o and z _h to z gives the log-linear bulk equations:

It should be noted that the roughness lengths are the height at which the extrapolated wind and temperature profiles would reach their surface values.

The above equations require a minimum of two profile heights for wind and temperature to derive the four unknowns u _*, θ _*, z _o and z _h assuming the surface temperature Θ_s to be known. The accuracy of the profile method will increase when more profile levels are used and these are fitted using least-squares techniques. For simplicity, in this paper we will limit ourselves to the case where only two profile heights are available; this defines the profile method for determining surface fluxes. The bulk method, on the other hand, makes use of only one height to derive u _* and θ _* assuming the roughness lengths (z _o, z _h) and the surface temperature Θ_s to be known. The application of these methods has been quite successful in the stable ABL, which is the reason for its popularity under most conditions.

We will now take bulk- and profile-derived fluxes calculated from the profile-mast A1 data and compare these to model simulations. We present the results in Figures 3 and 4 for both and where we have used the non-dimensional height z/H (H is the height of the wind-speed maximum) as the scaled vertical axis. In these figures all mast levels have been used to calculate the bulk- and profile-derived fluxes. For the profile method each point represents the profile-derived flux taken from two adjoining mast levels at a height determined by their average. For the bulk method each point indicates the bulk-derived fluxes for a particular mast level. Both datasets have been normalized by the bulk-derived value taken from the 0.7 m mast level using previously derived roughness lengths (Reference Denby and SmeetsDenby and Smeets, in press). This is considered to give the best estimate of the surface flux since we expect that when H ≫ z ≫ z_o, M-O theory will still be valid and so the normalized fluxes should approach unity in this range.

Fig. 3. Profile-derived normalized vertical momentum flux (a) and sensible-heat flux (b) as function of the non-dimensional height z/H (H is height of the wind-speed maximum) derived from observations (dots) and one-dimensional model simulations (continuous line). Both the heat and momentum fluxes have been normalized by the bulk-derived surface flux determined at 0.7 m height. See text for details.

Fig. 4. Bulk-derived normalized vertical momentum flux (a) and sensible-heat flux (b) as function of the non-dimensional height z/H (H is height of the wind-speed maximum) derived from observations (dots) and one-dimensional model simulations (continuous line). Both the heat and momentum fluxes have been normalized by the bulk-derived surface flux determined at 0.7 m height. See text for details.

The model results are calculated by a continuous run of the one-dimensional model where the surface temperature perturbation was allowed to vary slowly by 0.1 K h⁻¹ from 0 to −20 K. In this way a continuous range of simulated katabatic profiles was generated by varying the katabatic forcing. Fixed heights of 2 and 0.5 m within the model were then used to calculate the profile-derived fluxes, and a height of 2 m was used to calculate the bulk-derived fluxes. Both the bulk- and profile-derived model fluxes have been normalized by the model-generated surface flux value. The local slope, α = 3.5°, was used as well as the average lapse rate of γ_θ = 3.5 K km⁻¹.

Even though the simulations are one-dimensional, the basic turbulent structure of the katabatic flow remains the same even with the introduction of advection terms in Equations (1) and (2). As mentioned in section 3, the momentum budget is dominated by the katabatic-forcing and flux-divergence terms even when horizontal and vertical advection are present. The temperature budget, on the other hand, is more strongly affected by the horizontal and vertical advection of temperature perturbation. To test the sensitivity of the results to temperature advection, an entrainment velocity was introduced into the temperature budget, as was carried out in Reference DenbyDenby (1999), which allows the representation of vertical advection within the one-dimensional model. Even under conditions of strong entrainment the normalized profiles shown in Figures 3 and 4 remain essentially unaltered by the introduction of advection. In this regard the one-dimensional results obtained here are quite robust.

The profile-derived fluxes show a large amount of scatter in the observations, which results from the sensitivity of this method to measurement errors and variability. What is clear, though, is that turbulent fluxes are severely underestimated throughout the region below the wind-speed maximum. This is the result of reduced shear and increasing temperature gradients as the wind-speed maximum is approached. This ensures that the local-gradient Richardson number (Ri = g/Θ_o ∂/Θ/∂z/(∂U/∂z)² + (∂V/∂z)²) approaches its critical value (Ri_c = 0.23 in this model), above which all turbulence is suppressed according to M-O theory. This occurs at around half the height of the wind-speed maximum (Fig. 5a).

Fig. 5. Gradient Richardson number (a) and bulk Richardson number (b) as function of the non-dimensional height z/H (H is height of the wind-speed maximum) derived from observations (dots) and one-dimensional model simulations (continuous line). The dotted line in (a) indicates the model’s critical gradient Richardson number of 0.23. The bulk Richardson number is defined as (Ri_B = g/Θ_o(Θ(z) − Θ_s)z/U(z)² + V(z)²

The bulk method, on the other hand, shows far less scatter since gradients need not be determined. Though is underestimated when measurements are made above the height of the wind maximum, this is not as severe as in the profile case. , on the other hand, is slightly overestimated when measurements are made below the wind-speed maximum but reduces quickly as the measurement height increases above the wind-speed maximum. Generally, wind maxima are above a normal measuring height of 2 m, so we are not usually concerned with this region.

How can we explain the ability of the the bulk method to estimate the surface fluxes even in the region of the wind-speed maximum? The bulk method is essentially an integrated form of the profile method and so is less sensitive to variations in gradients brought about by the presence of a wind-speed maximum. As turbulence decreases with decreasing shear near the wind maximum, wind speeds reduce relative to their log-linear form (Equation (9)). This effect is shown in Figure 6 where we have plotted the katabatic-wind and temperature profiles from Figure 2 against logarithmic and log-linear profiles with equivalent surface flux values. This then leads to a reduction in the bulk-derived values in Figure 4a. However, this same reduction in turbulence will increase the temperature gradient, and therefore the temperature, as we approach the wind-speed maximum relative to its log-linear form (Fig. 6b). These two effects tend to cancel each other out, reducing u _* and increasing θ _*, and the resulting bulk-determined heat flux remains fairly constant until the diminishing wind speed above the wind maximum leads to stability corrections severely reducing the apparent turbulence.

Fig. 6. Comparison of model-simulated katabatic-wind (a) and temperature (b) profiles with the logarithmic and log-linear profile forms. The logarithmic and log-linear profiles are calculated assuming surface fluxes equivalent to the katabatic surface value.

For the one-dimensional katabatic case the normalized model profiles generated for Figure 6 are not universal. The height of the katabatic wind-speed maximum (H) is not the only scaling length which determines the form of the normalized curves since the turbulent M-O scaling length L is of a similar order. The dependence on slope for the bulk-derived heat flux is illustrated in Figure 7a where three model runs with differing slopes are plotted. From a pragmatic point of view the slope dependence is of little consequence since the largest deviations from unity occur with weaker slopes where wind maxima are higher, so measurements will always be made well below the wind-speed maximum. It is therefore relieving to see that even for steep slopes, where the wind-speed maximum can be quite low, the bulk-derived heat flux appears to be quite representative of the surface value.

Fig. 7. (a) Variation in the normalized bulk-derived heat flux for three different surface slopes taken from one-dimensional simulations. (b) Maximum wind-speed height as function of maximum wind speed for the same simulations. Points indicate observations from the profile mast A1 on Pasterzenkees. Local slope is approximately 3.5°.

To illustrate this point further, the height of the wind-speed maximum has been plotted as a function of maximum wind speed for the three differing slope angles in Figure 7b. As in all the simulations, the surface temperature deficit has been allowed to vary slowly from 0 to −20 K, and the closed circles in Figure 7b indicate 2 K steps. For wind-speed maxima of <2 m, the temperature perturbation and wind speed are quite low, corresponding to small surface turbulent heat fluxes, especially for lesser slopes, so any errors introduced due to the proximity of the wind-speed maximum will be, on an absolute scale, negligible.

5. Conclusion and Discussion

The suitability of M-O similarity theory under conditions of katabatic flow where a low-level wind-speed maximum is present has been discussed. Though assumptions made in M-O theory are not valid in the presence of a wind-speed maximum, it has been shown that bulk estimates of turbulent heat fluxes give quite reasonable results in the entire region below the wind maximum, whereas profile-derived fluxes severely underestimate the surface fluxes. Only in the lowest region, z/H < 0.5, does the profile method give any value for the surface flux at all, due to stability corrections, and only in a narrower region, z/H < 0.3, do we find values comparable to low-level bulk estimates. For a wind-speed maximum at 5 m this would imply that profile measurements would need to be made below a height of approximately 1.5 m. However, profile fits need to be made well above the average surface-roughness element size which can easily be of the same order on some melting ice surfaces.

That bulk methods are effective in determining turbulent fluxes, even under conditions of katabatic flow, is good news for energy-balance studies since measurements of mean wind and temperature at one level are far easier to carry out than profile or eddy-correlation measurements. There remain, however, two problems in using the bulk method to determine turbulent heat fluxes. The first is the determination of the surface temperature. Though this can be an important point for non-melting surfaces, it is inconsequential when studying melting ice surfaces, as we have done here, since the surface is permanently at the melting point.

The second problem is the determination of surface roughness lengths for momentum, temperature and water vapor which are essential parameters when using the bulk method. Normally roughness lengths are determined using profile methods since this is how they are defined, being the height at which the extrapolated profiles in the surface layer would reach their surface value. In essence, roughness lengths are a parameterization for the inertial sub-layer, a region affected by both viscous and inertial forces, close to the surface. The depth of this layer, and also z _o, is correlated with the surface-roughness element form, height and spacing (Reference GarrattGarratt, 1992).

In practice, the determination of roughness lengths for momentum requires near-neutral conditions of long fetch and a well-developed deep quasi-steady-state boundary layer. These conditions are seldom met on a melting glacier where the strong temperature inversions and katabatic winds prohibit the use of profiles to determine z _o due to stability considerations and the presence of the low-level wind-speed maximum, as shown in section 4. For instance, Reference Smeets, Duynkerke and VugtsSmeets and others (1998) could use just 38 of the 2600 half-hourly mean measured profiles from the same PASTEX dataset in order to estimate surface-roughness lengths. An alternative possibility for determining z _o, proposed by Reference Denby and SmeetsDenby and Smeets (in press), makes use of the dynamics of katabatic forcing to determine the vertical profile of momentum flux from wind and temperature profile data. This method, though shown to give quite good results, requires high-resolution profile masts in order to determine the height of the wind-speed maximum.

This is why various attempts have been made, both in the field and in laboratories, to relate surface-roughness element height and form with the aerodynamic roughness length, at least for momentum (Reference LettauLettau, 1969). The use of microtopographical surveys (e.g. Reference MunroMunro, 1989) to determine roughness lengths for momentum, and of surface renewal theories (e.g. Reference BrutsaertBrutsaert, 1975; Reference AndreasAndreas, 1987) to determine roughness lengths for temperature and water vapor, are much easier and appealing methods. Though these methods are promising, since measuring microtopography is much simpler than directly measuring aerodynamic roughness lengths, estimates based on roughness element data must be considered with care, as pointed out by Reference Smeets, Duynkerke and VugtsSmeets and others (1999), and should be treated as reasonable approximations when aerodynamical data are not available. Once more, confirmation of these methods on glacier surfaces is difficult due to the presence of katabatic flows.

This remains a problem, so best estimates will still have to be used with the available data. Luckily an order-of-magnitude error in z _o will only lead, when the bulk method is used, to roughly a 25% error in surface heat fluxes. This is still a far better alternative than the use of profile methods to determine the same fluxes under conditions of katabatic flow.