Near-surface internal melting: a substantial mass loss on Antarctic Dry Valley glaciers

Abstract

The McMurdo Dry Valleys, southern Victoria Land, East Antarctica, are a polar desert, and melt from glacial ice is the primary source of water to streams, lakes and associated ecosystems. Previous work found that to adequately model glacier ablation and subsurface ice temperatures with a surface energy-balance model required including the transmission of solar radiation into the ice. Here we investigate the contribution of subsurface melt to the mass balance of (and runoff from) Dry Valley glaciers by including a drainage process in the model and applying the model to three glacier sites using 13 years of hourly meteorological data. Model results for the smooth glacier surfaces common to many glaciers in the Dry Valleys showed that sublimation was typically the largest component of surface lowering, with rare episodes of surface melting, consistent with anecdotal field observations. Results also showed extensive internal melting 5–15 cm below the ice surface, the drainage of which accounted for ~50% of summer ablation. This is consistent with field observations of subsurface streams and formation of a weathering crust. We identify an annual cycle of weathering crust formation in summer and its removal during the 10 months of winter sublimation.

1. Introduction and Background

The 4500 km² (Levy, 2013) McMurdo Dry Valleys (MCM) is the largest ice-free area in Antarctica (Drewry and others, 1982). However, the primary water source to aquatic ecosystems in the streams and lakes of the valleys is runoff from small alpine glaciers and outlet glaciers of the East Antarctic ice sheet (Fig. 1) that flow into the valley bottoms (Fountain and others, 1999a). The motivation of this study was to develop a predictive model of glacial melt to support biological research investigating ecological structure and function in the polar desert of the MCM (Priscu, 1999 and references therein). Understanding meltwater production processes and variations in meltwater flux is particularly intriguing and challenging because the MCM summer climate is commonly a few degrees below freezing and the near-constant winds keep glacial surfaces frozen despite continuous solar radiation (Lewis and others, 1998; Hoffman and others, 2008). Under these conditions, small changes in the energy balance can result in profound changes in ablation and meltwater flux (Fountain and others, 1998; Hoffman and others, 2008). The purpose of this report is twofold: to explain the physical controls on meltwater production in a polar environment and to better elucidate near-surface glacier melting processes. The MCM environment provides an ideal natural laboratory for investigating melt processes that are otherwise overwhelmed in the high-energy and high-melt environments typical of most temperate glaciers.

Fig. 1. Map of Taylor Valley showing locations of the TAR, TAR2 and CAA sites used in this study. Glacier stakes used in the mass-balance calculation for each study site are indicated with black dots.

In environments where the air temperature hovers near the melting temperature of ice, the melting behavior becomes complex. Solar radiation penetrates to substantial depths in ice if the surface and internal scattering are low. When the ice surface is cooled by longwave radiation loss and a colder atmospheric boundary layer, penetrating solar radiation can generate a ‘solid-state greenhouse effect’ resulting in a maximum in the temperature profile at some depth beneath the surface (Brandt and Warren, 1993; Liston and others, 1999). Such subsurface heating can drive ice ‘weathering’ or ‘rotting’ through internal melting and partial drainage, lowering ice bulk density (Müller and Keeler, 1969; Larson, 1977; Braithwaite and Olesen, 1990; Hubbard and Glasser, 2005).

Subsurface ablation will only occur if meltwater is able to drain; refreezing in place results in no mass loss (Fig. 2). Because ice melts preferentially along grain boundaries (Müller and Keeler, 1969; Raymond and Harrison, 1975; Nye, 1991), a drainage network of channels can develop at low volumes of internal melt. Observations suggest inter-granular water flux is typically small (<0.1 m a^–1) prior to the initiation of weathering, particularly when bubbles are present that can block the flow through veins (Raymond and Harrison, 1975). However, once channels become large enough to overcome blocking bubbles and capillary tension, water can drain. Ice affected by internal melting appears to act like an unconfined aquifer, with surface melt percolating down to the water table (Wakahama and others, 1973; Larson, 1977). On temperate glaciers the water flux through the near-surface layer is negligible compared to that of supraglacial streams (Fountain and Walder, 1998). In contrast, near-surface drainage may be important on areas of MCM glaciers lacking well-developed supraglacial streams. For glaciers in the MCM, a near-surface drainage system has been observed in refrozen fractures that contain blue ice and melt preferentially (Fountain and others, 2004). However, the spatial extent, temporal evolution and connectivity within this system are unclear, as is its behavior in the bubbly, white ice found in these glaciers.

Fig. 2. Schematic of how drainage of subsurface melt can reduce ice density, causing subsequent ablation at the surface (sublimation or melt) to remove a greater thickness of ice. (a) shows the case of all subsurface melt refreezing in place, and (b) shows the case of some subsurface melt draining. The pole and flag represent an ablation stake, and the arrows represent ablation at the surface due to sublimation and/or melt. Features are not to scale.

On most temperate glaciers, the contribution of ice weathering processes to glacier surface mass balance is generally ignored since the effects are small relative to ablation at the surface (Müller and Keeler, 1969; Braithwaite and Olesen, 1990). However, the resulting ice density variations may be significant on polar glaciers where summer ablation is small. Müller and Keeler, (1969) observed the melting and draining of 1.3 cm w.e. of subsur-face ice over 12 hours in the upper 10 cm of White Glacier, Canadian Arctic. For the glaciers of the MCM, where summer ablation is <15 cm w.e. (Fountain and others, 2006), the drainage of 1.3 cm w.e. of subsurface melt over an entire summer yields an underestimate of ~10% in ablation measured by the glaciological method (the difference in the glacier surface level between two dates, typically measured on a stake drilled into the glacier; Fountain and others, 1999b and references therein; Kaser and others, 2003), assuming a constant ice density. On these glaciers, subsurface melt is observed in cryoconite holes and in the blue ice of refrozen cracks (Fountain and others, 2004). Additionally, modeling predicts subsurface melt in the bubbly, white ice that makes up the bulk of the glacier surfaces (Hoffman and others, 2008), suggesting ice weathering associated with internal melting may be an important ablative process on these glaciers.

In a previous study, we explored the year-round surface energy balance and surface melt production at daily time-steps for one site on Taylor Glacier in the MCM (Hoffman and others, 2008). Results showed that it was necessary to model penetration of solar radiation into the ice to accurately predict ablation. Here we investigate the mass loss due to drainage of internal melt, a process not included in the previous study, and which alters the subsurface heat balance. In addition, the model is applied to three sites on two different glaciers within the MCM to test the model across a wider range of surface conditions. Finally, we also assess parameter sensitivity through a detailed calibration procedure. We found that including the drainage of subsur-face melt improved model skill in predicting surface lowering, ice density and near-surface ice temperatures and that near-surface internal melting leads to substantial mass loss on MCM glaciers.

2. Site Description and Observations

Mean annual air temperature in the MCM ranges from –14°C to –30°C (Doran and others, 2002), but mean summer air temperatures are near 0°C in the valley bottoms (Nylen and others, 2004), allowing for an average 10 week glacier melt season (McKnight and others, 1999). Precipitation falls as snow, is <50 mm w.e. a^–1 and shows little seasonality (Keys, 1980; Fountain and others, 2010). Annual net mass loss in the ablation zones of the glaciers is on the order of 10^–1 m, with 70–90% in the form of sublimation, 15–30% in the form of melt, and 1–3% in the form of calving from terminal cliffs (Fountain and others, 1998). Runoff only occurs in the generally snow-free ablation zones; little melt occurs in the snow of the accumulation zones and any that does refreezes at depth (Fountain and others, 1998). The equilibrium-line altitudes of the glaciers rise rapidly up-valley by 30 m km^–1, due to gradients in precipitation and wind speed within the valleys (Fountain and others, 1999c). The alpine glaciers of the MCM are polar glaciers with interior and basal temperatures much lower than melting, and are frozen to their beds (Fountain and others, 1998; Cuffey and others, 2000; Fitzsimons and others, 2008). The cold ice prevents the development of englacial or subglacial passages for conveying meltwater, and runoff is restricted to the glacier surfaces (Fountain and others, 1998).

Within Taylor Valley, the McMurdo Dry Valley Long Term Ecological Research project (mcmlter.org) has maintained meteorological stations (Doran and others, 1995) and mass-balance stake networks (Fountain and others, 2006) on four glaciers continuously since 1994. Continuous measurements are collected of air temperature, relative humidity, wind speed (all at 3 m), incoming and outgoing shortwave (solar) radiation, incoming longwave (thermal) radiation, subsurface ice temperature, and distance to the ice surface for ablation and snowfall (Table 1). Most of the sensors collect data every 30 s (every 4 s for wind speed), and 15 min averages are stored on a solid-state data logger. Ice ablation causes the instrument height to increase over time and the station is periodically reset. Further details on these observations are presented in Hoffman and others, (2008) and Hoffman, (2011), and at http://www.mcmlter.org/data_home.htm. Mass-balance observations are collected in early November and late January using the glaciological method to provide seasonal estimates of mass change (Fountain and others, 2006). Because these observations ignore internal melt, we use the original measurements of surface lowering rather than calculations of ablation.

Table 1. Summary of study sites. Climate and mass-balance statistics are for the period December 2004 through January 2009. Summer is defined as December and January, and winter is the remainder of the year. Meteorological stations are located at TAR and CAA whereas the weather at TAR2 is calculated. n/a: information not available

To investigate ice density variations due to internal melting, six shallow (~50–100 cm) ice cores from Taylor, Canada, and Commonwealth glaciers were collected during the summers of 2007/08 and 2008/09. Each core was cut into segments between 4 and 43 cm long. The mass of each segment was measured, and the volume was calculated from measurements of length and diameter, to yield density variations with depth similar conceptually to the procedure described for snow by LaChapelle, (1959). We calculated uncertainty in ice density by propagating the standard error in three measurement replicates of the length, diameter and mass of each ice-core section through the density calculation. Ice below 40 cm depth had a uniform, unaltered appearance, and density was calculated to be 870 ± 30 kg m^–3 (Hoffman, 2011). Ice density results are presented in Section 5.3.

Three sites on glaciers in the MCM were studied including two sites coincident with meterological stations and subsurface ice temperature measurements (Table 1; Fig. 1): Taylor Glacier at the TAR station and Canada Glacier at the CAA station. Because ablation measured with stakes is subject to local variations (Braithwaite and others, 1998), we use the mean ablation of three to five stakes within 2.5 km of each study site which have a standard deviation of 2.1 cm. To consider a higher-ablation environment that is common in the low-elevation regions near the termini of MCM glaciers, an additional site 2 km away and 136 m lower near the terminus of Taylor Glacier (TAR2) is included. That site has no meteorological station, but conditions were calculated using the MicroMet model (Liston and Elder, 2006) for spatial interpolation between stations (detailed in Hoffman, 2011) and the measured albedo at the TAR station. The ice at all three sites is debris-free and representative of the clean, flat surfaces that cover the majority of these glaciers; we do not consider here the influence of debris on melt and energy balance of the ice, which can increase melt rates locally on the MCM glaciers by orders of magnitude (Fountain and others, 2004; Johnston and others, 2005; MacDonell and others, 2013).

3. Model Description

We apply the one-dimensional model based on Liston and others, (1999) and modified by Hoffman and others, (2008) because it accounts for both the surface energy balance and subsurface solar heating from absorbed solar radiation. The surface energy balance takes the form

(1)

where χ allocates the total solar radiation between the surface and that which penetrates into the ice acting as a local heat source, α is albedo, Q _si is incoming shortwave solar radiation, Q _li is incoming longwave radiation, Q _le is emitted thermal radiation, Q _h is sensible heat flux, Q _e is latent heat flux, Q _c is heat conduction in the ice, and Q _m is the energy available for melt, calculated as a residual. All heat flux terms (W m^–2) are positive toward the surface. Q _si, α and Q _li are measured directly, and the terms that cannot be directly measured (Q _le, Q _h, Q _e, Q _c) are cast in a form that leaves surface temperature, T ₀, as the only unknown (Liston and others, 1999). T ₀ is solved iteratively using Eqn (1), and if its value exceeds 0ºC, T ₀ is reset to 0ºC, and the temperature difference is used to calculate melt energy and the corresponding volume of meltwater produced.

The turbulent terms of sensible and latent heat, Q _h and Q _e, respectively, are estimated using a bulk energy method with a correction for stability based on Monin–Obukhov theory (Brutsaert, 1982):

(2)

(3)

where ρ _a is the density of air, c_p is the specific heat of ice, L _v is the latent heat of fusion of ice, and ζ is a stability function based on the Richardson number. T denotes temperature, e denotes vapor pressure and p denotes pressure, with subscript r indicating conditions at a reference height (3 m at TAR, 2 m at CAA) and subscript 0 indicating conditions at the ice surface. The turbulent exchange coefficients are a function of the surface roughness length, z ₀, which is the only adjustable parameter in the turbulent flux formulation:

(4)

where κ is von Kármán’s constant (0.4) and u _r is the wind speed at the reference height, z _r.

The conductive heat flux, Q _c, is calculated by a one-dimensional heat-transfer equation,

(5)

where T _i (K) is the ice temperature, z (m) is the vertical coordinate, t (s) is time,_i is the density of glacier ice (kg m^–3), C_p (J kg^–1 K^–1) is the specific heat of the ice, and q (W m^–2) is the net solar radiative flux. The thermal conductivity of the ice is dependent on the water fraction and ice temperature and density (Hoffman and others, 2008).

The spectrally dependent solar-radiation source term, ∂q/∂z, is computed using a two-stream approximation (Schlatter, 1972) for the radiation flux penetrating the ice, q. The absorption and reflection coefficients in the two-stream approximation are determined from the solar spectrum and the spectral-flux extinction coefficient, ηλ:

(6)

where ρ _ice is the density of pure ice and r _eff is the effective ice grain radius (Brandt and Warren, 1993; Liston and others, 1999). The extinction efficiency, Q _ext, the single-scattering albedo, 1 –ω, and the asymmetry factor, g, are all a function of wavelength and r _eff (Wiscombe and Warren, 1980). The net effect is for ηλ to vary approximately inversely with r _eff. Warren and Brandt, (2008) presented revised estimates of the imaginary part of the complex refractive index for ice that are substantially lower than previous estimates (Warren, 1984) for wavelengths between 300 and 500 nm. However, recalculating Q _ext, (1 – ω) and g based on this update requires Mie scattering algorithms beyond the scope of this study. We estimated that the revisions do not change the order of magnitude of the solar-radiation source term, and we mitigated errors associated with using the older optical constants by calibrating our model to observations of subsurface ice temperature (Section 4.1).

The surface energy balance Eqn (1) represents the upper boundary condition for Eqn (5), and the bottom boundary condition is no heat flux at depth 15 m, approximately the depth of attenuation of the seasonal temperature variations and where the ice temperature equals the mean annual temperature (Paterson, 1994). Equation (5) is applied on a vertical grid with 1 cm resolution in the upper 50 cm and exponentially increasing layer thickness below that for a total of 170 layers. The thickness of the upper layer of the vertical grid is adjusted to change . Ice is allowed to melt when energy is input beyond that needed to raise the temperature to 08C. The model is run with an hourly time-step and forced by hourly averaged meteorological data and daily averaged albedo (Hoffman, 2011).

3.1. Treatment of drainage of subsurface melt

The original model formulation assumed the water fraction within the subsurface ice was in steady state: any water that drained out of a gridcell was replaced by water flowing in (Liston and others, 1999; Hoffman and others, 2008). In this application, we investigate drainage of subsurface melt by comparing three increasingly sophisticated versions of the model. For reference to more traditional surface energy-balance models, we first consider an ‘Opaque Ice Model’ that does not allow the penetration of solar radiation into the ice column (dq/dz = 0 in Eqn (5) and χ = 1.0 in Eqn (1)). In the ‘Refreezing Model’, penetration of solar radiation into the ice is modeled and all subsurface melt refreezes in place (Fig. 2a). In the ‘Drainage Model’ a portion of the modeled subsurface melt is removed from the ice column (Fig. 2b); any subsurface melt fraction that exceeds a 10% threshold in each time-step is assumed to drain and water fraction is set back to 10%. We do not account for an air fraction within the gridcell, and therefore ice density (and associated variations in ice thermal conductivity) does not change during model runs, but we do account for time-evolving density due to drainage when calculating surface lowering. We choose a threshold of 10% water fraction for drainage to be consistent with observations of the irreducible water fraction (the volume fraction of the media occupied by water held in capillary tension and unable to drain) in snow of 5–15% (Coléou and Lesaffre, 1998) and observations of seaice permeability reducing to 0 when porosity is less than 5–10% due to capillary tension (Golden and others, 2007). In reality, the drainage of subsurface melt is likely to be dependent on slope and the extent of cavity formation (from previous drainage). Lacking observations with which to quantify these effects, we impose a constant drainage threshold in the model. In the Drainage Model all drained water is assumed to instantaneously become runoff, ignoring variations in permeability that are expected to occur with changes in porosity. All surface melt is assumed to become runoff for all three models.

3.2. Treatment of snow

Because snowfall in Taylor Valley is infrequent, accumulation is small (<10 cm) and snow presence is of relatively short duration, of the order of days (Fountain and others, 2010), we treat snow simply. Snow cover is estimated from an albedo threshold of 0.10 above the observed daily average bare ice albedo which is dependent on day of year due to dependence on solar elevation angle (Hoffman, 2011). This threshold was determined by comparing sonic ranger observations of snowfall with albedo variations. When albedo cannot be measured because the sun is below the horizon (~April–September), snow cover is ignored. Although this may be a limitation, our goal is to predict melt, and no melt occurs during the dark winter days. Furthermore, snow typically ablates through sublimation or wind erosion after days to weeks following initial accumulation (Fountain and others, 2010). During the sunlit hours of the year, ice ablation is assumed zero in the presence of a snow cover.

The ablation of snow is not modeled, but we instead rely on the albedo record to indicate when bare ice is once again exposed. Thus, the model does not allow for the possibility of snowmelt, which is of minor importance in the Dry Valleys (Chinn, 1981; Fountain and others, 1999a; Gooseff and others, 2003, 2006). This does not affect calculations of melt because the high albedo of snow and the sub-freezing air temperatures common in summer preclude melting in all but extreme and infrequent conditions. For simplicity, thermal and optical properties and surface roughness are left unchanged in the model when snow is present, but net solar radiation is reduced due to the increased albedo. Snow insulates the ice below from temperature changes at the surface, but also reduces the incoming solar radiation, and the effect of snow cover on subsurface ice temperatures depends on the relative magnitude of these processes (Hoffman, 2011). The influence of this simple snow accounting (or not accounting) on the simulations and results presented herein is thought to be negligible.

4. Model Calibration

The model results are sensitive to three parameters: surface roughness, z ₀, effective ice grain radius, r _eff, and solar radiation surface fraction, χ (Table 2). Both z ₀ and r _eff have physical significance and can be constrained by values from the literature, whereas χ is a numerical consequence of a discrete approximation of a continuous process. The surface roughness modifies the magnitude of the turbulent exchange of heat and mass (Eqn (4)), and therefore controls surface melt and sublimation. The magnitude of the turbulent fluxes more than doubles with an order-of-magnitude increase in z ₀ (Brock and others, 2000), and z ₀ varies over nearly five orders of magnitude for glacier surfaces (Brock and others, 2006; Smeets and Van den Broeke, 2008). The effective ice grain radius affects internal scattering of solar radiation and consequently affects both internal ice temperature and melt.

Table 2. Summary of model parameters used for calibration. Numbers in parentheses indicate number of values sampled for each parameter. For determining , 14 values were used for the thickness of the surface layer. T =subsurface ice temperatures, L =summer surface lowering

The inclusion of χ is required to satisfy conservation of energy when radiative processes are included both at the surface (Eqn (1)) and internally as a subsurface source term (Eqn (5)). Its value is a function of the effective ice grain radius and a specified surface layer thickness, within which no solar radiation is absorbed. The surface layer thickness can be thought of as the thickness of ice beneath the air/ice interface that interacts with the atmosphere over the duration of the model time-step (Hoffman and others, 2008). Due to the greater energy required to sublimate vs melt ice (~8 times greater) and the fact that solar radiation is the dominant energy source in summer and a large contributor of energy to melt events (Paterson, 1994; Hock, 2005; Hoffman and others, 2008), the value of χ strongly affects modeled melt.

We first used measurements of ice temperature to constrain a likely range of r _eff, and then used measurements of surface lowering to find optimal combinations of all three parameters. The calibration process considers results from 1680 simulations using both the Refreezing and Drainage models and eight simulations of the Opaque Ice Model at each of the three study locations for the period July 1996– January 2009.

4.1. Ice temperature

Along with ice density and albedo, the effective ice grain radius (r _eff) determines the distribution of absorbed solar radiation with depth through the extinction coefficient (Eqn (6)) (Brandt and Warren, 1993; Liston and others, 1999). Density and albedo are well constrained through direct measurements, but effective ice grain radius cannot be measured directly because it includes the optical effects of bubbles and impurities in addition to the size of physical grains (Brandt and Warren, 1993; Liston and others, 1999; Hoffman and others, 2008). Glacier ice in the Dry Valleys has high bubble content, which causes r _eff to be substantially lower than the measured grain radius (Mullen and Warren, 1988).

The Refreezing and Drainage models were run with 25 different values for r _eff varying between 0.005 and 1.0 mm. The effective ice grain radius was calibrated to the mean absolute error of summer ice temperatures in the upper 1 m of ice. Temperatures used in the calibration were measured by thermistors frozen within the ice for periods of 12–50 contiguous days during summer when the thermistor depth was well known and the glacier surface was snow-free. A separate calculation of mean absolute error in hourly ice temperatures was made for each thermistor and time period, yielding 24 separate calibration values for r _eff (Fig. 3). The average optimal r _eff was ~0.1 mm, and the range was 0.01–0.2 mm. Because of uncertainty in the observations and calibration process, we used the ice temperature calibration to constrain the range of permissible r _eff values used in Section 4.2, rather than selecting a single value.

Fig. 3. Optimal effective ice grain radii, r _eff, for the Refreezing Model at TAR (black) and CAA (gray). Square markers indicate periods of time from early summer, circle markers midsummer, and triangle markers late summer. Markers connected by lines indicate thermistor strings with measurements at multiple depths during the same time period. Solid lines are from 2004/05 and dotted lines are from 2002/03. Optimal r _eff for the Drainage Model are similar.

Optimal r _eff is largely independent of z ₀ and χ because those parameters specify properties of the surface (Hoffman, 2011), and the Refreezing and Drainage models give similar optimal values of r _eff. Hoffman and others, (2008) calibrated r _eff using measurements of the downward solar flux within the ice and obtained a value of 0.24 mm. This is at the high end of the range of values found here through calibration to measured ice temperatures, but should be considered consistent, given the uncertainty of the two calibration methods. Hoffman and others, (2008) also noted bubbles of 0.05–0.5 mm radius in ice thin sections, which is very similar to the range of r _eff values found here. This is consistent with a modeling study that found that scattering within lake ice is controlled by the number and size of air bubbles trapped within the ice (Mullen and Warren, 1988).

4.2. Surface lowering

Doran and others, (2008) identified extreme high (2001/02) and low (2000/01) glacier runoff summers which exhibited a 10³-fold difference in annual streamflow throughout the MCM. Peak temperatures in December of each year averaged across 12 meteorological stations were 9.2°C (2001/02) and 4.4ºC (2000/01). To develop a model that has skill across the range of melt conditions observed in the Dry Valleys, we calibrated r _eff, z ₀ and for each model using the model root-mean-square error (RMSE) in summer surface lowering for these two extreme summers and the summer with the median observed surface lowering at TAR (2006/07). The range of eight z ₀ values sampled (Table 2) spans the range from smooth blue ice to melt-dominated surfaces (Brock and others, 2000). Sampled values of χ (Table 2) span 11–95% and are based on 14 values of the thickness of the surface layer. Additionally, r _eff is sampled using 15 values (Table 2) spanning the range constrained in Section 4.1 from ice temperature observations (Fig. 3). Modeled surface lowering in the Drainage Model is sensitive to r _eff because subsurface runoff lowers ice density, causing subsequent surface lowering to be more rapid (Fig. 2b). The Opaque Ice Model has only one free parameter, z ₀. After finding optimal parameter values for each model for the calibration period, we test the optimal models using the remaining eight to ten summers of surface lowering measurements available at each site.

Over the parameter space sampled, the Drainage Model can achieve substantially lower model error than the Refreezing Model (Fig. 4; Table 3). Additionally, as shown for TAR in Figure 4, the contours of zero error for the three calibration summers nearly overlap with the Drainage Model (for z ₀ = 0.01 mm), but with the Refreezing Model they remain far apart in parameter space, even for the lowest-error parameter combination. Similar conditions exist for the other two study sites. Finally, the Drainage Model performs better at validation of the remaining summers (Table 3). As expected, ablation from the Refreezing Model is largely independent of r _eff, while ablation from the Drainage Model is sensitive to r _eff. The Opaque Ice Model has about an order-of-magnitude greater error than the other models.

Fig. 4. RMSE in summer surface lowering at TAR for the three calibration summers using the Refreezing Model (a) and the Drainage Model (b). Vertical white lines indicate range of r _eff found from calibration using ice temperatures (Fig. 3). Gray contours are zero error for the three calibration years, and red stars show the optimal parameter combination for each model. The red square is the parameter combination in light blue in Figure 5. Not all values of z ₀ sampled are shown.

Table 3. Summary of optimal model configurations. For each column, the multiple values correspond to values for TAR /TAR2 /CAA

The variation in optimal parameters between models and study sites is likely to be due to a combination of physical differences and model uncertainty. Because the Drainage Model includes all the physical processes that are likely to be important, parameter differences between sites are more likely to be due to physical differences for this model. The optimal r _eff obtained through calibration of the Drainage Model to surface lowering at CAA is lower than that for the two sites on Taylor Glacier. This is consistent with r _eff differences between the glaciers obtained through calibration to subsurface ice temperatures (Fig. 3) and may be due to differences in ice history between the alpine Canada Glacier and the ice-sheet outlet Taylor Glacier. The optimal z ₀ with the Drainage Model at TAR2 is 1.5 orders of magnitude larger than at the other sites. This is consistent with the up to 20 m deep channels and basins found on the lower part of Taylor Glacier (Johnston and others, 2005), which would be expected to produce orders-of-magnitude increases in z ₀ (Smeets and Van den Broeke, 2008). Finally, the close agreement in χ between the three sites for the Drainage Model (20–24%) inspires confidence in the calibration, because there is not a physical reason for this parameter to vary between sites.

5. Results

5.1. Subsurface ice temperature

Modeled subsurface ice temperatures for both the Refreezing and Drainage models are within 1–2°C of observations for all summers with observations, which is expected because r _eff was calibrated using the observations. However, subsurface ice temperatures from the Opaque Ice Model are consistently ~3ºC lower than summer measurements (Fig. 5a), highlighting the importance of radiation penetration in warming the ice.

Fig. 5. (a) Measured (black) and modeled subsurface ice temperatures at TAR for 15 December 2005 to 20 February 2006. The configurations of the Refreezing (blue), Drainage (red) and Opaque Ice (gray) models are those defined in Table 3. Also shown in light blue is the Refreezing Model using a r _eff value that best reproduces ice temperatures for this summer (see Fig. 4a), which demonstrates how extra latent heat in the Refreezing Model can delay the fall freeze-up. (b) Contours of modeled melt extent. No subsurface melt was predicted by the Opaque Ice Model. The black line indicates the approximate depth of the thermistor shown in (a).

For most time periods, there is little difference in modeled ice temperature between the Refreezing and Drainage models when the same value of r _eff is used. For summer 2005/06 the calibrated value at TAR of r _eff = 0.12 mm for the Drainage Model allows the model to reproduce seasonal and diurnal temperatures very well, with error less than 0.5°C for most of the season (Fig. 5). The lower calibrated value of r _eff = 0.05 mm for the Refreezing Model produces temperatures ~1°C too low for most of this particular season. However, if the same r _eff value of 0.12 mm is used for the Refreezing Model, ice temperatures match observations very well until the end of summer (light blue line in Fig. 5a), with a minimal increase in error in surface lowering (Fig. 4). Despite the good match of the Refreezing Model to temperature observations in summer when using the larger grain radius, the model retains too much latent heat in the form of water in late summer when observations show the ice cooling rapidly; the Refreezing Model with r _eff = 0.12 mm requires an additional ~2 weeks longer than the Drainage Model to refreeze subsurface water (Fig. 5b).

These results on the sensitivity of retained latent heat are similar to a modeling study of Semtner, (1984) where the length of the melt season was found to be largely governed by latent heat stored in brine pockets in sea ice. Though both the Refreezing and Drainage models can predict ice temperatures during midsummer well, the Drainage Model is able to more accurately reproduce the timing of freeze-up. The Refreezing Model performs well in midsummer, but this is probably more related to the lack of large swings in air temperature during this period than to adequacy of process modeling.

5.2. Surface lowering and ablation

Using the optimal model configurations (Table 3), we calculated surface lowering for all summers at each site (Fig. 6). Both the Refreezing and Drainage models perform reasonably well at each site, but the Drainage Model is better able to reproduce the range of surface lowering observed at all three sites. The Opaque Ice Model substantially overestimates surface lowering in all summers. The Refreezing Model tends to overestimate surface lowering in summers when surface lowering is small, whereas the Drainage Model predicts more accurately (Fig. 6; Table 3). At TAR2 there is considerably more error for both the Refreezing and Drainage models, and both under-predict the high-ablation year of 2001/02 by close to 10 cm w.e.

Fig. 6. Measured surface lowering (black) and modeled surface lowering for the Opaque Ice (gray), Refreezing (blue) and Drainage (red) models for optimal model configurations at (a) TAR, (b) TAR2 and (c) CAA. The vertical lines indicate the component of surface lowering each summer from surface melt (solid lines) and subsurface drainage (dashed lines), with the remainder of surface lowering from sublimation. Summers are identified by the year in which they started, with bold text indicating the three summers used for calibration.

Both the Refreezing and Drainage models indicate that surface melt is a small fraction of summer surface lowering at TAR, with the majority of summer surface lowering due to sublimation (Fig. 6; Table 4). However, drainage of subsur-face melt becomes increasingly important as the amount of surface lowering increases when using the Drainage Model (Table 4). Furthermore, the importance of subsurface drainage increases if total ablation is considered, including mass loss remaining below the end-of-summer surface rather than just surface lowering; subsurface drainage accounts for half of all modeled mass loss at TAR for mean summer conditions and is the largest ablation component in the highest-ablation summer (2001/02) (Table 5). Thus, while we confirm the dominance of sublimation in surface lowering of MCM glaciers and the minor contribution from surface melt even in summer (Lewis and others, 1998; Hoffman and others, 2008), our model results indicate that drainage of subsurface melt may be the dominant component of surface mass balance in some locations.

Table 4. Summary of modeled surface lowering and ablation at TAR using the Drainage Model. Entries are given for the summers with minimum (2000/01) and maximum (2001/02) measured ablation and for the mean of all summers modeled. Surface lowering due to subsurface drainage is the near-surface mass loss from ice density reduction prior to subsequent surface ablation (sublimation, melt). Total ablation includes surface loss and all subsurface melt drainage including that beneath the end-of-summer surface. Mass loss values (cm w.e.) assume an ice density of 870 kg m^–3

Table 5. Summary of values used by selected glacier energy-balance studies. dz (1) is the thickness of the surface layer. G89=Greuell and Oerlemans, (1989); B95=Bintanja and Van den Broeke, (1995); vdB08=Van den Broeke and others, (2008). n/a=information not available

5.3. Ice density

The measured density of ice on glaciers in the Dry Valleys unaffected by ablation (early in the summer or at depths greater than 40 cm) was ~870 kg m^–3 (Fig. 7a). This is lower than the 900 kg m^–3 used by Fountain and others, (2006) and commonly assumed for glacier ice (e.g. Paterson, 1994; Braithwaite and others, 1998) but comparable to values reported for Antarctic blue ice (Bintanja, 1999). The two cores from early summer (Taylor: 1 December 2008; and Commonwealth: 18 November 2008) show no significant change in density with depth, while the other cores show lower values in the upper 10–20 cm (Fig. 7a). The core obtained from a rapidly melting block of calved ice at the base of Canada Glacier had a density of 740 kg m^–3 in the upper 7 cm. This is an upper bound as the shallowest few cm were crumbly and could not be sampled, which was typical of the upper few cm of cores obtained in mid- to late summer.

Fig. 7. (a) Ice density measurements. Vertical bars indicate the range of depths that each core sample represents, with the data point located at the average depth. Horizontal bars indicate uncertainty in density calculation. Core T1 was taken on Taylor Glacier, cores C1 and C2 were taken on Canada Glacier and cores M1–M3 were taken on Commonwealth Glacier. (b) Modeled ice density evolution at CAA in summer 2007/08. Colored lines are modeled density profiles at 10 day intervals. The thick blue line is measured ice density from core C1 (17 December 2007), with the vertical bars indicating the range of depths represented by each sample. The gray area corresponds to depths above which measurements were not possible due to the presence of crumbly, low-density ice, and the model results are represented with dashed lines in this region.

Because the Drainage Model allows ice mass loss beneath the surface, we are able to calculate the amount by which modeled subsurface drainage reduces ice density. In most summers the Drainage Model predicts one or more cycles of density lowering in the upper 10–15 cm of ice followed by excavation of the low-density layer by surface ablation (e.g. Fig. 7b). Modeled bulk ice densities from the Drainage Model compare favorably with the limited observations available. For example, on 17 December 2007, ice density at CAA was measured to range from 780 kg m^–3 near the surface to 870 kg m^–3 at depth (Fig. 7a and b). This corresponds to a mass loss of 0.8 cm w.e. below 3 cm depth, with an unknown mass loss above that depth because the ice density was too low to sample. The Drainage Model predicts a mass loss of 0.9 cm w.e. below 3 cm. However, the modeled density profile is biased with too low density at shallow depths and too high density at deeper depths. Despite the poor prediction of the density profile, the Drainage Model reproduces the depth-integrated subsurface mass loss well and is more realistic than the Refreezing Model, which does not allow density to lower.

5.4. Surface energy balance

The three study sites have similar surface energy budgets (Fig. 8). Net shortwave radiation makes the largest contribution of energy to the glaciers, but most of it penetrates into the ice, warming the subsurface (Fig. 8a). On average during summer, conduction brings ~80 W m^–2 of energy to the surface, making up for nearly all the shortwave radiation that penetrated the surface and absorbed in the interior. Net longwave radiation and latent heat in the form of sublimation are large energy losses from the surface during summer. Sensible heat is a small term in the summer energy balance that may be either positive or negative. This arises from both summer air temperature and ice surface temperature being only a few degrees below 0ºC throughout much of the summer, creating small temperature gradients between the surface and atmosphere. Melt is a small energy loss to the surface, even if considering only times when melt occurs.

Fig. 8. Surface energy-balance terms from the Drainage Model for (a) mean summer (December and January) conditions, (b) hours during which surface melt occurred and (c) hours during which subsurface melt occurred. For each time period, the three study sites are labeled. SW_pen is penetrated net solar radiation.

Surface melt occurs when net shortwave radiation is large and net longwave radiation is small (Fig. 8b). Additionally, latent heat loss tends to be small and sensible heat tends to be large. Considering the conditions required for subsurface melt, the only substantial difference from summer average conditions (cold ice surface) is higher net solar radiation (Fig. 8c), which highlights the role of subsurface solar heating in driving subsurface melt.

The surface-energy balances at TAR and CAA are similar, but TAR2 gains more energy from sensible heat, particularly during times when the surface is melting, reflecting air temperatures exceeding 0°C more frequently at this low-elevation site (Fig. 8). TAR2 has larger turbulent fluxes than the other two sites because of the larger surface roughness lengths found and used there (Table 3).

6. Discussion

6.1. Energy-balance and melt characteristics

The summer energy balance at Taylor Glacier calculated with our hourly time-step model is consistent with the energy balance calculated using a daily time-step (Hoffman and others, 2008). Both show that net shortwave radiation is the primary source of energy to the surface, and net longwave radiation and latent heat are the primary energy losses. The hourly model uses a much smaller solar radiation surface fraction (χ) value (and smaller surface thickness of <1 cm), as expected, so most of the net solar radiation penetrates the surface layer. At daily time-steps the larger χ, 82%, (and larger surface thickness, 13 cm) results in more solar radiation absorbed within the surface layer. Accordingly, conduction to the surface is an order of magnitude larger with the hourly model, because much of the solar radiation absorbed internally is returned to the surface.

The hourly (this study) and daily (Hoffman and others, 2008) models calculate very different magnitudes of surface and subsurface melt at TAR. The daily model, which uses a Refreezing assumption for subsurface melt, calculated 1.7 cm w.e. of surface melt and a seasonal maximum water column depth of 0.3 cm w.e., averaged over the 11 summers from 1995 to 2006. In contrast, the hourly time-step Refreezing Model calculates less surface melt (1.1 cm w.e.) and more subsurface melt (maximum water column depth 2.2 cm), averaged over 13 summers. The hourly Drainage Model calculates even less surface melt (0.2 cm w.e.) but more subsurface melt (a maximum water column depth of 2.0 cm and 2.6 cm of drained subsurface melt contributing to ablation during the average summer).

Despite the difference in surface layer thickness between the daily (13 cm) and hourly (1 cm) Refreezing Models, both models predict similar fractions of subsurface melt in the upper 13 cm (92% of all melt) and exhibit similar errors: 2.3 cm w.e. in summer ablation for the daily model (Hoffman and others, 2008) and 1.9 cm w.e. calibration error in summer surface lowering for the hourly Refreezing Model (Table 3). In contrast, the hourly Drainage Model (surface layer thickness 0.25 cm) predicts a different subsurface melt profile, with only 51% of subsurface melt in the upper 13 cm, and has lower error (0.2 cm w.e. calibration error). Based on these comparisons, we find that including the drainage process in the model is more important for improving model quality than increasing the temporal resolution from daily to hourly time-steps alone. We also infer that if one wishes to adjust the vertical resolution of melt prediction near the surface, the time-step must be adjusted accordingly, and it is not possible to resolve fine vertical structure with a coarse time-step.

Overall, as discussed in Hoffman and others, (2008), the summer energy balance calculated here is similar to results from elsewhere in Taylor Valley (Lewis and others, 1998; Johnston and others, 2005), to Antarctic blue ice areas (Bintanja and Van den Broeke, 1995) and to glaciers on Mount Kilimanjaro, Tanzania, at ~5800 m a.s.l. (Mölg and Hardy, 2004). The difference between energy-balance conditions at these polar and high-elevation locations and conditions at temperate glaciers is due to the magnitude of sensible heat and, to a lesser extent, latent heat.

6.2. Drainage of subsurface melt

Our study demonstrates the importance of solar radiation penetration for the surface energy balance of glacier ice and the importance of subsequent subsurface melt and drainage for the surface mass balance of glaciers with low ablation rates. Assuming that all net solar radiation is absorbed at the ice surface (Opaque Ice Model) results in modeled surface lowering that is 6–15 cm too large (Table 3; Fig. 6) and modeled subsurface ice temperatures that are 3°C too low (Fig. 5). Inclusion of transmission of solar radiation into the ice allows the model to accurately reproduce both surface lowering and ice temperatures.

Including the drainage of subsurface melt in the model improves modeled surface lowering, subsurface ice temperatures and subsurface ice density. The Drainage Model exhibits smaller error in surface lowering than the Refreezing Model (Table 3), better matches observations in extreme seasons (Fig. 6) and better matches subsurface ice temperatures, particularly in late summer. Ignoring the drainage of subsurface melt can result in overestimating late-summer temperatures by 3°C (Fig. 5). Additionally, the Refreezing Model predicts unrealistically high subsurface water fractions of 100% in extreme years (e.g. 2001/02). These model comparisons show strong evidence for a drain-cooling of the ice that is analogous to the refreeze-warming of melting snowpacks (Pfeffer and Humphrey, 1998). Lastly, only the Drainage Model can reproduce the reductions in ice density over summer that are observed (Fig. 7a), and it is able to estimate the approximate magnitude of this mass loss (Fig. 7b). We do note, however, that both models perform poorly for the extreme melt summer of 2001/02 at TAR2 (Fig. 6b), suggesting one or more of: problems with that observation, inadequacy of both models in the transition to a high-melt environment, or a change in model parameters not accounted for in our simulations (e.g. a change in ice optical properties (Buckley and Trodahl, 1987; Mullen and Warren, 1988; Gardner and Sharp, 2010) or surface roughness (Smeets and Van den Broeke, 2008) associated with the high ablation rates during this summer).

Based on the differences between the Refreezing and Drainage models, the modeling indicates an annual cycle of ice density variations in the near-surface layer (upper ~20 cm) of the ice on Dry Valley glaciers (Fig. 9). Ice far below the surface has a density of 870 kg m^–3, and, at the beginning of summer, ice at the surface has this density. During summer, drainage of subsurface melt gradually lowers the density in the near-surface layer. This mass loss is centered on a depth 5–10 cm below the surface where subsurface melting is greatest. As the surface ablates, internal melt occurs at greater depths and the density profile is maintained. At summer’s end, melt ceases, leaving a surficial layer of reduced density. During winter, sublimation ablates this layer away, revealing denser ice at the surface at the start of the following summer. This process is comparable to the cycles of weathering crust development and elimination observed over the course of days by Müller and Keeler, (1969) on White Glacier, but it occurs over an annual cycle in Taylor Valley due to lower ablation rates.

Fig. 9. Schematic of annual cycle of ablation and ice density on MCM glacier ablation zones. The orange, italicized numbers are representative ice densities (kg m^–3). Figure is not to scale; the magnitude of summer and winter surface lowering is typically 10–20 cm each.

Our modeling suggests that a substantial fraction of ablation, and much of the meltwater, from MCM glaciers is subsurface. Therefore surface mass-balance measurements using traditional stake methods and other surface lowering methods (e.g. sonic rangers) will underestimate mass loss during summer and overestimate it during winter, with seasonal errors in mass balance on the order of 10–50% depending on ablation rate. Of course, these errors will not affect annual mass-balance measurements or calculations of cumulative mass balance (e.g. Fountain and others, 2006). Ice weathering is a process common to all ice surfaces subject to melting conditions in the presence of solar radiation (Müller and Keeler, 1969; Larson, 1977; Braithwaite and Olesen, 1990; Hubbard and Glasser, 2005). However, in climates yielding higher melt rates, the importance of this process to surface mass balance is reduced (Müller and Keeler, 1969; Van den Broeke and others, 2008).

6.3. Partitioning of net shortwave radiation

Little agreement exists on the proper partitioning of net shortwave radiation between the surface and the subsurface. Many glacier energy-balance studies assume that all net shortwave radiation is absorbed at the surface (i.e. χ = 100%). For temperate glaciers with high melt rates this may be a reasonable simplification. Rapid ablation quickly exhumes ice with density reduced due to internal melting, so there is little net difference whether radiation penetration is included or not (e.g. Van den Broeke and others, 2008).

When the near-surface ice is below the melting point, many energy-balance studies divide the net shortwave radiation into a surface fraction and a penetrated fraction, as done here. However, the partitioning used varies significantly, as do the surface layer thickness and time-steps in the associated models. Many published studies use solar penetration methods based on Greuell and Oerlemans, (1989) who divided the net shortwave radiation into two wavelength bands. The infrared portion (wavelengths >800 nm) is assumed to be completely absorbed in a constant 6 cm thick surface layer while the visible portion completely penetrates the surface layer, yielding χ = 36%. The penetrating portion is assumed to be absorbed exponentially with depth according to Bier’s law using a broadband extinction coefficient. This χ value of 36% has been adopted by other studies using a range of surface layer thicknesses from 3 to 12 cm and time-steps ranging from 1 min to 1 hour, while other studies have used the same approach but with a χ value of 80% (Table 5). More recently, Van den Broeke and others, (2008) and Kuipers Munneke and others, (2009) implemented explicit modeling of solar radiation penetration using spectrally dependent extinction coefficients (Brandt and Warren, 1993) for ice and snow in Greenland (Table 5), similar to the approach presented here.

The use of χ can be compared to the practice commonly used in sea-ice models of specifying the flux of radiative energy that passes through the surface into the ice, commonly represented as I ₀, and removing that flux from the surface energy balance (Maykut and Untersteiner, 1971; Launiainen and Cheng, 1998). In sea-ice modeling, the surface layer is often defined to be 10 cm thick, as this is the depth by which most near-infrared energy has been absorbed, leaving only visible light, which has a spectrally more uniform extinction coefficient, below this depth. This approach is similar to the method used by Greuell and Oerlemans, (1989).

Little discussion exists in the glacier or sea-ice literature regarding the uncertainty or variability of χ or I ₀, despite the conclusion by the authors of the original I ₀ model formulation that more careful studies of I ₀ were needed based on a sensitivity analysis (Maykut and Untersteiner, 1971). Maykut and Untersteiner, (1971) and Semtner, (1984) found that modeled sea-ice thickness increases as more solar radiation is partitioned into the subsurface because less energy is available at the surface for melting. We also find surface melt to be highly sensitive to the partitioning of solar radiation (Fig. 5). The length of the sea-ice melt season is largely governed by latent heat stored in brine pockets in the ice (Semtner, 1984), an effect we see for some summers (Fig. 5).

However, most energy-balance studies perform well despite the simplification of complete radiation absorption at the surface, and this suggests that partitioning of short-wave radiation penetration is not crucial for many applications. Most of the studies reviewed consider conditions where melt is substantial (e.g. Van den Broeke and others, 2008) or absent (e.g. Bintanja and Van den Broeke, 1995). However, we show that for ice conditions near the melting point small differences in energy receipt and distribution can result in large changes in melt. Consequently, for glaciers in the MCM, results are very sensitive to χ. For the best-fit hourly Refreezing Model, χ at TAR is 40%, half the value needed for the daily time-step (82%; Hoffman and others, 2008). As discussed previously, the melt estimates of those two models are broadly consistent, suggesting that the appropriate χ value is dependent on the chosen time-step. For the Drainage Model with an hourly time-step, the calibrated values of χ for all three sites studied here have a narrow range of 21–24% (Table 3).

7. Conclusion

In a previous study (Hoffman and others, 2008), we demonstrated the importance of including the partitioning of solar radiation into surface and subsurface pathways for the climate conditions found in the MCM. Here we show that subsurface melt, caused by penetration of solar radiation, is substantial and its drainage is a significant factor in the ablation of these glaciers. By including subsurface drainage, the model better reproduces observations of surface lowering, subsurface ice temperature, and ice density. Compared with the previous study, we find that increasing the temporal resolution of the model from daily to hourly time-steps provides minimal improvement in model skill at the seasonal timescales, but including the drainage process improves overall model skill considerably. These results indicate this approach is critical to predicting meltwater runoff from the glaciers in the Dry Valleys and similar environments.

Calculations show that drainage of subsurface melt can contribute 10–50% of observed surface lowering, with additional mass loss beneath the surface each summer. It can be the largest component of total ablation in summer if mass loss from weathered ice beneath the end-of-summer surface is included, and is typically an order of magnitude larger than melt generated on the surface. Subsurface melt and drainage play an important role in the annual cycle of surface mass balance by lowering ice density in the upper 20 cm. The low-density ice is ablated by sublimation each winter, leaving unweathered ice at the surface each spring. This annual process of weathering crust development and elimination on MCM glaciers is observed on any glaciers experiencing melt, but in higher-melt environments the cycle takes days to weeks instead of a full year (Müller and Keeler, 1969; Van den Broeke and others, 2008). Measurements of surface mass balance using the glaciological method in low-melt environments need to account for cycles of ice weathering to avoid biasing seasonal observations.

The low-melt environment of the MCM emphasizes the importance of properly partitioning net solar radiation between the surface and the subsurface when modeling the energy balance. Variations in partitioning of 10% can result in modeled differences in ablation of 5 cm w.e. or more in some cases, which for MCM glaciers can be 25– 100% of summer ablation. We found the partitioning appears to be dependent on model time-step. The hourly and daily versions of the model illustrated the importance of model time-step to the physical understanding gained. While both highlighted the importance of solar radiation as an energy source, and longwave radiation and latent heat as the important heat losses, the lessons of subsurface physical processes and controls differed in important ways. Numerical constraints dictated a thicker ‘surface’ for the daily time-step compared to the hourly time-step. The shorter time-step allowed a finer discrimination of internal melt in the near surface. Model results showed more internal melt for the hourly timescale and less surface melt, which is broadly consistent with the anectodal experience of rare observations of surface melt on broad, flat glacier surfaces in the MCM. Our review of published reports reveals that few studies make a detailed assessment of this parameter; however, higher-melt environments are probably less sensitive to its magnitude.