Introduction
Recent work on New Zealand glaciers has attempted to link fluctuations in their terminal positions with variations in precipitation and temperature. The findings, however, are sometimes contradictory. For example, Reference Salinger, Heine and BurrowsSalinger and others (1983), who have documented the retreat of the small Stocking Glacier to the east of the Main Divide of the Southern Alps in Mount Cook National Park, found a significant positive correlation with monthly temperature with a lag time of 2 years. There was no clear relationship between glacier behaviour and precipitation. They considered that the large retreat was more a consequence of higher temperatures than lower precipitation rates. By contrast,Reference Hessell Hessell (1983) found that fluctuations of Franz Josef Glacier to the west of the Main Divide were significantly correlated with precipitation but there was no relationship with temperature.
The rational way to resolve this divergence of opinion is to consider glacier mass balance with emphasis on long-term studies (Reference Gellatly and NortonGellatley and Norton, 1984). Estimation of glacier ablation then becomes a major problem and generally requires the use of long-term climatological data.Reference Marcus, Moore and Owens Marcus and others (1985) also stressed the need for further studies of ablation and accumulation processes and their relationships with synoptic climatology.
Within New Zealand there have been several studies that describe mass balance of glaciers (e.g. Reference Goldthwait, McKellar, Wexler, Rubin and CaskeyGoldthwaite and McKellar, 1962;Reference Thompson and Keils Thompson and Kells, 1973; Reference Anderton and Chinn Anderton and Chinn, 1978; Dickson, unpublished; Harding, unpublished). Recent work has made short-term estimates of the energy fluxes that drive the melting process over Franz Joseph Glacier Reference Owens, Marcus and Moore(Owens and others, 1984; Reference Marcus, Moore and OwensMarcus and others, 1985).Energy balances for melting snow have also been inferred or measured by Reference Fitzharris, Stewart and HarrisonFitzharris and others (1980), Reference Prowse and OwensProwse and Owens (1982),Reference Moore Moore (1983[a]), andReference Moore and Owens Moore and Owens (1984).
This paper presents and evaluates energy-balance and bulk-aerodynamic approaches for estimating the ablation using measured, or inferred, meteorological data. The models are tested with data obtained from Ivory Glacier(Fig. 1) during the International Hydrological Decade. A further paper (Hay and Fitzharris, in press) uses the models in an investigation of the relationships between glacier melt and synoptic climatology.
Fig. 1. Ivory glacier, new zealand, showing location and topography.
The Models
The surface-energy balance of a melting glacier can be written
where QM is energy used in ice melt, Q* is net radiation, QR is energy from sensible heat of rain, QH is sensible-heat transfer, and QE is latent heat of vaporization transfer. It is usual to express these terms as flux densities (W m−2 or MJ m−2 d−1). Energy directed toward the surface is considered positive and that away is negative.
QM can be related to the depth of ice melted by
or in water units by
where M i is ice melt (e.g. md−1, depth), Mw is ice melt (e.g. m d−1, water equivalent), L f is latent heat of fusion (0.333 MJ kg−1), ρi is density of glacier ice (905 kg m−3), and ρw is density of water (1000 kg m−3).
Once all the terms on the right-hand side of Equation (1) are measured or well estimated, glacier ablation can be determined. Q* can be measured directly. QR is calculated, following Reference PatersonPaterson (1981), as
where cw is specific heat of water (J kg−1°K−1), R is rainfall rate (kg m−2s−1 or mm s−1), TR is temperature of the rain (°K), and TS is temperature of the ice surface (°K).
QH and QE are obtained using the following equations (see Оkе, 1970;Reference Granger and Male Granger and Male, 1978;Reference Paterson Paterson, 1981, for details). In the turbulent boundary layer above a glacier surface:
where z is measurement height above the ice surface (m), ρa is density of air (kg m−3), c P is specific heat at constant pressure (J kg−1 °K−1), K H is coefficient of eddy conductivity (m2 s−1), T is air temperature above the ice surface (°K), Lv is latent heat of vaporization (J kg−1), K W is coefficient of eddy diffusivity for water vapour (m2 s−1), q is specific humidity above the ice surface (g kg−1), Τ is vertical flux of horizontal momentum (kg m−1 s−2), K M is coefficient of eddy viscosity (m2 s−1), and u is wind speed above the ice surface (m s−1).
K M, K H, and K W vary with height but it can be shown (e.g.Reference Deardorff Deardorff, 1968; Reference Price and DunnePrice and Dunne, 1976;Reference Moore Moore, 1983[b]) that
where K̄ H and K̄ W are bulk-exchange coefficients (dimensionless) for the layer (z 2 - z 1). The principle of similarity is usually invoked in order to evaluate these exchange coefficients. That is
Two approaches will be used to estimate K, with subscripts B and E distinguishing the resulting values.
Bulk-aerodynamic approach
Following the derivation presented by Reference MooreMoore (1983[b])
where k = 0.4 = von Karman’s constant, z0 = is the roughness length (m),
and Rb is the bulk form of the Richardson number given by
where g is acceleration due to gravity (9.8 m s−2).
A common assumption is that T0 equals Ts, which for a melting glacier equals 273.2 °K, thereby enabling K B to be calculated from measurements of T and u at a single height above the glacier surface.
Energy-balance approach
If independent measurements of ice ablation are available, QM may be obtained from Equation (2) and K E can then be calculated as follows. From Equation (1), for the convention of positive fluxes toward the surface
but from Equations (7), (8), and (9)
Therefore, combining Equations (13) and (14),
and
For a melting glacier surface, it is asumed that q0 = qS (= 4.5 g kg−1 at a pressure of 850 mbar), thereby enabling K E to be calculated using QM, Q*, QR, and measurements of T, q, and u at a single height above the glacier surface.
For both approaches, initial estimates can be made for periods of about 1 d, based on either hourly or daily averaged data. The estimates may then be averaged for intervals of a week or more to obtain K for time periods commensurate with the typical needs for estimates of glacier ablation.
Assumptions
The bulk-aerodynamic approach assumes the following:
(a) Log-linear profiles for temperature, humidity, and wind speed, and hence the turbulent exchanges QH and QΕ) are constant with height. During light wind conditions, the profiles can deviate markedly from log-linear forms due to radiative-flux divergence and a lack of vertical mixing Reference Halberstam and Schieldge(Halberstam and Schieldge, 1981), but in these circumstances QH and QE tend to be small. Turbulent fluxes are constant with height over extensive, homogeneous surfaces. However, on small glaciers, such as Ivory Glacier where the fetch for all directions except the north is less then 250 m, horizontal advection may occur.Reference Marcus, Moore and Owens Male and Gray (1981) suggested that the height of the constant-flux layer is one-hundredth of the fetch. This indicates that advection will not be a problem for Ivory Glacier measurements provided the sensors are less than 2 m above the ice surface.
(b) Equality of the scale heights for temperature, humidity, and wind speed. In a detailed analysis, Reference MooreMoore (1983[b]) showed that this assumption is often unjustified, but noted that any resulting bias might “average on the order of 10%, and could be masked by instrument error or uncertainty in the value of z 0”.
(c) Similarity of exchange coefficients (Equation (9)).Reference Grainger and Lister Grainger and Lister (1966) believed this is “applicable to observations in the lowest 2 m above a melting glacier surface”, but this has been questioned by some workers (e.g. Reference Marcus, Moore and OwensMale and Gray, 1981).Reference Moore Moore (1983[b]) reviewed the issue at some length and concluded that “the evidence against equality of the eddy diffusivities in the stable range is inconclusive, so the assumption that KH = KW = KM is as reasonable as any other, but carries some degree of uncertainty”.
(d) Steady-state conditions apply. To meet this assumption, measurements are made frequently during the day but are time-averaged to determine K . The effects of determining K over increasing time periods is discussed further in the results.
(e) No heat is exchanged with the underlying ice. This term is excluded from Equation (1) and the assumption is valid in summer when the glacier ice in the ablation zone is isothermal. For example, in their ablation measurements on Franz Joseph Glacier, Reference Marcus, Moore and OwensMarcus and others (1985) found that over periods of 1–3 d heat conduction into the ice was “identically zero or negligible”.
(f) The glacier surface is melting. During the summer this is valid for most of the time on Ivory Glacier, except for occasional nocturnal freezing.
(g) There is no mass exchange due to wind action (e.g. blowing snow) or horizontal movement of melt water.
Reference MooreMoore (1983[b]) noted that the few studies which have examined the validity of these assumptions for a melting snow surface in non-ideal field conditions confirm the appropriateness of the bulk-aerodynamic approach, but he also stressed the need for further investigations.
The energy-balance approach is less restrictive in its assumptions in that (a) and (b) need not apply. It is also more robust because it implicitly accounts for the “average” influence of such features as nocturnal freezing, advection, stability, and surface roughness.
Methods
Measurements to test Equation (1) are available for Ivory Glacier, on the west of the South Island (Fig. 1). It is a high input-output cirque glacier of area 0.8 km2 and annual precipitation exceeding 9000 mm. It was chosen for the IHD programme for the reasons discussed byReference Chinn and Bellamy Chinn and Bellamy (1970). The glacier faces south and its catchment ranges in elevation from 1400 to 2080 m. A detailed description of the glacier, including photographs, and a summary of the IHD measurements have been given byReference Anderton and Chinn Anderton and Chinn (1978).
First, estimates of K are obtained using the bulk-aerodynamic approach (Equation (10)) and the energy-balance approach (Equation (16)) together with data from the “calibration period” (CP) which ran from 29 January to 14 February 1973. Secondly, the models and the estimated values of K are evaluated using independent ablation data from the “validation period” (VP) which extended from 5 January to 14 February 1972.
Ice ablation (Mi) was measured by observing the change in ice level at several wooden poles placed in an approximately 100m2 regular grid at elevation 1500 m in the ablation zone (Fig. 1). During the CP, nine poles were used while during the VP there was one less. Q* was obtained with a net radiometer (CSIRO-Middleton) located 1 m above the ice surface. To evaluate QR, it is necessary to know R, TS, and TR (Equation (3)). Values of R were provided by a recording rain gauge. T S was assumed to be at 0°C during summer rainfall. A thermistor inserted into the stem of the rain gauge’s funnel was used to measure TRduring parts of two storms which collectively resulted in 48.4 mm of rain. Comparisons with the air temperature (TA) yielded the relationship: TR = 2.42 + 0.78TA. Since the air was saturated during the precipitation events, the air and wet-bulb temperatures were equal. Hence, this approach presumably has an advantage over equating TR and the wet-bulb temperature, as suggested by Reference AndersonAnderson (1976) and others.
To implement Equations (7), (8), and (10) or (16), measurements of u 1, T 1 and q1 are required. T 1 and q 1 were measured using screened thermistor and dew-point temperature probes (Atkins), respectively. These were located 1.0 m above the ice surface during the CP and at 1.2 m during the VP. Wind speed was obtained from a cup anemometer (Casella) at 1.0 m during the CP but at 2.8 m during the VP. Data from the VP were adjusted down to 1.2 m using the power law
For the weak or moderately stable conditions typical of the atmosphere above a melting glacier surface,Reference Marcus, Moore and Owens Male and Gray (1981) assumed a value of 0.167 for a and this value was adopted in Equation (17). Since wind speeds were also measured at 4.7 and 6.9 m during the CP, it is possible to study the relationship between K and the wind speed and the wind gradient although, as discussed earlier, measurements at these heights can be influenced by advection.
Measurements were made hourly from 06.00 to 18.00 NZST and every 3 h at night. All fluxes, and the values of MI and estimates of K were usually referenced to 24 h periods but occasional inclement weather resulted in measurement periods as long as 38 h and as short as 11 h.
Errors
Probable measurement errors for daily values are listed in Table I. Although ice ablation (Mi) can be measured to ±0.005 m at each pole, over short (hourly) time intervals this error can be a large proportion of the total error in estimating K However, at typical summer-melt rates, it becomes small over periods of a day or more. Over time, the errors in short-term ablation measurements are largely compensatory (Reference Müller and KeelerMüller and Keeler, 1969) and can be further reduced by pole replication. For example. Figure 2 shows melt for the CP at six of the nine poles. While systematic differences between the poles do exist, averaging the daily ice-ablation readings over the pole grid reduces the error associated with such measurements.
Table I. Estimated measurement errors over a 24 h period
If the errors of Table I are maximized and minimized in the same direction (i.e. no compensation of the individual error terms), the errors in estimating K B and K E from Equations (10) and (16), and thence Mi or Mw from Equations (1) and (2), cover the range ±18%. In practice, measurement errors will typically contribute less than this since the signs of the individual measurement errors are unlikely to be always the same; hence there will be some compensation.Reference Moore Moore (1983[b]) showed that values of the turbulent fluxes calculated using the bulk-aerodynamic method will be within about 20% only if the roughness length (z0) is estimated to within 1 or 2 mm. This is a considerable challenge given the errors associated with estimating z0 from either published values or wind-profile data. In the latter case, the calculated values are highly sensitive to errors in the wind speed.
Fig. 2. Accumulated ablation as measured at six of the nine poles during the calibration period (CP).
In both approaches, other errors will arise because some or all of the assumptions are not met, particularly those related to advection. The size of these errors is difficult to quantify and even order-of-magnitude estimates of the errors resulting from departures from the assumptions of the bulk-aerodynamic method are unavailable (Reference MooreMoore, 1983[b]).
Results
Energy Sources and Measured Ablation Rates
Attempts to model turbulent-energy exchanges and subsequently the glacier-ablation rates are best assessed in the context of the magnitude of these processes. Measured glacier melt (Mw) averaged 38 mm d−1 over the two summer observation periods (Tables II and V). There was considerable daily variation from less than 10mmd−1 to over 70 mm d−1 Energy for melt was almost equally shared between Q* (52%) and the turbulent fluxes QH and QE (46%). Energy from rain made a small contribution (2%), but again there was substantial day-to-day variability (Tables II and V). Our later paper will contain a more detailed analysis of the energy sources for ablation and the associated synoptic climatology. However, by way of a preliminary comparison, the short-term measurements ofReference Marcus, Moore and Owens Marcus and others (1985) at 600 m elevation on Franz Josef Glacier gave Mw as 65 mm d−1 in December 1981. Q* contributed about one-third and QH + QE) about two-thirds of the energy required to support this melt rate (Reference Owens, Marcus and MooreOwens and others, 1984).
Table II Bulk.-exchange coefficients (K̄) and energy- and mass-balance data for the calibration period (cp)
Determination and Behaviour of K̄
During the CP it was possible to calculate K for 17 periods ranging in duration from 11 to 38 h, but typically of 23–25 h duration. Table II presents estimates of K B and K E for the CP. Determination of the former required a value for z0 using Equation (12). Estimates were made for each hour that conditions were neutral (−0.012 < Rb < 0.012), resulting in only 15 values out of a possible 205. Individual values of z 0 ranged from 0.0 to 0.116 but the median value of 0.014 m was adopted. This is larger than values usually reported in the literature (see Reference MooreMoore, 1983[b]) but is reasonable given the heavily sun-cusped and rain-scoured nature of the Ivory Glacier ablation surface.
For the entire CP, the mean daily values of K B and K E were 0.0037 (range 0.0–0.0075; standard deviation ±0.0020) and 0.0039 (range 0.0004–0.0093; standard deviation ±0.0024), respectively. Despite the similarity in the means and variances of these estimates of K , the correlation coefficient for the two time series was only 0.43. This poor correlation may result from measurement errors and failure to meet the assumptions. Both possible sources have been discussed earlier but, as an example, measurement of ablation using changes in height of surface ice at fixed poles results in a measurement error over one period being associated with a compensating error in the following or subsequent periods. This, at least, partly explains the occurrence of the highest and lowest values of K E on consecutive days (11 and 12 February; Table II). Table III demonstrates how variation in K E over these two days is considerably reduced as a result of averaging over intervals which are increased progressively from one-half, to one, and finally to 2 d. This analysis indicates why in this paper results are presented for averaging periods of 24 h despite the fact that the twice-daily ablation measurements allow K to be estimated and the models validated for shorter time intervals.
In an attempt to explain some of the other variability in K the exchange coefficient was correlated with a variety of wind and stability statistics. Table IV shows that K E has a strong dependency on both the wind speed and the wind gradient, and a weaker dependency on stability. In contrast, K B shows weak dependency on the wind gradient but, as is obvious from Equation (10), it exhibits a strong association with the Richardson number and therefore wind speed.
Table III Values of the bulk-exchange coefficient ( K E) over increasing time intervals
Table IV. Correlation coefficients for relationships between the bulk-exchange coefficients K B and K E, and various wind and stability parameters for 17 daily observations
Given these results, and the fewer assumptions required, it is postulated that K E is the more realistic short-term estimate of the exchange coefficient. The independent data analysed in the next section will substantiate this assertion.
Model Validation
Data gathered in 1972 (the VP) provided an opportunity to make an independent test of the model values of K E obtained during the CP. Ablation rates for the VP were calculated using Equation (1). Estimates of QH and QE were made with Equations (7) and (8), and the mean value of K E 0.0039) as determined in the CP. Values of ablation, measured with an 8-pole array, were converted to their energy equivalents using Equation (2a). The results are presented in Table V and Figure 3. In Figure 3a the least-squares line forced through the origin has a slope of 0.98, a correlation of 0.79, and a standard error of ±9 mm d−1 water equivalent. The longer-term means are in close agreement, thus illustrating that short-term error can be substantially reduced by temporal averaging. Most of the reduction in error occurs with only a doubling of the time interval (Table VI and Fig. 3b).
The VP also provides a second opportunity to assess the method used to estimate K B With z 0 again equal to 0.014 m, the model exhibits large short- and long-term errors (Table VI) when calculations are performed on an hour-by-hour basis. However, use of daily averaged input data effectively eliminates the bias and substantially reduces the short-term error.
Table VI also shows that, despite the use of K B values that vary with wind speed and stability, short-term errors are always larger than those associated with the use of a constant value of K E determined from independent data using the energy-balance approach. This even applies to hourly calculations, confirming the superiority of the energy-balance approach.
Table V. calculated and measured ablation and other energy- (in mj m−2) and mass-balance (m w in mm; m i in cm) data for the validation period (vp)
Sensitivity Studies
The energy balance requires the selection of an appropriate constant value for K E, while the bulk-aerodynamic approach requires an estimate of z0 in order to calculate K B. Figure 4 illustrates the sensitivity of the validation statistics to changes in K E and z0 in the respective approaches. To minimize errors, the energy-balance calculations were performed on an hourly basis while the bulk aerodynamic method used daily averages of the input data (for justification see Table VI).
The systematic error is considerably less sensitive to errors in z0 than in K E. Specification of the roughness length to within ±0.002 m of its “true value” results in systematic errors of only a few per cent. This finding contrasts with the conclusion ofReference Moore Moore (1983[b]) that errors of around ±20% would result, probably because the sensitivity decreases with increasing values of the roughness length. The large aerodynamic roughness that is typical of the melting Ivory Glacier surface therefore facilitates the use of the bulk-aerodynamic approach.
A somewhat surprising result is that an appropriate constant value for the bulk-exchange coefficient K E can actually provide improved estimates of the ablation on individual days in comparison to those derived using a coefficient K B that varies on a daily basis with wind speed and stability. However, departures from the appropriate value are associated with a rapid deterioration in the short-term performance to the extent that any advantage soon disappears.
Summary and Conclusions
A useful way to explain the past behaviour of a glacier is to model its mass balance, rather than consider the roles of temperature and precipitation in isolation. This study shows that an important aspect of the mass balance, namely ablation, can be successfully calculated through the energy balance, with the turbulent fluxes evaluated using a bulk-exchange coefficient ( K E). This model only applies when the glacier is melting. On Ivory Glacier in summer, the mean value of this coefficient is 0.0039. A similar value was obtained using the bulk-aerodynamic approach and a measured value of 0.014 m for z 0. In an independent (different summer, same glacier) validation and, despite the use of a constant exchange coefficient, the energy-balance approach consistently out-performed the bulk-aerodynamic approach on a day-to-day basis.
Fig. 3. Calculated (using the energy-balance approach) and measured values of ablation (Qm and M i) during the validation period (VP), a. One-day totals; b. Two-day totals.
Table VIa. Validation statistics for increasing time intervals with K E = 0.0039
Table VIb. Per cent mean bias (mbe) and root-mean-square (rmse) errors for the energy-balance ( K E = 0.0039) and bulk-aerodynamic (variable K B models applied to data averaged over the specified time intervals
* Refers to line of least-squares best fit between measured and calculated ablation.
Fig. 4. Sensitivity of estimates of validation-period (VP) ablation rates to variations in the exchange coefficient ( K E) for the energy-balance approach (continuous line) and to variations in roughness length (z0) for the bulk-aerodynamic approach (dashed line).
Although, in short-term, K E varies with wind speed and its vertical gradient, and with atmospheric stability, the average value of 0.0039 provides good estimates of glacier ablation over periods of two or more days. Under typical field conditions, the short-term estimates of the exchange coefficient provided by the bulk-aerodynamic approach were found to be inferior to the use of this average value, despite the fact that wind-speed and stability effects are accounted for in such calculations.
We do not know how robust the empirically derived bulk-exchange coefficient will be, but we expect that it varies from place to place, and with the season, because it implicitly accounts for the “average” influence of such factors as advection, stability, and aerodynamic roughness on turbulent exchange. On the other hand, combining the present results with those ofReference Kuhn Kuhn (1979) for Hintereisferner in the Ötzal Alps of Austria suggests that the exchange coefficient may be a relatively conservative variable. Kuhn determined a bulk-transfer coefficient for sensible heat from the residual of the energy balance and from its vertical gradient over the entire ablation period as well as for short periods. The exchange coefficient (K̄) is related to his transfer coefficient (α) by
From data in Reference KuhnKuhn (1979) (α = 1.68 MJ m−2 d−1 K−1 and u 1 = 4 m s−1 at 2 m), K equals 0.0045. Given the assumptions and approximations Kuhn used in deriving his four values of the transfer coefficient, this K is remarkably consistent with the Ivory Glacier value of 0.0039.
Clearly, more measurements are required on other glaciers in New Zealand and elsewhere before the value of 0.0039 can be used with any assurance that it provides short-term estimates that are superior to those obtained using the bulk-aerodynamic approach. The latter invokes more restrictive assumptions, some of which may not apply to an alpine glacier environment. However, this study and the results ofReference Moore and Owens Moore and Owens (1984) show that the bulk-aerodynamic approach can provide useful estimates of the turbulent exchanges, given an appropriate value of the surface roughness. The short-term errors can be reduced by using daily averaged as opposed to hourly meteorological data. Nevertheless, they remain larger than those associated with the use of an appropriate constant value for the bulk-exchange coefficient determined using the energy balance.
Reference KuhnKuhn (1979) noted that it is desirable to relate the factors determining the energy balance of a glacier to the large-scale atmospheric processes, thereby minimizing the number of parameters referring to the glacier surface and relying more on synoptic weather data. The two approaches investigated in the present paper can be used to calculate ablation with relatively simple measurements of radiation, wind speed, temperature, and humidity. If the objective is to model past glacier behaviour, such measurements are normally not available. However, with limitations and some difficulty, the necessary long-term data can often be estimated. For example, equivalent data could be interpolated from historic upper-level weather maps and radio-sonde soundings, or alternatively, nearby valley stations could be lapsed to the higher levels of the glacier Reference Tangborn(Tangborn, 1980[a], [b]). By following these procedures, it may be possible to isolate the specific factors that have influenced past glacier behaviour.
Acknowledgements
The support of the New Zealand Ministry of Works, and in particular the assistance of P. Anderton and T. Chinn, are gratefully acknowledged. The field data are the result of the conscientious efforts of D. Harrowfield, F. Harding (VP), and B. Dickson (CP). D. Moore suggested a number of significant improvements to an early version of this paper. P. Jance drew the diagrams.
