Sensitivity of a distributed temperature-radiation index melt model based on AWS observations and surface energy balance fluxes , Hurd Peninsula glaciers , Livingston Island , Antarctica

We use an automatic weather station and surface mass balance dataset spanning four melt seasons collected on Hurd Peninsula Glaciers, South Shetland Islands, to investigate the point surface energy balance, to determine the absolute and relative contribution of the various energy fluxes acting on the glacier surface and to estimate the sensitivity of melt to ambient temperature changes. Long-wave incoming radiation is the main energy source for melt, while short-wave radiation is the most important flux controlling the variation of both seasonal and daily mean surface energy balance. Short-wave and long-wave radiation fluxes do, in general, balance each other, resulting in a high correspondence between daily mean net radiation flux and available melt energy flux. We calibrate a distributed melt model driven by air temperature and an expression for the incoming short-wave radiation. The model is calibrated with the data from one of the melt seasons and validated with the data of the three remaining seasons. The model results deviate at most 140 mm w.e. from the corresponding observations using the glaciological method. The model is very sensitive to changes in ambient temperature: a 0.5 C increase results in 56 % higher melt rates.


Background
Retreating and thinning glaciers have come into sharp focus in relation to increased atmospheric temperatures attributed to anthropogenic greenhouse emissions.Changes in meteorological conditions and the associated changes in air temperature will, in different ways and to different extent, affect the various fluxes providing energy for heating and melting a glacier surface as part of an intricate system involving several feedback mechanisms (Ohmura, 2001).Surface energy balance (SEB) calculations using automatic weather station (AWS) data from glaciers aim to separate and quantify the contributing fluxes, and form a basis for understanding the coupling between meteorological conditions and glacier melt.As a prognostic tool for the response in energy fluxes, and eventually melt rates, to perturbations in meteorological conditions, SEB models have the advantage of being physically based but the disadvantage of involving a complicated extrapolation procedure to distribute the fluxes over the glacier surface (Hock, 2005).To overcome this complexity, simpler temperature-index models, based on the coupling between energy fluxes and temperature, are widely used.They perform best in environments where long-wave radiation and sensible heat are the dominating energy sources, as those fluxes are strongly coupled to temperature (Ohmura, 2001), while in environments dominated by solar radiation the model performance is reduced (Sicart et al., 2008).In many glaciated areas of the world, long-wave radiation is the major energy source for melt, but short-wave radiation plays a dominant role for the diurnal, daily and inter-seasonal variation of energy available for heat and melt (Greuell and Genthon, 2004).Therefore models including representations of a combination of temperature and short-wave radiation are commonly used (e.g.Hock, 1999;Pellicciotti, et al., 2005;Schuler et al., 2005;Möller et al., 2011).
An increasing number of studies based on AWS data combined with glaciological surface mass balance data are continuously improving the understanding of glacier SEB and its geographical differences (e.g.Andreasssen et al., 2008;Giesen et al., 2008;Hulth et al., 2010).In this paper we present an AWS record extending over four melt seasons and its associated SEB from a region scarcely represented in the literature, namely the Hurd Peninsula on Livingston Island, South Shetland Islands, an archipelago parallel to the northwestern extreme of the Antarctic Peninsula (AP) (Fig. 1).The AP region has experienced a remarkable warming during the last decades (Vaughan et al., 2003;Turner et al., 2009).The 41-year long temperature record from Bellinghausen research station on King George Island (KGI), located about 90 km to the NE of Hurd Peninsula, shows a 0.25 • C per decade warming (http://www.nerc-bas.ac.uk/icd/ gjma/temps.html).The impact of this regional warming on surface mass balance and melt rates has been analysed by several authors (e.g.Braun and Hock, 2004;van den Broeke, 2005;Turner et al., 2005Turner et al., , 2009;;Vaughan, 2006).Additionally, atmospheric warming, together with warmer ocean temperatures, have been pointed out as drivers, through different physical processes, of the disintegration of some ice shelves on the northeastern coast of the AP (MacAyeal et al., 2003;Shepherd et al., 2003;van den Broeke, 2005;Cook and Vaughan, 2010), with subsequent acceleration of the inland glaciers feeding the ice shelves (Rott et al., 1996;Rignot et al., 2004;Scambos et al., 2004), and also of the widespread retreat of marine glacier fronts of the AP over the past halfcentury (Cook et al., 2005).An overall tendency of retreating ice fronts has also been observed in studies analysing, over the period 1986-2002, both marine-terminating and landterminating glaciers in the region (Rau et al., 2004).A widespread acceleration trend of glaciers on the AP west coast has been observed as well from repeated flow rate measurements within 1992-2005, and attributed to a dynamic response to frontal thinning (Pritchard and Vaughan, 2007).On the northeastern side of the AP, however, the rate of recession of ice-shelf tributary glaciers has slowed down markedly during the last decade as compared to the previous one (Davies et al., 2011).All of these observations, and their underlying physical processes, are key for understanding the contribution of grounded ice losses from this region to sea level rise that has been currently estimated as 0.22 ± 0.16 mm a −1 (Hock et al., 2009).The Hurd Peninsula glaciers surface mass balance record (Navarro et al., 2012) is the only cur-rently ongoing mass balance programme from the AP with both summer and winter balances listed at the World Glacier Monitoring Service, thus providing a particularly suitable field data source for studying the coupling between meteorological conditions and mass balance.
We use the AWS record to briefly outline the meteorological conditions at the site and to quantify the absolute fluxes and their relative contribution to melt.We also apply a distributed temperature-radiation index melt model, which is calibrated with in-situ surface mass balance data obtained using the glaciological method (Navarro et al., 2012).We compare the melt rates obtained from SEB and temperatureindex model at the point scale, and use the four years of AWS data to compare the distributed modelled and measured melts.We use the calibrated model to analyse its sensitivity to changes in meteorological conditions.Furthermore, we present a novel approach to handle sub-daily albedo data that reduces the diurnal variations that arise due to measurement deficiencies but keeps the diurnal variations due to changes in cloud cover.

Physical setting
The Hurd Peninsula ice cap (62 • 39-42 S, 60 • 19-25 W; Fig. 1) covers an area of about 13.5 km 2 and spans an altitude range from sea level to about 330 m a.s.l.It is partly surrounded by mountains that reach elevations between 250 and 400 m.Hurd Peninsula ice cap can be subdivided into two main glacier systems: Johnsons glacier, a tidewater glacier, flowing north-westwards that terminates on an ice cliff about 50 m in height a.s.l.extending approximately 500 m along the coast, and Hurd glacier, terminating on land and mainly flowing south-westwards with ice thickness tapering to zero in the snouts of Sally Rocks, Argentina and Las Palmas.Two other basins are also part of the ice cap, both flowing eastwards towards False Bay.However, these are not included in this study because they are heavily crevassed icefalls preventing field measurements.The local ice divide separating Johnsons and Hurd lies between 250 and 330 m a.s.l.Hurd glacier has an average surface slope of about 3 • , though the small westward flowing glacier tongues Argentina and Las Palmas are much steeper (approx.13 • ).Typical slopes for Johnsons glacier range between 10 • in its northern areas and 6 • in the southern ones.The average ice thickness of the ice cap, determined from groundpenetrating radar data in 2000/2001, was 93.6 ± 2.5 m, with maximum values about 200 m in the accumulation area of Hurd glacier (Navarro et al., 2009).The ice surface velocities of Johnsons glacier increase downstream from the ice divide, reaching annual-averaged values up to 65 m a −1 at the fastest part of the calving front (Otero et al., 2010), while the largest ice velocities for Hurd glacier are typically about 5 m a −1 (Otero, 2008).The average accumulation area ratio 2002-2011 was 44 ± 24 % for Hurd Glacier and 61 ± 21 % for Johnsons Glacier, with the range being equal to one standard deviation.Similar figures for the equilibrium line altitude were 228 ± 57 m a.s.l. and 187 ± 37 m a.s.l., respectively (Navarro et al., 2012).The average annual surface mass balance over the same time period has been slightly negative for Hurd glacier (-150 mm w.e.) and slightly positive for Johnsons glacier (50 mm w.e.).This is explained by lower accumulation rates on Hurd glacier attributed to snow redistribution by wind, together with slightly higher ablation rates, due to Hurd's hypsometry, which shows a much larger share of area at the lowermost altitudes (<100 m) as compared to Johnsons (Navarro et al., 2012).The Hurd Peninsula ice cap is a polythermal ice mass, showing an upper layer of cold ice, several tens of metres thick, in the ablation area (Molina et al., 2007;Navarro et al., 2009).The cold layer is rather uniform in Hurd glacier while in Johnsons glacier the layer is more patchy.
The Spanish Antarctic Station Juan Carlos I (JCI) is located very close to Hurd Peninsula ice cap (Fig. 1).Meteorological measurements are maintained in JCI all year round by an AWS, and are complemented by manual meteorological observations during the summer.In December 2006, an AWS was installed on the upper ablation area of Johnsons Glacier (Fig. 1).The latter provides the main meteorological data source used in this paper and will be further described in the next section.
The Hurd Peninsula ice cap is subjected to the maritime climate of the western AP region, with some particularities due to local conditions.The orography of Livingston Island alters the regional prevailing northwesterly dominating wind direction.At JCI and Johnsons Glacier the predominant wind direction is from NNE, followed by SSW (JCI), ENE (Johnsons Glacier).Average wind speeds are about 4 m s −1 , but The cloudiness is high, with an average of 6/8 and, consequently, insolation is small, with 2 h day −1 of average insolation during summer and spring, though the cloud-free days during such seasons show a high irradiance.The relative humidity is very high, with average values above 80 % at JCI and 90 % on the glacier (unpublished data from Agencia Estatal de Meteorología, AEMET).
On the glacier, mass gain is dominated by direct snowfall and wind redistribution of snow, without any contribution from snow avalanches.The glacier ice hardly receives any debris from the surrounding mountains, except at the lowest elevations of its outlets.Tephra layers from the recent eruptions of the neighbouring Deception Island, however, are a common feature of these glaciers (Pallàs et al., 2001).

AWS data record
The meteorological data used for SEB analysis and as input to the temperature-radiation index melt model was obtained from the AWS located on Johnsons glacier (62 The data sampling and averaging intervals of the AWS logger differ between the summer seasons, when regular visits to the AWS are done from the nearby JCI research station, and the unmanned winter seasons.The majority of the data corresponds to the summer seasons, but still a portion was collected during winter mode of operation.Thus, temporal averages of the meteorological parameters are based on different number of readings.In case of missing data due to gaps or preposterous registrations, hourly averages are calculated down to the representation of at least one reading during the time interval; otherwise linear interpolation between the adjacent time steps was made for gaps up to four hours.A gap in the temperature and relative humidity records due to sensor malfunctioning, spanning the period 16-24 December 2009, was replaced with temperature data from a T107 thermistor installed on the AWS and with RH data from JCI. Offsets calculated from data two days before and after the gap were applied.The temperature sensor was unventilated, which at high short-wave radiation fluxes and low wind speeds exaggerates temperatures and we therefore applied the correction to overheating suggested by Smeets (2006) on all air temperature data.

Correction of short-wave radiation and albedo
Large pseudo-diurnal variations in albedo (Fig. 2), α, may arise from the poor cosine response of the pyranometeres at high zenith angles and from an instrument set-up non-parallel to the glacier surface.The bias is mainly introduced by the upward facing sensor as incoming solar radiation flux is less diffuse compared to the outgoing (van den Broeke, 2004).As changes of surface conditions affecting α are of a gradual nature, with the exception of snow falls, van den Broeke ( 2004) presented an approach to improve the accuracy of the incom- ing short-wave fluxes by filtering out diurnal variations using a 24-hour running mean "accumulated" albedo, α d (Eqs. 1 and 2).This approach was subsequently adapted by others (Andreassen et al., 2008;Giesen et al., 2008).α d is calculated as the ratio between the sums of the instant (measured) outgoing flux, S r , and the measured incoming solar radiation flux, S i , with the timing of measurement as midpoint.A corrected incoming short-wave flux, S d , is thereafter calculated.
where t is time in hours (h) Equation ( 2) produces a smooth diurnal signal of α d that leaves out the extreme values arisen from measurement shortcomings, which normally occur at high zenith angles.As α d is based on the sums of the short-wave fluxes, the low fluxes occurring at high zenith angles reduce the impact of the measurement deficiencies on α d and S d .However, this approach filters out the rapid variations in α associated with changes in cloud cover also during midday hours when shortwave fluxes are high, which may have significant impact on SEB.The rapid variations are an effect of clouds absorbing radiation to a higher degree in the near-infrared spectra, leaving a higher portion of the visible wave lengths to reach the ground, for which the snow reflectance is higher.This effect is further enhanced by multiple reflections between the surface and the cloud base (Warren, 1982;Cutler and Munro, 1996).With Eqs. ( 3) and ( 4) we present and use an extension of Eq. ( 1) which yields a corrected albedo, α c , and corrected short-wave flux, S c , that restore the effect of clouds reasonably well where λ is the slope of the linear relation between the change in instant albedo, α i , and the instant fraction of potential top of atmosphere radiation that reaches the ground, θ i , between two subsequent recording time steps (Fig. 3).The term restoring the effect of clouds is proportional to the difference between θ i and its corresponding 24-hour running mean, θ d , and scaled by λ. θ i is calculated as where the top of atmosphere radiation, I toa , is given by I S is the solar constant (1366 W m −2 ), d and d m are the instant and mean Sun to Earth distance, and Z is the local solar zenith angle.A freely available Matlab code based on Reda and Andreas (2008) was used to calculate the required parameters.An iterative procedure is needed for solving Eq. ( 3) and subsequently Eq. (4) as Eq. ( 1) uses S c as input.The result normally stabilizes after less than five iterations.Regression analysis to obtain the linear slope λ was applied to each of the four seasonal data sets on all data when zenith angle was 80 • or less.The resulting four λ values were all within -0.09 ± 0.01 with corresponding correlation coefficients between -0.59 and -0.48.Consequently, λ was statically set to -0.09 in Eq. ( 3).The daily mean α will be the same regardless of using S i , S d or S c as incoming short-wave flux.Jonsell et al. (2003) found only a limited effect of clouds on α during ice conditions, probably because the multiple reflection of radiation will be limited due to lower α.The location of the AWS was snow covered during the period of data collection, with the exception of a few days at the end of the melt season 2006/07, but still with α exceeding 0.55, and therefore we applied Eqs. ( 3) and (4) to the full data set.The rapid shifts in α c due to sudden changes in S c are well represented.During early morning and late afternoon hours an offset towards higher α c is present.This is due to reduced θ i caused by the longer atmospheric pathway and increased scattering of solar radiation at high zenith angles, and thus it is not primarily an effect of clouds.However, during the morning and afternoon hours S c is small and the impact on SEB is limited, while the correction by Eqs. ( 3) and (4) during midday hours, when S c is large, will have greater importance for SEB.

Energy balance calculation
We calculate SEB, using hourly mean AWS data, to investigate the absolute and relative contribution of the energy fluxes involved in glacier heating and melt.The energy balance condition at the glacier surface is represented by the equation where S n is the net short-wave flux, L in and L out are the incoming and outgoing long-wave fluxes, respectively, H is the sensible heat flux, and E is the latent heat flux, which summed up are expressed as the atmospheric energy fluxes, A, and must balance the sum of ground heat flux, G, and melt energy flux, M. A flux is considered to be positive when directed towards the surface.We do not quantify the energy from the sensible heat of rain fall due to lack of precipitation measurements on the glacier, but this term is generally small (Hock, 2005).The long-wave radiation components are, in contrast to their short-wave counterparts, treated separately.This division mirrors better their different physical nature and their different response to atmospheric and surface conditions (Ohmura, 2001).We use the bulk aerodynamic method based on Monin-Obukhov length following Hock and Holmgren (2005).No field measurements of the roughness length for momentum (Z 0w ) are available for the site.Although Z 0w is expected to vary with time (Brock et al., 2000;Hock, 2005), we use Z 0w as a constant-in-time model parameter to fit the SEB to the ablation rates estimated from the surface lowering registered by the ultra-sonic ranger.Conversion from surface lowering to mass loss was made by using depth-density relations from snow pits dug in the area.The roughness lengths for heat, Z 0T , and moisture, Z 0e , are treated as functions of Z 0w following Andreas (1987).For stable stratification we use the stability functions by Beljaars and Holtslag (1991), and for unstable conditions the expression by Panofsky and Dutton (1984).The best fit was obtained using Z 0w =1 × 10 −4 m, yielding typical seasonal averages and standard deviations of (1.3 ± 0.06) ×10 −4 m for Z 0T and (1.6 ± 0.10) ×10 −4 m for Z 0e .
Surface temperature, T s , in • C, is calculated from L out assuming emissivity of unity www.the-cryosphere.net/6/539/2012/The Cryosphere, 6, 539-552, 2012 where σ is the Stefan-Boltzmann constant (5.67 When E is negative, the phase change is considered as sublimation, while when it is positive it is considered as condensation or as resublimation, depending on whether T s equals or is lower than 0 • C. G could not be captured by any of the sensors of the AWS.On a perfectly temperate glacier G is by definition zero.The Hurd Peninsula glaciers are of polythermal type (Navarro at al., 2009) and during the melt season G will be directed to heat the cold surface layer.From other studies (see table in Ohmura, 2001) we conclude that the average G will be limited to a few W m −2 , but will occasionally be high due to the release of latent heat of refreezing of melt water upon the first onset of melt in the summer, and during subsequent recommences of melt-refreeze cycles.During the ablation season, melting conditions at the Hurd Peninsula glaciers are repeatedly interrupted by refreezing events.In other modelling approaches, the quantification of the energy compensation needed to reach melting conditions after a period of freezing and cooling of the surface (negative A) differs from full compensation (van de Wal and Russel, 1994) to no compensation (Hock and Holmgren, 2005).The latter means that, as soon as A turns positive, melting is resumed.Hock and Holmgren (2005) argued that the positive compensation flux will just be a part of the accumulated negative flux, because the latent heat from percolating refreezing melt water will make a significant contribution to the reheating of the cooled volume when melt conditions resume at the surface.
We consider the surface to be melting when T s > −0.5 • C and we set G constant to 5 W m −2 during melting.The temperature offset from zero degrees is made to allow for melting at depths below the surface skin layer that T s is representing.Sub-surface melting at sub-zero T s can occur because of the ability of short-wave radiation to penetrate into the snow cover, where the thermal conductivity is low (Koh and Jordan, 1995;Kuipers et al., 2009).We consider that the surface is refreezing and cooling when A turns negative, and to be heated to resume to melt conditions when A is positive and T s < −0.5 • C. Thus, G balances A during non-melting conditions.In practice this implies that the refreezing flux (negative A) is compensated by the positive flux of A when T s < −0.5 • C. Melt rates calculated via the energy balance are converted to water equivalent (w.e.) using the latent heat of fusion (334 × 10 3 J kg −1 ).
Following Sicart et al. (2008), we describe the relative contribution of the energy fluxes to the variation of A as where r x is the correlation coefficient of flux x to A, and σ x, σ A are their respective standard deviations.The sum of ρ x for all fluxes contributing to A is one.

Temperature-radiation index model description
We performed the surface melt modelling and the sensitivity analysis to perturbations in meteorological conditions using a grid-based distributed temperature-radiation index model following Hock (1999) where M xy is the modelled melt at a specific grid cell.The melt factor, m, and the surface specific radiation factor, δ snow/ice , are model parameters obtained via calibration.The temperature at a specific grid cell, T xy , is given by the AWS temperature record with an offset based on altitude difference and air temperature slope lapse rate (dT/dZ).Braun and Hock (2004) showed the dependence of dT/dZ and melt rates with the synoptic weather pattern and recommended to avoid the use of a constant lapse rate.We did not find conclusive correspondence between dT/dZ and wind direction, wind speed or θ d that could be indicative of different weather patterns and which could be used to differentiate dT/dZ among our model runs for distinct seasons.Consequently, we applied a constant dT/dZ based on the temperature difference between JCI and AWS during season 2009/10, that yielded -7.0 ± 0.5 K km −1 , with the range being one standard deviation.This value is within the normal range of reported lapse rates in the AP region (Braun and Hock, 2004).In the evaluation of the calibration run of the model, we analyse the model sensitivity to perturbations of a constant dT/dZ.D is a representation of the direct solar radiation flux and is distributed over the glacier surface depending on the angle of incidence of the solar beam, and scaled from clear-sky conditions using θ i where Z is the solar zenith angle, is the solar azimuth, β is the surface slope angle and is the surface aspect.We considered Z and to be spatially constant in the study area, while we calculated β and from a digital elevation model (DEM).The DEM has a grid cell resolution of 25 m and was obtained by kriging interpolation from 852 points surveyed in 2000/01.The surrounding rock topography was digitized from Servicio Geográfico del Ejército 1: 25000 map for Hurd Peninsula (SGE, 1991).
A shaded grid cell implies that D is set to zero.We calculated the shading of grid cells by applying standard geometry on the DEM.When the AWS is shaded in the sunset, the preceding non-shaded value of θ i is used to calculate nonshaded grid cells, while the subsequent non-shaded value is used in the sunrise.
We initialized the model with a distributed snow cover thickness which is continuously reduced with the daily The Cryosphere, 6, 539-552, 2012 www.the-cryosphere.net/6/539/2012/modelled melt and increased in case snow fall was registered by the ultra-sonic ranger.Snow fall was added only to grid cells with negative temperature and rain was treated as runoff, i.e. not contributing to accumulation.We distributed the snow fall following the mean accumulation pattern obtained from determining the surface mass balance by the glaciological method for the period 2001-2010 (Navarro et al., 2012).We obtained the initial snow cover grid of each season by interpolating snow depth measurements conducted simultaneously to the first stake height reading, and converting it into mass loss using the density measurements from three snow pits conducted at the beginning of the melt season.

Temperature-radiation index model calibration
The model calibration is based on comparing the surface mass balance, b sfc , of a number of stakes obtained by the glaciological method (Navarro et al., 2012) and the model results for the grid cells corresponding to the stake positions.As Eq. ( 10) is set to compute melt during a specific period and not b sfc , snow fall rates registered by the ultrasonic ranger (Fig. 4) were distributed over the glacier and subtracted from the modelled melt to produce comparable figures.We assumed a snow density of 300 kg m −3 for fresh snow.
We converted the stake height differences into water equivalents using the density measurements from three snow pits at the beginning and the end of the melt season.As densification of the snow pack, due to refreezing of melt water, is expected during the melt season, the total mass loss at each stake at the end of the calibration period is given by the beginning of the season mass minus the end of the season mass.On average, the end of melt season densities were 14 % higher.Snow density is however expected to show local variations and introduces a source of error.Based on extensive density measurements during 2004-2010 on Johnsons and Hurd glaciers (Navarro et al., 2012), the relative double standard deviation of all measured snow densities was 13 % and provides an indication on the error involved.Where the model indicates ice at the surface, we used a density of 900 kg m −3 to convert into water equivalents.
Among the four seasons with data available from the AWS, we chose the period 29/11/2008-10/02/2009 as calibration period, as it provided the best combination of the following criteria: number of stakes, length of calibration period, high ablation rates, and number of stakes on ice surface at the end of the period.For the calibration period 2008/09 the total snowfall registered by the ultra-sonic ranger was 50 mm w.e.
We calibrated the model parameters to obtain the least root mean square (rms) error between modelled and measured b sfc (Fig. 5), with the criteria that σ ice >σ snow to reflect the higher absorption of short-wave radiation on an ice surface.As pointed out by others (Hock, 1999;Schuler et al., 2007), the calibration of the current formulation of the model is not unambiguous as it is possible to find different sets of tuning parameters (Table 2) giving equally good fits to the observed data.In addition to having least rms, the ratio between δ snow or δ ice should be similar to the typical ratio of snow to ice albedo as reflecting the physical basis for using a dual δ, and give the average best fit (lowest rms) for the three validation periods (Fig. 5).However, this latter condition does not involve any tuning of the parameters in the validation runs, only a decision on acceptance or rejection.

Temperature and energy fluxes during the melt seasons
For simple comparison, we set the periods of SEB analysis to a standard melt season, 15 December-15 March, which corresponds to the period with most surface melt.Melt seasons are characterized by a mean temperature close to zero (Table 3).Hourly mean temperature rarely deviates more than ±5 • C from zero (Fig. 6), due to the maritime setting of the glaciers.Maximum (minimum) hourly mean temperature during the four seasons was 6.8For all seasons the main energy source (sink) was L in (L out ) and the resulting mean long-wave net flux was negative.The net short-wave flux, S n , together with L in , show the highest inter-seasonal variations, but are opposing each other yielding a reduced variation in mean seasonal A to a maximum of 5 W m −2 , which is equivalent to 110 mm w.e. of melt for a full melt season.The inter-seasonal variation in S n is driven, with similar weights, by the variations in α and S c .S n was for all seasons the individual source with best fit to A (Table 4) and was also, expressed as ρ, the flux having the highest influence on the variation in A. In 2009/10 ρ was 120 %, as the variation was counteracted by L in .The generally negative coupling between L in and S n is mainly an effect of their different response to cloud cover, as in general clouds absorb and scatter a great portion of the incoming short-wave radiation, preventing it to reach the ground, but emit long-wave radiation to the ground as a function of their temperature.ρ for L in varies between the individual seasons from positive, through insignificant, to negative as a result of the balancing of the energy fluxes as exemplified in Fig. 7 that shows the two seasons with most contrasting standard deviations, r and ρ for S n and L in .The lower r and ρ for L in in 2009/10 is mainly explained by the occurrence of more extreme low daily mean L in values that are not associated with low daily mean A, as they are being balanced by high daily mean S n .Furthermore, the days with low mean A are distributed towards higher daily mean L in during 2009/10, while the days with high daily mean A are similarly distributed for both seasons.
A result of the negative correlation between the major energy fluxes L in and S n , particularly marked during season 2009/10, is the strong coupling between net radiation, R, and   A. R explains up to 82 % of the variation of A for an individual season.The turbulent fluxes H and E, have, considering their low mean absolute values, a proportionally large impact on the variation of A that was largest in 2008/09, when they together explain a third of the variance of A while the influence of S n was reduced.The coupling between T air and the energy fluxes constituting A is the physical basis of the temperature-index models.T air was best correlated with A during the seasons 2006/07 and 2008/09, coinciding with the highest ρ of L in .The generally low correspondence between A and T air , in particular for 2007/08 and 2009/10, is discouraging for using a solely temperature-driven ablation model, both at point scale and distributed over the glacier.T air and S n were slightly negatively or insignificantly correlated, with the highest correlation factor (-0.47) occurring in 2007/08.
The degree of explanation for the variation of A when using a combination of T air and S n becomes similar for all seasons and is just slightly lower than the corresponding degree of explanation for R. In a combination where S n is substituted by D, thus corresponding to the input parameters in the temperature radiation-index model (Fig. 8), the degree of explanation is only slightly lower.Using a distributed model driven by D instead of S n reduces the complexity and avoids the source of error associated with introducing a spatial and temporal modelling of α that is required for a model driven by S n .
The high peaks in melt (Fig. 8) coincide with high turbulent fluxes, driven by high relative humidity and positive T air .Backward trajectory analyses using the NOAA Hysplit trajectory model (http://ready.arl.noaa.gov/HYSPLITtraj.php)show that these events are associated with air masses rapidly moving in from the NW.The indication of higher melt from the ultra-sonic ranger is possibly at least partly an effect of the compaction of the snow during the wet conditions, and partly an effect of the SEB and the temperature-index models not accounting for the additional melt effect of rain and fog.

Model validation and sensitivity analysis
The calibrated model was run on three periods from other melt seasons for which corresponding stake measurements and data from AWS were available (Table 5).The period from season 2006/07 is significantly shorter than the others because the installation of the AWS took place several weeks after the first stake reading of the season and thus the period of correspondence starts with the second set of stake readings of the season.
The modelled B sfc generally agrees well with B sfc obtained by the glaciological method (Table 5).The difference is largest for season 2007/08, when the model underestimates B sfc by 140 mm w.e., which is slightly larger than the typical error of the glaciological method (ca. 100 mm w.e.; Jansson, 1999).The scatter from a linear relation between observed and modelled b sfc is large for the short validation period 2006/07.This might be explained by the fact that during the validation period a large part of the surface at the position of the stakes changed from snow to ice conditions, and thus the actual more gradual transition between snow and ice compared to the model's binary set-up may reduce the model performance.B sfc in 2009/10 was small because of a combination of little melt and extensive snow fall.The discrepancies between the modelled and the observed distributions of the snow fall possibly introduces a significantly larger error compared to the other years, which is exemplified by an outlier that was falsely modelled as ice surface for much of the period.The difference in seasonal mean T air and D is modest.A higher inter-seasonal variation, reflecting more contrasting weather conditions between the periods, would possibly have enhanced the validity test of the model.
We analysed the sensitivity of the model by applying perturbations to the calibrated model parameters, and to the meteorological parameters by manipulating the original AWS data corresponding to the calibration period.
Changes in m by ± 10 % produce a change in melt of a few percent (Fig. 9), while a similar change in δ snow/ice varies the total modelled melt by about ± 10 % (or ± 50 mm w.e.).A change by ± 50 % in snow cover thickness over the Snow cover depth, ± 50% q i , ± 0.1 T air , ± 0. entire glacier surface, which influences the temporal and spatial distribution of snow and ice, represented in the model through variations in the use of δ snow or δ ice , changes the melt rates by less than 10 %.
A change in dT/dZ is equivalent, for the model results, to a change in ablation gradient.A change of dT/dZ by ± 0.2 K km −1 in the calibration run changed B sfc by ± 5 %.As the average T air is close to zero and the AWS is located close to the ELA, the area-integrated change in melt due to a more negative dT/dZ will to some extent balance out the larger changes in individual grid cells, but will result in changed rms of the difference between measured and modelled melt for the stakes.A 0.2 K km −1 increase of dT/dZ raised the rms of the calibration run to 160 mm w.e., while a corresponding decrease in fact lowered the rms slightly.
We run the model with a step increase of 0.5 • C in T air , which is of the same magnitude as one standard deviation of the summer mean temperature for the last 30 years from the nearest long-term temperature instrument record (Bellinghausen Station on KGI), and close to half a standard variation of the daily mean T air for the calibration period.The corresponding modelled melt rates increased by 56 % to 805 mm w.e.A temperature decrease of the same magnitude resulted in a reduction of melt by 44 % to 296 mm w.e.
A similar sensitivity was obtained when perturbating analogously the threshold temperature for onset of melt (Eq.10).There is a strong temperature threshold effect which is illustrated by the change in number of grid cells integrated over all the time steps where the model indicates melting conditions.For the temperature increase (decrease) scenario the change was +52 % (-54 %).
Mean T air for the validation period 2009/10 was 0.5 • C colder than that of the calibration season 2008/09, while D and S n registered at the AWS were similar.The periods are of similar lengths and we note that the measured melt rate in 2009/10 was 38 % lower than the melt measured in 2008/09, which is close to the modelled response to a 0.5 • C temperature decrease on the 2008/09 data.
Because of the limited elevation range of the glacier, the entire glacier area experiences temperatures flickering around zero during summer time.This leads to a high sensitivity to temperature changes during the melt seasons.
We perturbed the radiation climate represented in the model by changing θ i , as the other inputs to obtain D are not related to the meteorological conditions on the glacier.We run the model with a step change in θ i that, just as in the T air scenario, corresponds to half a standard deviation of the daily means (0.06).An increase in θ i (implying more direct solar radiation) leads to 15 % more melt (80 mm w.e.), and a similar decrease in melt for a corresponding decrease in θ i .
T air and θ i are anti-correlated in the AWS record, meaning that under present climate a temperature increase is in general associated with a cloudier sky.Extrapolating this to a temperature increase driven by climate change implies that an increased melt due to higher temperature will to some extent be balanced by reduced direct radiation.However, climate change at these latitudes will probably be mainly driven by the associated changes in cyclonic activity and pathways, hence changes will neither be static nor solely impact a single meteorological parameter, and the model sensitivity as a predictive tool provides only a first level of understanding.
A high sensitivity of the mass balance to changes in T air for glaciers on the nearby KGI was pointed out by Knap et al. (1996) from the results of running a simple ice-flow model forced by an energy-balance model, after perturbating T air , until a new equilibrium state was reached.An energybalance model based on single-point measurements on Ecology glacier, KGI, produced an increase in ablation by 15 % in response to a 1 K temperature increase (Bintanja, 1995).This considerably lower sensitivity as compared to our results can be explained by almost constantly positive temperatures over the 30-day period of Bintanja's study.Thus, the effect of the 0 • C threshold will be considerably lower.A contrasting situation is probably behind the lower sensitivity -as compared to ours -(27 % increase in ablation as a response to a 1 K air temperature increase) obtained by Braun and Hock (2004) when applying a distributed energy balance model to the western part of KGI ice cap.In this case, the hypsometry of the ice cap indicates that a great part of the area was for most of the six week study period well below 0 • C.An additional explanation of the higher sensitivity found in our study is that the two KGI studies mentioned above were performed during limited time periods in December and early January, which is before the usual period of strongest seasonal melt.

Conclusions
We have used an AWS record located on Hurd Peninsula glaciers, Livingston Island, to analyse the SEB for four melt seasons (2006/07-2009/10).Further, we set up and run a www.the-cryosphere.net/6/539/2012/The Cryosphere, 6, 539-552, 2012 temperature-radiation index model.The model was calibrated and validated with b sfc obtained with the glaciological method.A novel correction method of α i , which adds the diurnal variation of clouds on a 24-hour running mean α, is presented and adopted.The advantage of the method is that the rapid and large variations during the high fluxes at midday in the short-wave balance are better reproduced.The following main conclusions can be drawn from our analysis: 1. S n is the most important individual energy flux impacting on the variation in A, but using the combined expression of the radiation fluxes, R, the degree of explanation increases and can account for 76-82 % of the variation in A. This high degree of explanation is due to a balancing of the generally anti-correlated fluxes S n and L in ,which is mainly an effect of their different response to cloud cover.
2. The seasonal means of T air at the AWS site are close to zero during the four melt seasons, and due to the small altitude range of the glaciers their whole area shows temperatures flickering around zero degrees during the summer season.The poor correlation factors between daily mean T air and A questions the performance of a solely temperature-based melt model, both at point scale and distributed over the glacier, and supports the use of a temperature radiation-index model.A combination of T air with S n or D increases the correlation factor at the point scale.
3. The modelled surface mass balance is in good agreement with that obtained by the glaciological method.Differences between model results and observations were generally below the typical error of the glaciological method (ca. 100 mm w.e., Jansson, 1999), with no significant bias.
4. The model results show that these glaciers are very sensitive to air temperature changes.An increase (decrease) in temperature of 0.5 • C implies an increase (decrease) of the melt rates by about 56 % (44 %), which is an effect of the strong zero degree threshold for onset of melt.The high model sensitivity of these glaciers to temperature change are indicative, but it must be noted that it only provides a first level of understanding of the response of the glacier mass balance to real climatic changes.
5. An increase (decrease) in the fraction of potential top of atmosphere radiation that reaches the ground, θ i , by half a standard deviation of its daily mean leads to an increase (decrease) of melt by 15 %.

Fig. 1 .
Fig. 1.Location of Hurd Peninsula, Livingston Island, with the position of Juan Carlos I Station (JCI), the AWS and the mass balance stakes in 2008/09.
gusts with wind speeds above 30 m s −1 are frequent.The highest wind speeds come from East and North-East directions and are associated with deep westward facing low pressure systems passing north of Livingston Island.The annual average temperature at JCI during the period 2000-2010 was -1.1 • C, with average summer (DJF) and winter (JJA) temperatures of 2.8 • C and -4.4 • C, respectively, with maximum/minimum registered temperatures of 15.5 • C and -22.1 • C.

Fig. 3 .
Fig. 3. Change in θ i versus change in α i between subsequent AWS recording time steps (10 min) for the melt season 2007/08.The linear fit has slope -0.10.

Fig. 4 .
Fig. 4. Distance to snow surface registered by the ultra-sonic ranger.

Fig. 5 .
Fig. 5. Measured and modelled point b sfc of the individual ablation stakes for the calibration and validation runs (mm w.e.).Blue dots indicate that the model handled the surface of the position of the stake as ice at the end of the modelled period.

Fig. 7 .
Fig. 7. Daily mean offset to the seasonal mean L in versus corresponding S n for melt seasons 2008/09 and 2009/10.Black dots indicate days when A was larger than 150 % of seasonal mean A, grey dots when lower than 50 % and white dots when A was between 50 % and 150 %.

Fig. 8 .
Fig. 8. Surface melt at the position of the AWS according to the surface lowering registered by the ultra-sonic ranger (SR50), SEB calculations and temperature-radiation index model for part of the model calibration period.

Fig. 9 .
Fig. 9. Change in melt rate for the Hurd Peninsula glaciers resulting from perturbation in various parameters relative to the calibrated run of the 2008/09 AWS data.

Table 1 .
Instruments of the AWS.

Table 2 .
Temperature-radiation index model parameters.

Table 3 .
Mean daily T air ( • C), energy fluxes (W m −2 ), α c and θ i with their standard deviations, for the four melt seasons.

Table 4 .
Correlation coefficient (r), and corresponding ρ to A, for daily mean energy fluxes, T air and combinations of T air , S n , and D.

Table 5 .
Characteristics of model calibration and validation runs.Figures in brackets in column No. of stakes refer to the number of stakes observed in bare ice at the end of the period.B sfc values are negative (representing ablation), but in the table we display -B sfc for easy comparison with melt, which is defined as a positive quantity.
List of symbolsA atmospheric energy fluxes = R + H + E b sfc point surface mass balance for the specific time period considered for each season B sfc glacier-wide surface mass balance for the specific time period considered for each season d and d m instant and mean Sun to Earth distance The Cryosphere, 6, 539-552, 2012 www.the-cryosphere.net/6/539/2012/