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 (Reference Liston, Winther, Bruland, Elvehøy and SandListon and others, 1999; Reference Hoffman, Fountain and ListonHoffman 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% (Reference Coléou and LesaffreColéou and Lesaffre, 1998) and observations of seaice permeability reducing to 0 when porosity is less than 5–10% due to capillary tension (Reference Golden, Eicken, Heaton, Miner, Pringle and ZhuGolden 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 (Reference Fountain, Nylen, Monaghan, Basagic and BromwichFountain 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 (Reference HoffmanHoffman, 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 (Reference Fountain, Nylen, Monaghan, Basagic and BromwichFountain 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 (Reference ChinnChinn, 1981; Reference FountainFountain and others, 1999a; Reference GooseffGooseff 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 (Reference HoffmanHoffman, 2011). The influence of this simple snow accounting (or not accounting) on the simulations and results presented herein is thought to be negligible.

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)) (Reference Brandt and WarrenBrandt and Warren, 1993; Reference Liston, Winther, Bruland, Elvehøy and SandListon 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 (Reference Brandt and WarrenBrandt and Warren, 1993; Reference Liston, Winther, Bruland, Elvehøy and SandListon and others, 1999; Reference Hoffman, Fountain and ListonHoffman 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 (Reference Mullen and WarrenMullen 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 (Reference HoffmanHoffman, 2011), and the Refreezing and Drainage models give similar optimal values of r _eff. Reference Hoffman, Fountain and ListonHoffman 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. Reference Hoffman, Fountain and ListonHoffman 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 (Reference Mullen and WarrenMullen and Warren, 1988).

4.2. Surface lowering

Reference DoranDoran 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 (Reference Brock, Willis, Sharp and ArnoldBrock 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 (Reference Johnston, Fountain and NylenJohnston and others, 2005), which would be expected to produce orders-of-magnitude increases in z ₀ (Reference Smeets and Van den BroekeSmeets 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.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 Reference SemtnerSemtner, (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 (Reference Lewis, Fountain and DanaLewis and others, 1998; Reference Hoffman, Fountain and ListonHoffman 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=Reference Greuell, Oerlemans and OerlemansGreuell and Oerlemans, (1989); B95=Reference Bintanja and Van den BroekeBintanja and Van den Broeke, (1995); vdB08=Reference Van den Broeke, Smeets, Ettema, Van der Veen, Van de Wal and OerlemansVan 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 Reference Fountain, Nylen, MacClune and DanaFountain and others, (2006) and commonly assumed for glacier ice (e.g. Reference PatersonPaterson, 1994; Reference Braithwaite, Konzelmann, Marty and OlesenBraithwaite and others, 1998) but comparable to values reported for Antarctic blue ice (Reference BintanjaBintanja, 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.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 (Reference Hoffman, Fountain and ListonHoffman 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 (Reference Hoffman, Fountain and ListonHoffman 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 (Reference Hoffman, Fountain and ListonHoffman 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 Reference Hoffman, Fountain and ListonHoffman and others, (2008), the summer energy balance calculated here is similar to results from elsewhere in Taylor Valley (Reference Lewis, Fountain and DanaLewis and others, 1998; Reference Johnston, Fountain and NylenJohnston and others, 2005), to Antarctic blue ice areas (Reference Bintanja and Van den BroekeBintanja and Van den Broeke, 1995) and to glaciers on Mount Kilimanjaro, Tanzania, at ~5800 m a.s.l. (Reference Mölg and HardyMö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 (Reference Pfeffer and HumphreyPfeffer 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 (Reference Buckley and TrodahlBuckley and Trodahl, 1987; Reference Mullen and WarrenMullen and Warren, 1988; Reference Gardner and SharpGardner and Sharp, 2010) or surface roughness (Reference Smeets and Van den BroekeSmeets 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 Reference Müller and KeelerMü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. Reference Fountain, Nylen, MacClune and DanaFountain and others, 2006). Ice weathering is a process common to all ice surfaces subject to melting conditions in the presence of solar radiation (Reference Müller and KeelerMüller and Keeler, 1969; Reference LarsonLarson, 1977; Reference Braithwaite and OlesenBraithwaite and Olesen, 1990; Reference Hubbard and GlasserHubbard and Glasser, 2005). However, in climates yielding higher melt rates, the importance of this process to surface mass balance is reduced (Reference Müller and KeelerMüller and Keeler, 1969; Reference Van den Broeke, Smeets, Ettema, Van der Veen, Van de Wal and OerlemansVan 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. Reference Van den Broeke, Smeets, Ettema, Van der Veen, Van de Wal and OerlemansVan 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 Reference Greuell, Oerlemans and OerlemansGreuell 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, Reference Van den Broeke, Smeets, Ettema, Van der Veen, Van de Wal and OerlemansVan den Broeke and others, (2008) and Kuipers Reference MunnekeMunneke and others, (2009) implemented explicit modeling of solar radiation penetration using spectrally dependent extinction coefficients (Reference Brandt and WarrenBrandt 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 (Reference Maykut and UntersteinerMaykut and Untersteiner, 1971; Reference Launiainen and ChengLauniainen 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 Reference Greuell, Oerlemans and OerlemansGreuell 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 (Reference Maykut and UntersteinerMaykut and Untersteiner, 1971). Reference Maykut and UntersteinerMaykut and Untersteiner, (1971) and Reference SemtnerSemtner, (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 (Reference SemtnerSemtner, 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. Reference Van den Broeke, Smeets, Ettema, Van der Veen, Van de Wal and OerlemansVan den Broeke and others, 2008) or absent (e.g. Reference Bintanja and Van den BroekeBintanja 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%; Reference Hoffman, Fountain and ListonHoffman 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).