Representing moisture fluxes and phase changes in glacier debris cover using a reservoir approach

Due to the complexity of treating moisture in supraglacial debris, surface energy balance models to date have neglected moisture infiltration and phase changes in the debris layer. The latent heat flux (QL) is also often excluded due to the uncertainty in determining the surface vapour pressure. To quantify the importance of moisture on the surface energy and climatic mass balance (CMB) of debris-covered glaciers, we developed a simple reservoir parameterization for the debris ice and water content, as well as an estimation of the latent heat flux. The parameterization was incorporated into a CMB model adapted for debris-covered glaciers. We present the results of two point simulations, using both our new "moist" and the conventional "dry" approaches, on the Miage Glacier, Italy, during summer 2008 and fall 2011. The former year coincides with available in situ glaciological and meteorological measurements, including the first eddy-covariance measurements of the turbulent fluxes over supraglacial debris, while the latter contains two refreeze events that permit evaluation of the influence of phase changes. The simulations demonstrate a clear influence of moisture on the glacier energy and mass-balance dynamics. When water and ice are considered, heat transmission to the underlying glacier ice is lower, as the effective thermal diffusivity of the saturated debris layers is reduced by increases in both the density and the specific heat capacity of the layers. In combination with surface heat extraction by QL, subdebris ice melt is reduced by 3.1% in 2008 and by 7.0% in 2011 when moisture effects are included. However, the influence of the parameterization on the total accumulated mass balance varies seasonally. In summer 2008, mass loss due to surface vapour fluxes more than compensates for the reduction in ice melt, such that the total ablation increases by 4.0%. Conversely, in fall 2011, the modulation of basal debris temperature by debris ice results in a decrease in total ablation of 2.1%. Although the parameterization is a simplified representation of the moist physics of glacier debris, it is a novel attempt at including moisture in a numerical model of debris-covered glaciers and one that opens up additional avenues for future research.


Introduction
Numerical modelling of debris-covered glaciers has received renewed scientific interest in recent years, because their contribution to changes in ice mass and water resources in many regions remains poorly understood (e.g. Kääb et al., 2012) and because the proportion of debris-covered glacier area is rising as glaciers recede (e.g. Stokes et al., 2007;Bolch et al., 2008;Bhambri et al., 2011).
It is well established that supraglacial debris exerts an important control on glacier melt rates. Subdebris ice melt is strongly enhanced when the debris thickness is less than a few centimetres, due to a reduction in surface albedo, an increase in absorption of shortwave radiation, and the rapid transfer of energy to the underlying ice. Melt decreases exponentially as the thickness increases, as a result of insulation of the underlying glacier ice from the overlying atmosphere (e.g. Østrem, 1959;Loomis, 1970;Fujii, 1977;Inoue and Yoshida, 1980;Mattson et al., 1993). The presence of debris also alters the glacier surface energy balance, by permitting surface temperatures to rise above the melting point and by altering surface heat and moisture exchanges with the atmosphere (e.g. Brock et al., 2010).
Numerous point models of the surface energy balance of debris-covered glaciers have been developed to simulate subdebris ice melt (e.g. Kraus, 1975;Nakawo and Young, 1982;Han et al., 2006;Nicholson and Benn, 2006;Reid and Brock, 2010). In recent years, models of debris cover have been extended to distributed simulations (Zhang et al., 2011), and to include both explicit calculation of heat conduction through debris layers resolved into multiple levels and snow accumulation on top of the debris (Reid et al., 2012;Lejeune et al., 2013;Fyffe et al., 2014).
However, due to the complexity of treating moisture in supraglacial debris cover, surface energy balance models to date have neglected the latent heat and surface moisture flux components, with the exception of (1) testing the two endmember cases of completely dry or completely saturated debris layers (e.g. Nakawo and Young, 1981;Kayastha et al., 2000;Nicholson and Benn, 2006), and (2) using measurements of surface relative humidity to calculate the flux when the surface is saturated (Reid and Brock, 2010;Reid et al., 2012). In addition, moisture inputs to the debris layer -by percolation of snowmelt and rainfall, or from the underlying melting ice via capillary action -and their phase changes have not been taken into account. Rather, any water is assumed to run off immediately, without influencing the thermal properties of the debris (e.g. Reid and Brock, 2010;Reid et al., 2012;Lejeune et al., 2013;Fyffe et al., 2014).
Both field observations and laboratory experiments indicate that debris covers can be partially or entirely saturated at times during the ablation season, depending on their thickness and the environmental conditions, with a minimum of a saturated region adjacent to the interface if the underlying ice is at the melting point (e.g. Nakawo and Young, 1981;Conway and Rasmussen, 2000;Kayastha et al., 2000;Reznichenko et al., 2010;Nicholson and Benn, 2012). The presence of interstitial water and ice modifies the thermal properties of the debris layer, particularly during transition seasons (e.g. Conway and Rasmussen, 2000;Nicholson and Benn, 2012). In addition, percolation of rain through a debris layer, which can reach as high as 75 % of the total rainfall at the surface (Sakai et al., 2004), and other inputs of moisture can influence the thermal regime by heat advection (Reznichenko et al., 2010), and by providing a source of moisture for evaporation that cools the debris and therefore reduces heat transmission to the ice below.
Surface vapour exchanges between the debris and the overlying atmosphere influence the surface energy balance and have been observed to be non-negligible at times. Sakai et al. (2004) estimated that the ablation calculated by an energy balance approach that neglects the latent heat flux, QL, would provide an overestimate of up to 100 %, since its lowering effect on surface temperature would not be captured. During the ablation season on the Miage Glacier in the Italian Alps, Brock et al. (2010) calculated large spikes in QL, of up to −800 W m −2 , that coincided with daytime rainfall events on the heated debris surface. Furthermore, while they estimated that energy inputs due to condensation and deposition were negligible, ground frosts were observed on a weekly to biweekly basis in the upper parts of the glacier, which may have slowed early daytime heating of the debris layer. Given the clear influence of moisture on the surface energy balance and the subsurface thermal regime, there is a need to develop a treatment for moisture fluxes into and within the debris layer, as well as for phase changes, that would allow for a variation in the thermal properties and energy sources and sinks of the debris layer with depth and time.
In this paper, we explore the utility of a reservoir scheme for parameterizing moisture fluxes and phase changes in a glacier debris layer that has been incorporated into a glacier climatic mass balance model. We exploit a short period of available in situ measurements over supraglacial debris to evaluate the model performance during an ablation season, with a second simulation of a fall season to fully demonstrate the capabilities of the model. Within the context of the simplified parameterization, we show the influence of moisture on heat transfer in the debris layer, its physical properties, and subdebris ice melt, as well as assess the scale of the impact of phase changes. The eventual goal of this work is to incorporate the debris modifications into an interactively coupled modelling system of the atmosphere and alpine glaciers at the regional scale (Collier et al., 2013). The inclusion of debris is essential for (1) accurately capturing surface conditions over debris-covered glaciers and, therefore, atmosphere-glacier feedbacks, and (2) rigorously assessing regional climatic influences on the CMB of debriscovered glaciers.

Debris-free glacier CMB model
The debris-free version of the glacier CMB model is described in detail by Mölg et al. (2008Mölg et al. ( , 2009Mölg et al. ( , 2012. The model has been applied to simulating glaciers in a wide variety of climatic settings (e.g. Mölg et al., 2012;Collier et al., 2013;Nicholson et al., 2013;MacDonell et al., 2013). The CMB model solves the surface energy balance equation to determine the energy available for melt and other mass fluxes, given by where the terms correspond to, from left to right, net shortand longwave radiation, turbulent fluxes of sensible and  Brock et al. (2010) latent heat, the ground heat flux (composed of conduction and penetrating shortwave radiation) and the heat flux from precipitation. Following the convention in mass balance modelling, fluxes are defined as positive when energy transfer is to the surface. The residual energy flux, F NET , constitutes the energy available for melt provided the surface temperature has reached the melting point. The specific mass balance is calculated from solid precipitation, surface vapour fluxes, surface and subsurface melt, and refreeze of liquid water in the snowpack. Surface vapour fluxes (M v ; i.e. sublimation or deposition (kg m −2 ) depending on the sign of QL) at each time step t are calculated according to where L H is the latent heat of sublimation (2.84×10 6 J kg −1 ) or vaporization (2.51 × 10 6 J kg −1 ), depending on the surface temperature. The CMB model treats numerous additional processes, including the evolution of surface albedo and roughness based on snow depth and age; snowpack compaction and densification by refreeze; and the influence of penetrating solar radiation, refreeze and conduction on the near-surface englacial temperature distribution. Physical parameter values for snow and ice are provided in Table 1.

Inclusion of debris
For this study, the glacier CMB model was modified to include a treatment for supraglacial debris according to two cases: (1) one with no treatment of moisture fluxes or phase changes in the debris layer, congruent with previous studies (CMB-DRY); and (2) one that introduces a reservoir to parameterize the moisture content of the debris layer and its phase, and also includes a latent heat flux calculation (CMB-RES). The simulations are performed as point simulations, due to the availability of both meteorological-forcing and evaluation data at a single location.

Surface temperature
Consistent with previous modelling studies of debris-covered glaciers (e.g. Nicholson and Benn, 2006;Reid and Brock, 2010;Reid et al., 2012;Zhang et al., 2011), the model employs an iterative approach to prognosing surface temperature, with the solution yielding zero residual flux in the surface energy balance (Eq. 1). The model uses the Newton-Raphson method to calculate T SFC at each time step, as implemented in Reid and Brock (2010), with a different termination criteria of |F NET | < 1 × 10 −3 . When snow or ice are exposed at the surface, the resulting T SFC is reset to the melting point if it exceeds this value, and energy balance closure is achieved by using the residual energy for surface melt.

Subsurface temperature
Both versions of the CMB model prognose the temperature distribution in the upper subsurface following the conservation of energy. The vertical levels selected for the case study in Sect. 2.3 are defined in Table 2, and are set at fixed depths in the subsurface, from 0.0 to 9.0 m, that track the glacier surface as it moves due to mass loss or gain. On this grid, the 1-D heat equation becomes where ρ is the density (kg m −3 ); c is the specific heat capacity (J kg −1 K −1 ); T is the englacial temperature (K); k is the thermal conductivity (W m −1 K −1 ); and Q is the heat flux due to non-conductive processes (penetrating shortwave radiation; W m −2 ). For these simulations, the numerical scheme used to solve Eq. (3) was updated from a centred-difference approach to a Crank-Nicolson scheme, which was solved following Smith (1985). The greater stability of the numerics permits the subsurface layer spacing throughout the debris to be decreased to 1 cm from ∼ 10 cm previously. The vertical grid spacing is thus consistent with the small number of previous studies that explicitly simulate heat conduction in the debris (Reid and Brock, 2010;Reid et al., 2012;Lejeune et al., 2013), rather than assume that the temperature gradient is  With the exception of Lejeune et al. (2013), the ice temperature in previous modelling studies has been assumed to be at the melting point, due to the focus on the ablation season (e.g. Nicholson and Benn, 2006;Reid and Brock, 2010). Although this assumption has been validated by field measurements (e.g. Conway and Rasmussen, 2000;Brock et al., 2010), it limits the temporal applicability of the model and may contribute to the overestimation of night-time surface temperatures when the overlying air temperature drops below the melting point (Reid and Brock, 2010). The CMB models explicitly simulate heat conduction throughout the glacier column. Therefore, the ice temperature is a prognostic variable at all levels except the bottom boundary, where a zero-flux condition is imposed. Finally, subsurface heating due to penetrating shortwave radiation is not considered when glacier debris is exposed at the surface (e.g. Reid and Brock, 2010).

Physical and thermal properties
The important physical properties of the glacier subsurface in Eq. (3) -density ρ, thermal conductivity k, and specific heat capacity c -are non-uniform with depth. Defining m S and m D as the levels corresponding to the bottom of the snowpack and debris layers (cf. Fig. 1), respectively, the column properties (generalized as f (z)) are specified as Standard values are selected for snow and glacial ice properties (Table 1), with the exception of snow density, which is a prognostic variable. Within the debris layer, the properties of each 1 cm layer are a weighted average of the depth-invariant whole-rock values, f wr , and the content of the pore space, f φ , as determined by an assumed linear porosity function, φ: For CMB-DRY, the debris pore space contains only air (f φ = f air ), while the weighted average in CMB-RES also considers the bulk water and ice content of the debris of saturated layers. The porosity function is discussed further in Sect. 2.3.

Moisture in the debris layer
For CMB-DRY, rainfall or other liquid water inputs are instantaneously removed as runoff from the debris layer and do not accumulate or contribute to vapour exchanges between the debris and the atmosphere, similar to previous modelling studies (e.g. Reid and Brock, 2010;Reid et al., 2012;Lejeune et al., 2013). For CMB-RES, a reservoir is introduced for moisture accumulation and phase changes ( Fig. 1). The reservoir depth for each column is calculated as the sum of the debris porosity over the debris thickness. Thus, the pore space in the debris is represented as a single reservoir, rather than the storage in each 1 cm layer being treated individually. Liquid water, from rainfall or melt of the overlying snowpack, instantly infiltrates the reservoir. The location of the water and/or ice in the debris is not prognosed; rather, moisture is assumed to occupy the lowest debris layers, adjacent to the glacier ice.
In addition, when the ice-debris interface reaches the melting point, a minimum debris-water content is imposed to reflect field observations of a basal saturated layer during the ablation season (e.g. Nakawo and Young, 1981;Conway and Rasmussen, 2000;Kayastha et al., 2000;Reznichenko et al., 2010;Nicholson and Benn, 2012). As the water content of glacier debris cover is poorly constrained and no measurements are available, the minimum value is set to the amount of water needed to saturate the lowest 1 cm layer in the debris, given its porosity and ice content. The horizontal drainage of debris water is accounted for using a simplistic representation of the runoff timescale, which is a linear function of terrain slope and varies from 1 to 0 h −1 between 0 • and 90 • (Reijmer and Hock, 2008).
Congruent with the simple nature of the reservoir parameterization, the heat flux from precipitation is only applied at the surface in CMB-RES, and subsurface heat transport by water percolation is not included. This treatment is consistent with the findings of Sakai et al. (2004), namely that the heat flux due to rainfall percolation contributes minimally to subdebris ice melt, although its influence may depend on debris permeability (Reznichenko et al., 2010).

Turbulent fluxes of latent and sensible heat
The turbulent fluxes of sensible heat (both models) and latent heat (CMB-RES) were computed using bulk aerodynamic formulae and corrected for atmospheric stability according to the bulk Richardson number, as is standard in glacier energy balance modelling (e.g. Braithwaite, 1995;Reid and Brock, 2010). The bulk Richardson number was constrained within reasonable limits following Fyffe et al. (2014), with the correction applied to fewer than 2.0 % of total time steps in the simulations. The latent heat flux in CMB-DRY was set to zero, to be consistent with previous studies of debris-covered glaciers, as no measurements of surface relative humidity were available. For CMB-RES, the surface vapour pressure was needed but unknown.
For the case study described in Sect. 2.3, an automatic weather station (AWS) measured relative humidity at a height of z air = 2.16 m, from which the partial vapour pressure was calculated. The partial density of water vapour was then obtained from where the symbols correspond to, from left to right, the air's water vapour partial pressure, the partial density of water vapour, the specific gas constant for water vapour (461.5 J kg −1 K −1 ), and the air temperature at a height of z air . In this study, we assumed that ρ vap air is constant between the sensor and the surface of the debris layer, i.e. that water vapour in the atmospheric surface layer is well mixed. The vapour pressure at the surface is therefore given by For a completely unsaturated glacier debris layer in CMB-RES, a latent heat flux would nonetheless arise due to the vapour pressure gradient that results from the temperature difference between the surface and z air . However, when water or ice are present in the debris, the final calculation of the surface vapour pressure, e sfc , includes a linear correction towards the saturation value e sfc sat at T sfc according to where e * sfc is the initial guess in Eq. (7); air is the void fraction of the bulk layer that is occupied by air; and φ bulk is the bulk debris porosity, which is invariant under different debris thicknesses due to the linear specification of φ (as described in Sect. 2.3). air is given by where m K is the level of the saturated horizon in the debris and N is the total number of layers in the debris. When the debris is completely unsaturated, air = φ bulk , and when it is completely saturated, air = 0. Therefore, the surface vapour pressure in CMB-RES is a linear function of the moisture content of the reservoir rather than a wetted debris surface: as the reservoir fills from infiltration of rainfall or snowmelt, the distance between the surface and the saturated horizon (represented by air ) decreases and e sfc approaches saturation.

Mass balance
The total mass balance calculation in CMB-DRY and CMB-RES accounts for the following mass fluxes (kg m −2 ) at each time step: solid precipitation, surface and vertically integrated subsurface melt, meltwater refreeze and formation of superimposed ice in the snowpack, changes in liquid water storage in the snowpack, and surface vapour fluxes. The contribution of surface vapour fluxes to or from the debris layer is zero when overlying snow cover is present and in CMB-DRY. In CMB-RES, these fluxes also contribute to changes in the debris water and ice content of the reservoir. For both models, subdebris ice melt is calculated as the vertical integral of melt in the ice column underlying the debris.
Liquid precipitation contributes indirectly to the mass balance in both CMB models through changes in storage in the snowpack, and contributes directly in CMB-RES, via reservoir storage. However, changes in the debris water and ice content in CMB-RES are not included in the mass balance calculation, so as to allow for a more direct comparison between CMB-RES and CMB-DRY of the influence of including the latent heat flux. The impact of changes in the storage of water and ice in the debris is quantified in Sect. 3 and has a negligible influence on the total accumulated mass balance.

Miage Glacier case study
The study area is the Miage Glacier in the Italian Alps (45 • 47 N, 6 • 52 E; Fig. 2). This glacier was selected due to the availability of meteorological data from an automatic weather station (AWS) located on the lower, debris-covered part of the glacier at an elevation of 2030 m a.s.l. The debris thickness was determined by a point measurement to be 23 cm. At the surface, the debris is composed mainly of coarse gravel and cobbles, ranging in size from a few centimetres to 25 cm, with occasional larger rocks, 1-2 m in size.   Figure 2. Map showing the location of the Miage Glacier. The AWS located on the glacier is denoted with a red circle and the AWS2 from which precipitation data were obtained is shown by a red triangle.
The AWS site was deliberately chosen to be upwind from any nearby large boulders. We performed two simulations, one for summer 2008 and one for fall 2011. The former covered the period of 25 June-11 August 2008, with the first 25 days discarded as model spin-up time. For much of the 2008 simulation, the AWS provided hourly values of air temperature, vapour pressure, wind speed, and incoming short-and longwave radiation (Fig. 3). However, during the spin-up period, wind speed and incoming longwave radiation were missing due to a programming error in the AWS. To provide this missing data, wind speed was generated synthetically using the hourly average from the measured data during the evaluation period. Incoming longwave radiation was obtained from the ERA-Interim reanalysis (0.75 • × 0.75 • resolution; Dee et al., 2011), using data from the closest model grid cell after interpolation from 12-hourly to hourly reference points. For the time period where both ERA-Interim and AWS data overlap (20 July-11 August 2008), the mean deviation (MD) and mean absolute deviation (MAD; ERA minus AWS) are ∼ 13 and ∼ 35 W m −2 , with the deviation likely arising due to the difference between modelled and real terrain height of −450 m. Lastly, a rain gauge was not installed at the AWS site in 2008. We therefore used input data from another AWS located 4 km away (denoted as AWS2 in Fig. 2) and assumed that they were representative of conditions at the AWS on the Miage Glacier. The 2008 simulation was intended to coincide with a supplementary field measurement program. Between 20 July and 11 August, surface temperature and the turbulent fluxes of latent and sensible heat were measured. The first field was measured with a CNR1 radiation sensor (Kipp & Zonen, Delft, the Netherlands), while the latter two fluxes were measured by an eddy covariance (EC) station. This comprised a CSAT three-dimensional sonic anemometer and KH2O Krypton Hygrometer (both Campbell Scientific Limited, Shepshed, UK), installed at a height of 2 m above the debris surface. These sensors measured the three components of turbulent wind velocity, virtual temperature and water vapour concentrations at an interval of 50 ms. Raw data were processed using Campbell Scientific OPEC software, which included a "WPL" (Webb-Pearman-Leuning) correction for density effects (Webb et al., 1980) and 30 min averages of the 50 ms scans were stored. The data were filtered for outliers using three times their standard deviation before being used for evaluation . Surface temperature was calculated from the upwelling longwave radiation recorded by the CNR1, using an emissivity of 0.94. The AWS tripod provided a stable platform on the slowly melting glacier surface, although the possibility of tilting of the instrument mast cannot be excluded. These measurements provide a unique data set with which to evaluate the CMB models using direct measurement of turbulence in the surface atmospheric layer above a debris-covered glacier.
However, the 2008 simulation does not contain any phase changes, since the air temperature remained above freezing (cf. Fig. 3a). In order to fully demonstrate the model capabilities, we performed a second simulation from 6 June to 11 October 2011, discarding all but the period of 14 September-11 October as model spin-up time, due to the focus on the influence of phase changes. We focused our analysis on two   freezing events, from 18 to 19 September and 7 to 9 October 2011. Incoming longwave radiation, precipitation and mean wind speed were available hourly from the AWS (forcing data not shown), and measured surface temperature data, estimated from the upwelling longwave radiation recorded by the CNR1, were available for model evaluation.
A final forcing variable for the calculation of the debris surface energy balance, surface pressure, was missing for both the 2008 and 2011 simulations. These data were obtained from the ERA-Interim reanalysis, at 6-hourly temporal resolution, and again from the closest grid cell. A correction was applied for the difference between the real and modelled terrain height using the hypsometric equation, assuming a linear temperature gradient calculated from the AWS and the air temperature on the first model level in the ERA-Interim. For both simulations, the same subsurface layer spacing was used and is provided in Table 2. The englacial temperature profile was initialized at the melting point, since both simulations began in June. Uncertainties in the temperature initialization were addressed by the inclusion of long spin-up periods.
For both CMB-DRY and CMB-RES, we assumed that the debris porosity was a linear function of depth in the debris, decreasing from 40 % at the surface down to 20 % at the debris-ice interface. A range of 19-60 % percent void space by volume was measured on the Miage Glacier, by placing a known volume of surface debris in a graduated bucket and measuring the volume of water required to fill the air spaces (Brock et al., 2006). For this study, we used an upper bound of 40 %, such that the bulk porosity (30 %) was consistent with other reported values for glacier debris (Nicholson and Benn, 2012 bound of 60 % showed that while subdebris ice melt was strongly affected (it decreased by ∼ 17 % in both simulations), the CMB model behaviour and the main results presented in Sect. 3 remained intact. Other physical and thermal properties of the column were either taken from field measurements or specified from values used in previous modelling studies of this glacier (e.g. Reid and Brock, 2010). The porosity value of 20 % in the lowest 1 cm layer in the debris gave a minimum water content of 2 kg m −2 that was imposed only when the subdebris ice was at the melting point. Subdebris ice melt changes by ±1.8 % if the minimum value is removed or doubled in the 2008 simulation. A slope of 7 • at the AWS gives a runoff timescale of 0.92 h −1 . This simple representation of runoff timescales does not consider contributions from upslope regions in the glacier; however, we feel that this is an appropriate first step given that horizontal transport of water within the debris is poorly constrained and no measurements are available. Varying the runoff timescale by ±4 % (equivalent to changing the slope from 4 • to 10 • ) results in small changes in total accumulated mass balance and subdebris ice melt during the summer 2008 simulation, of less than ±0.6 and ±0.4 %, respectively. The results in the transition season of fall 2011 are more sensitive, with changes in these variables of up to ±1.0 and ±2.0 %, respectively.
Finally, although the CMB models are evaluated against a short summer period in 2008 and in fall 2011, they are applicable throughout the annual cycle and to glaciers of any temperature regime, as illustrated in Fig. 1.

Comparison with in situ measurements
The surface temperatures (T sfc ) simulated by CMB-DRY and CMB-RES are in good agreement with measurements for both the 2008 and 2011 simulations (Fig. 4a, d; Table 3). The models tend to underestimate daily maximum temperatures in 2008 and night-time radiative cooling in 2011. However, for both simulations, the models reproduce the diurnal cycle and its variability well. The CMB models also capture the variability of the sensible heat flux (QS), but the simulated magnitude of heat transfer to the overlying atmosphere is greater than reported by the EC station (Fig. 4b).
The overestimation of QS for the CMB-DRY run is in part attributable to the lack of latent heat flux (QL), which means that an average energy loss of ∼ 24 W m −2 is not captured (Fig. 4c; cf. Table 3). CMB-RES has a greatly reduced but still non-negligible bias in QS, again, in part, because evaporative cooling is underestimated, by ∼ 6 W m −2 . The smaller simulated latent heat flux compared with the EC data results from the approach used to estimate surface vapour pressure (cf. Sect. 2.2), which produces an average gradient of only −0.5 hPa m −1 between the surface and overlying air (Fig. 5a).

Modelling insights from the 2008 simulation
In total, the influence of the reservoir parameterization on the accumulated mass balance between 20 July and 11 August 2008 is small, increasing from −241.0 kg m −2 in CMB-DRY to −250.6 kg m −2 in CMB-RES (Fig. 5b). These values are equivalent to an ablation rate of approximately 11 mm w.e. d −1 (w.e. water equivalent), which is in order of magnitude agreement with the value of 22 mm w.e. d −1 reported by Fyffe et al. (2012) for the Miage Glacier, based on the entire ablation seasons of 2010 and 2011. The mass fluxes underlying the simulated mass balance signal are determined by the surface energy balance, whose daily-mean components are shown in Fig. 6a for CMB-RES. Energy receipt, mainly through net shortwave radiation, is generally counteracted by energy losses though net longwave radiation, heat conduction (QC), and the turbulent fluxes QL and QS. The heat flux to the debris surface from precipitation (QPRC) has an average value of −12.5 W m −2 during rainfall events. However, since the precipitation temperature is assumed to be the same as T air , QPRC is a stronger energy sink for daytime rainfall. These energy fluxes produce ablation that is dominated by subdebris ice melt and evaporation over the evaluation period ( Fig. 6b; Table 4). Surface melt, refreeze, sublimation and deposition are zero, since there is no solid precipitation and both the debris surface and internal temperatures remain above the melting point.
Compared with CMB-DRY, CMB-RES simulates slightly lower daytime debris-surface temperatures, as a result of heat extraction by QL (cf. Fig. 4a, Table 3). Energy transfer to the debris-ice interface is therefore also lower, contributing to   a small reduction in subdebris ice melt, of 7.5 kg m −2 (Table 4). However, the reduction in melt is more than compensated by surface vapour fluxes, with a total of 17.3 kg m −2 of evaporation over the evaluation period. Evaporation dom-inates during the day (95 % of the total), while smaller amounts of condensation occur mainly at night (64 %) or in the early morning. Water accumulates in the supraglacial debris after rainfall events and is then removed, mainly by horizontal drainage but also by evaporation (Fig. 7). The total accumulated mass balance is negligibly altered if changes in debris-water content are considered in addition to surface vapour fluxes. Both models treat the physical properties of the debris layer -thermal conductivity, density, and specific heat capacity -as functions of depth. Figure 8a-c show their variation with depth for "dry" conditions, when there is no significant debris-water storage, and for "wet" conditions, when there is significant water present, as a result of rainfall. "Dry" conditions prevail, comprising 76 % of the evaluation period ( Fig. 8d-f), under which, as the porosity decreases with depth, the debris thermal conductivity and density increase while the specific heat capacity decreases. The debris physical properties in CMB-DRY and CMB-RES are the same, with the exception of the bottom layer adjacent to the debrisice interface, which remains fully saturated in CMB-RES, as a result of the moisture source term described in Sect. 2.2. Water present in this layer acts to increase all three properties compared with CMB-DRY. Rainfall events and the associated moisture storage extend this influence upwards through the debris layer, with a significant alteration to the fully saturated layers (spanning the depth between 20 and 23 cm for  Figure 6. CMB-RES values for (a) daily mean energy fluxes over the evaluation period (W m −2 ). The grey curve is net shortwave radiation, the black curve is net longwave radiation, and the grey dots show daily-mean surface albedo, which remains constant at the debris value because there is no solid precipitation. (b) Daily total mass fluxes (kg m −2 ). Maximum daily values of evaporation and condensation are 1.4 and 0.02 kg m −2 , respectively, although the latter flux is not visible. Note that while daily-accumulated rainfall is shown (purple asterisks), it is not technically a mass flux, since the mass balance calculation in CMB-RES does not account for debris-water storage. Rather, this field is plotted to show its correspondence with other fields, such as net shortwave radiation.
the "wet" sample time slice) and a smaller effect on the partially saturated layer (at a depth of 19 cm). The debrisspecific heat capacity is the most strongly affected physical property, since the value of water is approximately four times that of air (4181 vs. 1005 J kg −1 K −1 ).
The effective thermal diffusivity of the debris is inversely proportional to the specific heat capacity and the debris density. Increases in both these quantities, but particularly in the former, reduce heat diffusion over affected layers compared with CMB-DRY. Therefore, in combination with heat extraction by QL, the change in subsurface physical properties reduces the amplitude and depth penetration of the diurnal temperature cycle in the debris layer (Fig. 9). Fluctuations in the magnitude of QL have a correlation coefficient of 0.78 with the temperature difference between CMB-RES and CMB-DRY in the top 6 cm of the debris, while reductions in the effective thermal diffusivity have a correlation coefficient of 0.6 with the temperature difference in the bottom 6 cm.

Impact of phase changes in the 2011 simulation
Two freezing events occur during the 2011 simulation, between 18 September 23:00 LT (local time) and 19 September 14:00 LT and between 7 October 09:00 LT and 9 October 09:00 LT, at the end of two precipitation events with subzero air temperatures (cf. Fig. 4d). Energy and mass fluxes for this simulation are summarized in Table 5. Net longwave and shortwave radiation are reduced, due to cooler surface temperatures and to small amounts of snowfall that increase the surface albedo (Fig. 10a). Rapid melt of the thin overlying snow cover (< 0.5 cm) and infiltration of rainfall at the beginning of the precipitation events provide the source water for refreeze in the debris (Figs. 10b,  11a). During the first event, a maximum of 1.0 kg m −2 of ice is produced, which persists in the basal debris layer for a further 3 days after the last time step with refreeze. In the , and (c) specific heat capacity [ J kg −1 K −1 ], shown for CMB-DRY in grey-unfilled circles and for CMB-RES in both black-filled circles ("dry" time slice) and blue asterisks ("wet" time slice). Time series of bulk values for these same properties are shown in panels (d-f) for CMB-RES in blue and CMB-DRY in grey. The locations of the "dry" and "wet" time slices are indicated by the first (solid grey) and second (dashed grey) reference lines on the x-axis, respectively. , and (c) specific heat capacity (J kg −1 K −1 ), shown for CMB-DRY in grey-unfilled circles and for CMB-RES in both black-filled circles ("dry" time slice) and blue asterisks ("wet" time slice). Time series of bulk values for these same properties are shown in panels (d-f) for CMB-RES in blue and CMB-DRY in grey. The locations of the "dry" and "wet" time slices are indicated by the first (solid grey) and second (dashed grey) reference lines on the x axis, respectively. second event, the debris ice content reaches 1.4 kg m −2 , and does not melt away before the end of the simulation.
The bulk presence of liquid water and ice in the debris layer influences the vertical temperature profile in two competing ways (Fig. 11b-d). Latent heat release due to refreezing warms the subsurface, on average by 0.3 K but exceeding 0.7 K for the hourly time steps with the greatest refreeze. However, the presence of ice in saturated basal layers constrains the debris temperature to the melting point. In combination with a reduction in the effective thermal diffusivity of saturated layers, the modulation of debris temperature results in a decrease in subdebris ice melt of 7.0 % in CMB-RES compared to CMB-DRY.
The accumulated mass balance between 14 September-11 October 2011 is −172.4 kg m −2 for CMB-DRY and −168.8 kg m −2 for CMB-RES. Changes in water and ice storage again have a negligible impact on simulated mass balance, resulting in a further ablation of 0.2 kg m −2 . Thus, for the fall transition season, surface vapour fluxes do not compensate for the reduction in subdebris ice melt due to the thermodynamic influence of ice in the debris. However, considering the same summer period in 2011 as in 2008 (20 July-11 August), the percent changes in accumulated mass balance and subdebris ice melt are +4.0 and −3.2 %, respectively, consistent with the findings of the 2008 simulation. Therefore, the influence of the reservoir parameterization varies seasonally.  raindrops interfering with the path of the sonic anemometer (e.g. Aubinet et al., 2012). However, of the 15 occurrences of spikes greater than ±100 W m −2 in the EC data, only two occur during or within 1 hour of precipitation. Given previously reported large QL values, of up to −800 W m −2 during rainfall events on heated debris , neglecting QL in a surface energy balance calculation can be inappropriate, and under certain meteorological conditions is likely to have a significant impact on the calculated energy fluxes.
The difference in accumulated mass balance between CMB-RES and CMB-DRY is relatively small, for a point application in this configuration. However, the daily mean evaporation rate was ∼ 0.9 mm w.e. in 2008 (June-early August) and ∼ 0.6 mm w.e. in 2011 (June-September), which is comparable to values reported for clean glaciers (e.g. Kaser, 1982). Scaled up to a larger debris-covered area, evaporation would represent a significant mass flux. Furthermore, the presence of debris ice, even in small amounts, has an important thermodynamic influence by suppressing subdebris ice melt, with implications for dry simulations of debris-covered glaciers in or close to transition seasons.
The simulated QL and surface vapour fluxes depend on the estimate of the surface vapour pressure, which is an important source of uncertainty in the CMB-RES model. In unsaturated soil sciences, the relative humidity is often treated as an exponential function of the liquid water pressure in the pore space using the thermodynamic relationship of Edlefsen and Anderson (1943) (e.g. Wilson et al., 1994;Karra et al., 2014). However, testing an exponential relationship with the moisture content of the debris in CMB-RES resulted in strong biases in QL (MD = 28; MAD = 96 W m −2 ) and a shift from QL as an energy sink to a gain, which was inconsistent with the EC data. For simplicity, we employed a linear approach, and there may be some support for this treatment in coarser texture soil, as Yeh et al. (2008) found that the effective degree of saturation in sand decreased approximately linearly in the top 2 m above the water table.
In reality, water vapour fluxes occur at the saturated horizon, either at the surface or within the debris layer. However, in the 2008 simulation, the mean depth of the saturated horizon was 21.5 cm, where the proximity of glacier ice damped temperature fluctuations and constrained the mean temperature to ∼ 275 K. Therefore, computing vapour fluxes at this level produced a very small latent heat flux, of −3.1 W m −2 on average, that was also not in agreement with the EC data. CMB-RES likely provides an underestimate of the simulated location of the saturated horizon, since capillary action was not taken into account. For fine gravel soils (grain size of 2-5 mm), capillary rise is on the order of a few centimetres (Lohman, 1972), while for coarser, poorly sorted glacier debris, the effect may be smaller. Underestimation of the height of the saturated horizon, and therefore of both the debris temperature and the saturation vapour pressure, is consistent with the small latent heat flux when vapour fluxes are computed at this level. As a part of future work, there is a need to accurately compute the vapour fluxes at the level of the saturated horizon.
In addition to neglecting capillary action, CMB-RES also does not account for many internal physical processes that have been highlighted in unsaturated soil sciences, including   water vapour flow due to gradients in concentration and temperature, liquid water flow in response to hydraulic gradients, volume changes due to changes in the degree of saturation (e.g. Sheng, 2011), deposition of water vapour and its contribution to the formation of thin ice lenses (e.g. Karra et al., 2014), and heat or moisture advection as a result of airflow (e.g. Zeng et al., 2011). However, incorporation of these processes into CMB-RES is currently limited by a lack of appropriate evaluation data. Instead, we focus on including processes related to phase changes, which have been demonstrated to have an impact on the subsurface tempera- ture field and ablation rate (Reznichenko et al., 2010;Nicholson and Benn, 2012). As a part of future work, CMB-RES could be improved by distinguishing the location of debris ice and water separately within saturated layers, thus potentially improving the simulated debris temperature profiles, as the melting point constraint would only be applied to saturated layers containing ice. The magnitude of QS is sensitive to the choice of debris thickness, which was selected to be 0.23 m in this study based on a point measurement. However, the turbulent fluxes measured by the EC station respond to a larger area, with a variable and unknown debris thickness that likely ranges between 20 and 30 cm. The agreement between measured and modelled QS in 2008 is improved if the debris thickness in the models is reduced slightly. For example, using a thickness of 20 cm reduces the MD and the MAD by ∼ 7 W m −2 , for both model versions. Investigating additional causes of discrepancies between modelled QS and that measured by the EC is not directly related to the inclusion of moisture in CMB-RES and is reserved for future work.
There are no ablation measurements available for either of the two simulation periods. To examine the general behaviour of the CMB models, the 2008 simulation was repeated with debris thicknesses of 1-20 cm, holding the subdebris ice depth constant and scaling the minimum debris-water content as 3 % of the reservoir capacity (consistent with the 23 cm simulation; Fig. 12). Total column melt is suppressed for all debris thicknesses compared with the clean-ice-melt rate, with less melt in CMB-RES than CMB-DRY due to heat extraction by QL and the reduced thermal diffusivity discussed in Sect. 3.2. Therefore, the www.the-cryosphere.net/8/1429/2014/ The Cryosphere, 8, 1429-1444, 2014 CMB models do not reproduce the typical Østrem curve, wherein melt is enhanced below a critical debris thickness that ranges between 1.5 and 5 cm (e.g. Loomis, 1970;Fujii, 1977;Inoue and Yoshida, 1980;Mattson et al., 1993) and suppressed above this value. The rising limb of the Østrem curve is not reproduced for several reasons. First, in the clean-ice and thinly debris-covered simulations, lower nighttime air temperatures in the beginning of the evaluation period (20-24 July 2008; cf. Fig. 4a) produce freezing events that cool the subsurface. Averaged over the entire evaluation period, a non-negligible amount of energy is expended to warm the ice column as a result. For example, in the clean-ice simulation, this heat flux amounts to 3.7 W m −2 . For CMB-RES (CMB-DRY) with debris thicknesses of 1 and 2 cm, the average energy required is 4.4 (5.3) and 3.1 (3.5) W m −2 , respectively. In addition, subzero englacial temperatures in the clean-ice simulation are eradicated more quickly, since penetrating shortwave radiation is considered. Finally, other processes that are not treated in the CMB models may be important to fully reproduce the rising limb of the Østrem curve, such as (1) changes in the surface albedo as the debris cover becomes more continuous, as in the albedo "patchiness" scheme introduced by Reid and Brock (2010), and (2) wind-driven evaporation inside the debris layer (Evatt et al., working paper, 2014).

Conclusions
In this paper, we introduce a new model for the surface energy balance and CMB of debris-covered glaciers that includes surface vapour fluxes and a reservoir parameterization for moisture infiltration and phase changes. Although the parameterization is a simplification of the complex moist physics of debris, our model is a novel attempt to treat moisture within glacier debris cover, and one that permits two important advances: (1) it incorporates the effects of ice and water on the physical and thermal properties of the debris and therefore on ice ablation, and (2) it includes an estimate of the moisture exchanges between the surface and the atmosphere. The inclusion of the water vapour flux opens up avenues of future research. For example, distributed simulations are required to more rigorously investigate relevant scientific questions about debris-covered glaciers, such as projecting their behaviour and runoff under changing climate conditions. A key constraint in performing such simulations is obtaining forcing data, since the highly heterogeneous surface of debris-covered glaciers makes the spatial distribution of air temperature and winds uncertain. Current approaches, employing elevation-based extrapolation, appear to be inadequate (Reid et al., 2012). Interactive coupling with a high-resolution atmospheric model provides one solution; however, the conventional modelling approach would introduce errors due to the absence of moisture exchange be-tween the surface and the atmosphere. In incorporating that flux, CMB-RES is a step toward more precisely computing glacier-atmosphere feedbacks within coupled surface-andatmosphere modelling schemes and more accurately predicting alterations in freshwater budgets and other potential impacts of glacier change.