Introduction

High-resolution meteorological and glaciological data are needed to model non-linear glacier–atmosphere interactions (e.g. long-wave radiation) and threshold effects (e.g. phase changes and albedo feed-backs). Non-linear effects are particularly significant in maritime climates where temperatures often fluctuate around the melting point of ice. Accurate modeling of glacier mass balance has been hampered because it is difficult to determine the spatial and temporal variations of meteorological variables. Micro-meteorological measurements have been made on only a few glaciers, and most records extend back less than 30 a and cover only part of the year, Twice-daily radiosonde flights started at many locations in 1948, and here we investigate the feasibility of using these data to estimate glacier ablation. Our study is motivated by the possibility that radiosondes, or equivalent gridded atmospheric data, could be used to model the recent mass-balance history of any glacier.

The study has two stages: (i) extrapolation of meteorological data from the radiosonde to glacier surfaces; (ii) use of radiosonde data to drive physical models of mass and energy exchange. For the first, we compare measurements interpolated from a nearby radiosonde with measurements from Blue Glacier, western Washington, U. S.A. National Weather Service radiosonde flights in the region started in 1948 at Tatoosh Island (about 100 km northwest of Blue Glacier) and moved in 1966 to Quillayute (about 70 km west-northwest of the glacier). Radiosondes measure vertical profiles of pressure, temperature, relative humidity, wind speed and direction. With few exceptions, radiosondes have been launched twice daily since 1958 at 0000 and 1200 UTC, and we interpolate each variable linearly in time between flights. In this preliminary investigation, we interpolate variables linearly in altitude from measurements at just three standard pressure levels (850, 700 and 500 mbar). The 850 mbar pressure level is about 150 m higher than the terminus of Blue Glacier, and 700 mbar is about 600 m above the lop of the glacier. In the second part of the study we use the radiosonde data in established formulations of mass and energy exchange at glacier surfaces.

Extrapolating Radiosonde Data to the Glacier

Relationships between radiosonde and ground-based measurements are studied by comparing data from a measurement epoch during the summer of 1991. An automatic weather station, located on rocks at 2008 m beside the glacier, measured air temperature, the 3 m wind speed and solar radiation at 5 min intervals and recorded the average hourly value of each variable.

Figure 1 shows temperatures at the glacier and temperatures interpolated from the radiosonde at 2008 m. The radiosonde does not capture all details of the diurnal cycle, but the free-air temperature is closely similar to the average temperature at the glacier (rms = 2.3°C, r ² = 0.73). By comparison, measurements from a nearby low-altitude station at Quillayute (64 m) offer a poor estimate of glacier temperature (rms = 4.5°C, r ² = 0.19), partly because stratus clouds with lops just below the glacier often complicate the temperature profile in the lower atmosphere during summer.

Fig. 1. Temperatures at the glacier (averaged over the 12 h centered around the radiosonde times), and temperatures interpolated from the radiosonde to 2008 m. On the line, the temperatures are equal.

A realistic simulation of the local winds requires a three-dimensional physical model of interactions between the synoptic winds and topography. However, here we assume the speed at any time is constant over the entire glacier. It is calculated empirically from the speed V ₈₅₀ at 850 mbar (Fig. 2) with an rms of 1.3 m s⁻¹ (r ² = 0.44).

Fig. 2. Wind speeds at the glacier (averaged over the 12 h centered around the radiosonde times), and speeds at 850 mbar from the radiosonde. The line is the best fit to the measurements; on the line the speed at the glacier is 0,34V ₈₅₀ + 1.2.

A model of atmospheric transmissivity is needed to determine the flux of incoming solar radiation. The flux at the top of the atmosphere Q ₀ is calculated from the solar irradiance (about 1371 Wm⁻²; Reference Hickey, Alton, Kyle and Hoyt.Hickey and others, 1988) and ephemeris data (Reference ListList, 1966) which give the Earth–Sun distance through the year. The relative path length m and the angle of incidence Ψ the solar beam makes with the earth surface are simple functions of the surface geometry and the elevation angle θ and azimuth of the Sun. Here m is the ratio of the actual path length to the zenith path length. Values tor these variables are calculated as a function of the hour of day and day of year (Reference RasmussenRasmussen, 1974). The solar beam is attenuated by air molecules, aerosols, water vapor and clouds as it passes through the atmosphere, and we assume the transmissivity τ _SW can be factored into a clear-air component τ _air and a cloudy component τ _Cl.

(1)

The clear-air component is further factored into a dry-air component τ _dry (which includes the effects of aerosols) and a water-vapor component τ _w.

(2)

The vapor component (Reference Lacis and Hansen.Lacis and Hansen, 1974) is

(3)

in which w is the precipitable water (in mm) above the site. The dry-air factor is a function of p, the pressure at altitude relative to the pressure at sea level (Reference ListList, 1966).

(4)

The optical properties of clouds vary with cloud type and drop-size distribution, but here we assume the cloudy component τ _cl is a simple, linear fonction of the relative humidity H (in per cent) at the indicated pressure levels.

(5)

The five model parameters in Equations (4) and (5) were determined by optimizing the agreement between shortwave radiation modeled using Equations (1) and (7), and hourly measurements from the weather station.

A fully saturated atmosphere (H = 100% at all levels) implies τ _cl = 0. This condition is rare; a typical “wet” condition (H ₈₅₀ = 92%; H ₇₀₀ = 91%; H ₅₀₀ = 66%) yields τ _cl = 0.194. When H = 55% at all levels, τ _cl = 1; but τ _cl is bounded in Equation (1) so it never exceeds 1. The transmissivity τ _SW is continuous, but the slope of the function changes abruptly when the attenuation from the cloudy component sets in, at about 55%. In reality, an abrupt increase in short-wave attenuation occurs when atmospheric conditions cause clouds to form.

We ignore the effects of reflecting surfaces but we do partition the incoming beam into a direct-beam component Q _dir and a diffuse-sky component Q _dif , because the direct-beam component received by a surface depends on the surface slope. We modify a formulation suggested by Reference OerlemansOerlemans (1992),

(6)

where f = (0.2 + 0.65 τ _sw), which implies that the direct-beam component increases linearly with transmissivity to a maximum of 85% of the total incoming radiation when τ _sw = 1. We have no data to support this formulation because our radiation measurements were made on a horizontal surface and do not contain information on the partitioning. That is, when Ψ = θ,

(7)

Figure 3 shows that when the transmissivity model is used to estimate incoming short-wave radiation at 2008 m, the agreement with glacier measurements has an rms of 29 W m⁻¹ (r ² =0.86). The agreement is important because short-wave radiation is the dominant source of energy to the summer snow surface. Furthermore, the method is valuable because it eliminates the difficulty of obtaining- an explicit measure of cloudiness.

Fig. 3. Modeled and measured flux of short-wave radiation at 2008 m. The flux is estimated using the transmissivity model which makes use of the temperature and humidity profiles from the radiosonde. On the line, the modeled and measured fluxes are equal.

The automatic weather station did not measure relative humidity. However, wet- and dry-bulb temperatures, measured manually at 1600 UTC during the summers of 1985–87 at 2008 m. were used to calculate relative humidity. The agreement between glacier data and values interpolated in time and altitude from the radiosonde (Fig. 4) has an rms of 20%, (r ² = 0.44). Measurements at the glacier are about 10% higher than from the radiosonde, which might be expected because melting snow or ice is a potential source of moisture; but we have some reservations about the accuracy of the measurements. For example, the large number of days when H = 100% on the glacier (about 15% of all days) makes us suspect that the wet-bulb thermometer was not always fully saturated. This problem would result in artificially high values of H at the glacier. A different type of error is introduced because radiosonde sensors are not accurate at low humidities. Most radiosonde data are post-processed and any low H values are increased to 19% (Reference Elliot and Gaffen.Elliot and Gaffen, 1991), which effectively increases the average relative humidity slightly.

Fig. 4. Measurements of relative humidity made at 1600 UTC during the summers of 1985–87 at the glacier, and values interpolated in altitude and in time from the radiosonde. On the line, the relative humidities are equal.

Mass and Energy Exchange at the Glacier Surface

The energy exchanged at the surface of temperate glaciers such as Blue Glacier is much greater than the energy associated with ice flow or geothermal fluxes; here we ignore contributions from these sources. We also ignore any heat storage on the assumption that there is sufficient winter rainfall to provide latent heat to keep the temperature of the glacier at 0°C. For a melting surface, the energy exchange and ablation are coupled:

(8)

where dM/dt is the ablation rate (in m s⁻¹ w.e.), ρ _w the density of water in kg m⁻³, L _M the latent heat of melting (0.335 MJ kg⁻¹) and Q _SW, Q _LW, Q _SH and Q _LH are the net fluxes (in W m⁻²) of short-wave (solar) radiation, long-wave (terrestrial) radiation, sensible heat and latent heat, respectively, Fluxes directed towards the surface have a positive sign and those leaving are negative. We use meteorological measurements interpolated in time and altitude from radiosondes to derive physical models of mass and energy exchange at the glacier surface. Model parameters are derived by optimizing the fit to measurements of net short-wave and long-wave radiation, mass condensation, evaporation and melt that were made at 2050 m on Blue Glacier during 1958 (Reference LaChapelleLaChapelle, 1959a, Reference LaChapelle1960). The sensible-heat flux was estimated by calculating the residual of these measurements. The various models and optimization procedures are discussed below.

(i) Absorbed Short-Wave Radiation

The absorbed short-wave radiation depends on the effective albedo α:

(9)

where Q _dir and Q _dif are defined in Equations (6). Albedo varies as the surface properties evolve during the ablation season. Here we assume that aging processes (e.g. grain coarsening, accumulation of light-absorbing particulates, surface roughness) cause the surface albedo α′ to decrease exponentially with time from an initial value α ₀ at the start of the ablation season to a limiting value α _n . The formulation is similar to one proposed by Reference Oerlemans and Hoogendoorn.Oerlemans and Hoogendoorn (1989):

with a solution,

(10)

where t _e is the e-folding time-scale in days. Values for α ₀ and α _n , which were derived from measurements made by LaChapelle in 1958, differ for snow (0.85 and 0.59), firn (0.48 and 0.40) and ice (0.40 and 0.33). The time-scale t _e is taken to be 30 d for all three surfaces.

Clouds absorb light at near-infrared wavelengths, which effectively increases the spectrally integrated albedo at the surface (Reference WarrenWarren, 1982). In the absence of detailed information about clouds, we model the wavelength dependence of albedo implicitly and assume it varies with the attenuation due to clouds. The effective albedo is written as a linear combination of 1 – τ _cl and α′:

(11)

The two coefficients for the parameterization were obtained by optimizing the fit to daily measurements of snow albedo during 1958 (Reference LaChapelleLaChapelle, 1960). Figure 5 shows the model agrees very well with measurements (rms = 0.03, r ² = 0.89).

Fig. 5. Modeled and measured surface albedo at 2050 m dining 1958. The model parameterizes the effects of clouds and surface aging processes (Equation (11)). During the period of measurements, the surface remained snow-covered (α ₀ = 0.85, α _n = 0.59, t _e = 30 d).

(ii) Absorbed Long-Wave Radiation

The components of incoming long-wave radiation are from the clear-sky (Q _air) and from clouds (Q _cl). The emissivity of snow ε _snow is 0.99 (Reference WarrenWarren, 1982), so the resulting outgoing black-body radiation from melting-snow or ice is 312.5 Wm⁻². The net long-wave flux absorbed by the surface is

(12)

The clear-sky component is

(13)

The Stefan–Boltzmann constant σ is 5.67 × 10⁻⁸ W m⁻²K⁻⁴, T _air is the temperature (K) and ε _cs is the emissivity of the clear sky, which is calculated as a function of altitude z (meters) and the vapor pressure (Reference Kimball, Idso and Aase.Kimball and others, 1982; Reference OerlemansOerlemans, 1992).

(14)

The vapor pressure P _vap (Pa) is estimated from the relative humidity and an approximation for the saturation vapor pressure given by Reference KrausKraus (1972).

Q _cl depends on cloud type, height, temperature and cloud amount; but, lacking explicit measurements of cloudiness, we assume the flux varies with the attenuation due to clouds:

(15)

The cloud temperature T _cl is taken to be the temperature 1000 m above the site. Both the altitude coefficient in Equation (14) and the coefficient in Equation (15) were obtained by optimizing the agreement between Equation (12) and Reference LaChapelleLaChapelle’s (1959a) daily measurements, which has an rms of 12 W m⁻² (r ² = 0.43).

(iii) Sensible Heat

The flux of sensible heat ran be written (Reference DeardorffDeardorff, 1968)

(16)

in which v is the wind speed, and ΔT _sfc is the temperature difference between the air and the surface. A value for the mean bulk turbulent coefficient D _T = 1.65 J m⁻³ K⁻¹ was derived using concurrent measurements of wind speed and air temperature, and Reference LaChapelleLaChapelle’s (1959a) calculations of the daily average flux. LaChapelle’s values are questionable because they were calculated as the residual of the energy equation, and any errors in the measurements of the other fluxes or of ablation would accumulate in this term. This is particularly evident on occasions when the calculated flux was negative and yet ΔΤ _sfc. was positive. Other errors arise because of uncertainties in the wind speed; also, the assumptions implicit in Equation (16) may not be adequate (Reference MunroMunro, 1989). Given the difficulties, it is not surprising that the uncertainty in the estimate of the daily flux is high (rms = 31 W m⁻², r ² = 0.49).

(iv) Latent Heat

The flux of latent heat over a melting snow or ice surface can be written (Reference DeardorffDeardorff, 1968)

(17)

in which L _E is the latent heat of vaporization (2.51 MJ kg⁻¹), c _p the specific heat of air, v the wind speed, P _atm the atmospheric pressure. P _vap the vapor pressure of the air, and P _sat the saturation vapor pressure over a melting ice surface. In 1958, LaChapelle calculated the magnitude of the latent-heat flux directly from measurements of mass condensation or evaporation. The agreement between the model and measurements has an rms of 11 W m⁻² (r ² = 0.30).

Results of Ablation Model

We are now in a position to model ablation at any site on the glacier. The model is run using hourly time steps, primarily to resolve the daily solar cycle, fan also to represent temporal variations of meteorological variables through linear interpolation. Air temperature and relative humidity are interpolated directly in time and altitude from radiosondes. Wind speed is assumed to be a simple function of the speed at 850 mbar (see Fig. 2) and spatially invariant over the glacier. Although the influences of local topography and irregular drainage winds complicate the local wind field, we do not have data from other sites to support a more sophisticated formulation. Altitudinal variations of the clear-air transmissivity are modeled explicitly, but the cloud factor τ _cl is assumed to be constant over the altitude range of the glacier. This implies that the cloud base is always above the top of the glacier, which is not always true. The influences of topography on the magnitude and partitioning of short-wave radiation at each site are calculated using a 30 m digital elevation model. The parameters for the albedo model, α ₀ and α _n , are changed as the surface evolves from snow, to firn, to ice.

Stake networks have been used to measure summer ablation at Blue Glacier for most years since 1958. Variations in the near-surface density and sub-surface melting can introduce errors of up to 20% in stake measurements (Reference LaChapelleLaChapelle, 1959b). Here we convert measurements of surface lowering to water equivalent by using ρ _ice = 900 kg m⁻³; ρ _firn = 580 kg m⁻³; ρ _snow = 510 kg m⁻³, based on measurements of average density from snow pits at Blue Glacier. This transformation could introduce systematic errors, particularly when freshly deposited snow of low density exists at the glacier surface.

We test the ablation model using measurements from one site at 2050 m in the accumulation zone and another at 1540 m in the ablation zone. Figure 6 indicates the model does well over a wide range of altitudes, and the agreement (less than 16%) is within the errors of stake measurements. On average, there is a positive bias of 1.9 mm d⁻¹; that is, modeled ablation is higher than measurements. Differences in albedo (ice is exposed for longer periods at lower altitudes) and generally higher temperatures at lower altitudes cause higher ablation there; the ablation at 2050 m is about that at 1540 m (Fig. 6).

Fig. 6. Modeled and measured ablation in the accumulation zone (at 2050 m) and the ablation zone (at 1540 m) showing 53 seasonal measurements from 1958 to 1987. On the line, the modeled and the measured ablation ate equal. Ash from Mount St. Helens was deposited on the glacier in May 1980 .

The 1980 ablation season was anomalous because a thin layer of volcanic ash from Mount St. Helens (about 220 km south-southeast of Blue Glacier) remained at the surface throughout the season. The measured ablation at 2050 m (33 mm d⁻¹) is much higher than the model (24 mm d⁻¹) when using the albedo formulation from Equation (11). This suggests the ash lowered the albedo significantly at that site. The impact is less evident at 1540 m, probably because the albedo of ice does not change much when covered by ash.

Conclusions

In a maritime climate such as western Washington’s, atmospheric measurements from radiosondes provide a much better estimate of mountain conditions than do measurements from low-altitude stations. A method of using radiosonde data to model mass and energy exchange at glacier surfaces was tested at Blue Glacier. The model was not tuned explicitly to measurements of ablation; rather, model parameters were adjusted to optimize relationships between meteorological variables and measurements of the surface energy fluxes. Nevertheless, the modeled ablation is within the accuracy of stake measurements, which suggests the meteorological data from radiosondes are suitable for estimating seasonal ablation on Blue Glacier. The slight positive bias in the modeled ablation indicates that the model parameters (e.g. the coefficients in Equations (14)–(17)) are not optimally tuned. The results are highly sensitive to the numerical values of many model parameters, but which values need revising is not clear in the absence of additional information. Another possible improvement would be to make better use of the vertical resolution of the radiosondes.

The physical basis of the model gives us confidence that the methodology should apply to any glacier. Radiosonde or equivalent gridded data with wide spatial coverage extend back more than 45 a and offer an opportunity to model the recent ablation history of glaciers, including those where micro-meteorological measurements are not available. Furthermore, this type of approach could easily be coupled to a predictive regional climate model and used to assess the response of glaciers to possible future climate scenarios.