Albedo over rough snow and ice surfaces

Both satellite and ground-based broadband albedo measurements over rough and complex terrain show several limitations concerning feasibility and representativeness. To assess these limitations and understand the effect of surface roughness on albedo, firstly, an intrasurface radiative transfer (ISRT) model is combined with albedo measurements over different penitente surfaces on Glaciar Tapado in the semiarid Andes of northern Chile. Results of the ISRT model show effective albedo reductions over the penitentes up to 0.4 when comparing the rough surface albedo relative to the albedo of the flat surface. The magnitude of these reductions primarily depends on the opening angles of the penitentes, but the shape of the penitentes and spatial variability of the material albedo also play a major role. Secondly, the ISRT model is used to reveal the effect of using albedo measurements at a specific location (i.e., apparent albedo) to infer the true albedo of a penitente field (i.e., effective albedo). This effect is especially strong for narrow penitentes, resulting in sampling biases of up to ±0.05. The sampling biases are more pronounced when the sensor is low above the surface, but remain relatively constant throughout the day. Consequently, it is important to use a large number of samples at various places and/or to locate the sensor sufficiently high in order to avoid this sampling bias of surface albedo over rough surfaces. Thirdly, the temporal evolution of broadband albedo over a penitente-covered surface is analyzed to place the experiments and their uncertainty into a longer temporal context. Time series of albedo measurements at an automated weather station over two ablation seasons reveal that albedo decreases early in the ablation season. These decreases stabilize from February onwards with variations being caused by fresh snowfall events. The 2009/2010 and 2011/2012 seasons differ notably, where the latter shows lower albedo values caused by larger penitentes. Finally, a comparison of the ground-based albedo observations with Landsat and MODIS (Moderate Resolution Imaging Spectroradiometer)-derived albedo showed that both satellite albedo products capture the albedo evolution with root mean square errors of 0.08 and 0.15, respectively, but also illustrate their shortcomings related to temporal resolution and spatial heterogeneity over mall mountain glaciers.

The representativeness of surface albedo with respect to the footprint of the sensor is an important parameter to take into account when using albedo derived from field measurements or remote sensing data.One decisive factor for the representativeness is the "macroscopic" roughness of the surface, which has a strong effect on the surface albedo (Warren et al., 1998;Zhuravleva and Kokhanovsky, 2011).Although this effect has been extensively quantified over sastrugis and crevasses (Hudson and Warren, 2007;Kuchiki et al., 2011;Leroux and Fily, 1998;Pfeffer and Bretherton, 1987;Warren et al., 1998), its effect over penitentes (spike formations of snow and ice up to several meters high; Lliboutry, 1953), typical of several high altitude mountain glaciers and snow fields, remains less understood and limited to individual measurements.For example, Corripio and Purves (2005) and Kotlyakov and Lebedeva (1974) noted albedo reductions of 8-10 % over penitentes.Penitentes, however, have a surface roughness that is often much larger than sastrugis (Kuchiki et al., 2011) and which evolves over the ablation season (Cathles et al., 2014).Therefore, variable effects on the surface albedo can be expected and quantifying these is essential to model and understand the energy balance of glaciers with penitentes.Warren et al. (1998) reviewed the effect of surface roughness on albedo and mentions two causes for albedo reduction over a sastrugi field.Firstly, sastrugis lower the averaged incidence angle, which reduces the albedo due to the strong dependence of albedo on the incidence angle of incoming radiation (Warren, 1982).This effect depends on the sun's azimuth position relative to the sastrugi axis, as perpendicular insolation results in an albedo decrease between 2 and 4 % relative to parallel insolation (Carroll and Fitch, 1981;Kuhn, 1974).Secondly, multiple reflections between the walls cause light trapping in the trough.In this framework Pfeffer and Bretherton (1987) define the effective albedo that differs from the flat surface albedo (i.e., albedo of a flat surface with identical surface material properties and illumination conditions) due to light trapping within crevasses, whereas the flat surface albedo differs from the material albedo (i.e., albedo when the incident radiation has a incidence angle of 0 • ) due to changes in the zenith angle of incoming radiation.The concept of effective albedo is useful as it combines both the surface properties of the material and the light trapping due to multiple reflections.Pfeffer and Bretherton (1987) developed a radiative transfer model to simulate that, depending on the opening angle, crevasses reduce the effective albedo up to 0.4.Furthermore, they show that the opening angle of the crevasse determines the differences between effective albedo and flat surface albedo as smaller opening angles (i.e., steeper walls) result in stronger albedo reductions.Cathles et al. (2011) extended the radiative transfer model to differently shaped channels and crevasses and found a decrease in effective albedo over time due to changing morphologies of the roughness features.Alternatively, Fortuniak (2007) presents a radiation model to simulate effective albedo over urban canyons as a function of height-to-width ratios, whereas Helbig et al. (2009) did use a radiosity approach to estimate effective albedo over complex terrain.
Although the use of radiative transfer models (Cathles et al., 2011(Cathles et al., , 2014;;Fortuniak, 2007;Helbig et al., 2009;Pf-effer and Bretherton, 1987) allows quantifying the effect of surface topography on effective albedo, their use in energy balance models remains limited (e.g., Corripio and Purves, 2005) as the exact rough topography often remains unknown.Instead albedo measurements derived from hemispherical shortwave radiation sensors or remote sensing data are typically used as effective albedos in the energy balance models (Corripio and Purves, 2005;Pellicciotti et al., 2008;Winkler et al., 2009).However, the albedo measured over a rough surface may be quite different from the effective albedo depending on the position and footprint of the sensor, as penitente surfaces are heterogeneous in their incoming/outgoing radiation (Corripio and Purves, 2005).In this context, Pirazzini (2004) discusses the apparent albedo (i.e., the albedo measured under particular geometric conditions) and how it can differ from the "true" or effective albedo depending on the position of the sun/sensor with respect to the surface, and the shape, size, and orientation of the surface topography.This stresses the need for a comprehensive understanding of the differences between flat surface albedo, apparent albedo and effective albedo over a rough surface.This understanding is especially important when using albedo data for validation of remote sensing imagery, interpretation of automated weather station (AWS) radiation data or incorporation in energy balance models.
This paper aims to address the current need for a more thorough understanding of the effects of penitentes on surface albedo and how it can vary depending on the position of the sensor and size/shape of the penitentes.More specifically, the objectives are (i) to assess the effect of penitente size and shape on the outgoing radiation and effective albedo; (ii) to quantify the difference between flat surface albedo, apparent albedo and effective albedo measured by a sensor placed at different heights above a penitente surface; and (iii) to use the uncertainty related to the use of apparent albedo data for comparing albedo data from AWS measurements to satellite observations.Within this framework, an intrasurface radiative transfer model (ISRT) is used to simulate the incoming/outgoing radiation within a penitente trough and the apparent and effective albedo above a penitente surface.The simulated radiation and effective albedo data derived from the radiative transfer model are subsequently compared to radiation and apparent albedo measurements over a real penitente surface with varying geometrical/sun conditions.Moreover, the uncertainty due to apparent albedo is put into context by presenting albedo time series for two markedly differing ablation seasons and comparing them with satellitederived albedo values.

Study area
This study was performed on Glaciar Tapado (30 • 08 S, 69 • 56 W, Fig. 1), the largest glacier of the upper Elqui River catchment, close to the border between Chile and Argentina.The glacier is situated in the semi-arid Andes, south of the Arid Diagonal, its elevation range is between 4600 and 5536 m (Cerro Tapado) and its size 1.05 km 2 .The climate is characterized by predominantly clear skies, intense solar radiation, low air humidity and low precipitations.Higher peaks adjacent to Cerro Tapado, such as Cerro Olivares (30 • 17 S, 69 • 54 W, 6252 m), are currently free of glaciers, suggesting that the few glaciers existing in the area are atypical features and that local climatic conditions (e.g., excess precipitation due to wind redistribution of snow) play an important role (Gascoin et al., 2011(Gascoin et al., , 2013;;Ginot et al., 2006;Kull et al., 2002).

Penitente surface topography
Four different penitente surface topographies were sampled during individual experiments over the 2012/2013 ablation season (Fig. 2).For each of the experiments the penitente geometry (penitente height H and width W over one trough) and sun geometry (solar zenith angle θ, solar azimuth angle  φ) was assessed (Fig. 3; Table 1).Each experiment showed elongated penitentes with an east-west orientation of the ridges and troughs but with little or no tilt.

Radiation and albedo measurements vs. sensor height
To quantify the variation in outgoing radiation and apparent albedo due to the changes in the position of the sensor, a tripod made of 6 m long aluminum stakes was installed on the glacier surface over a penitente trough (Fig. 2) during the four experiments.A downward looking pyranometer (see Table 2 for details) was mounted to a weight hanging on a cord from the tripod top.During each experiment the outgoing radiation (S app out ) was recorded at different heights (h) above the penitente tips in 0.5 m steps and at different depths (d) within the troughs in 0.25 m steps (Fig. 3).Simultaneously the incoming radiation (S in ) was measured on a second, fixed tripod ca. 1 m above the surface in an open area where penitentes were small enough to have no effect on the upward-looking sensor.The distance between the second tripod and the downward-looking sensor was set large enough to avoid mutual influence.All experiments were performed under cloud-free conditions with no or negligible wind to avoid movement of the downward looking pyranometer.

ISRT model
A two-dimensional (2-D) ISRT model similar to the models of Cathles et al. (2011Cathles et al. ( , 2014) ) and Pfeffer and Bretherton (1987) was run additionally to the radiation and albedo measurements, as it allows quantifying the effect of penitente surface topography on effective albedo and assessing the differences between flat surface albedo, apparent albedo and effective albedo.The ISRT model consists of (i) representing the penitentes and sensor as 2-D geometric shapes composed of small segments (Sect.3.3.1)and (ii) numerically solving the shortwave radiation received/reflected by each segment (Sect.3.3.2).

2-D penitente field
Penitentes are well suited to be described in 2-D given their elongated shape in the east-west direction (Cathles et al., 2014).Therefore, the penitente field of each experiment was represented by simulating a statistical population of 2-D representative penitente surfaces (i.e., 75 samples per experiment).This was done based on the measured size parameters and for different shapes.Firstly, representative samples of penitentes with similar heights and widths were created over a 40 m transect, corresponding to the diameter of the sensor footprint when the sensor is at h = 4 m (i.e., 99 % of the signal is coming from a viewing angle of ±84 • ).For these samples the size of the central penitente trough below the sensor was defined based on the measured size parameters H and W .The size of the neighboring penitentes in the north-south direction was defined as H = H ± h and W = W ± w (Fig. 3) based on the assumption that the neighboring penitentes have randomly varying geometry (i.e., higher/lower/wider/narrower).h and w are random samples from a normal distribution N (µ = 0, σ = 7.5 cm) due to the lack of measurements for the neighboring penitentes.Secondly, different penitente shapes were simulated for each representative sample as the actual shape of the penitentes was not assessed during the experiments.These shapes range from triangular to convex-, concave-and cosine-shaped penitentes (Fig. 4) that were modeled according to Table 3.None of the modeled shapes corresponds to real observed penitente geometries, which often show more complex shapes that evolve over time (e.g., Betterton, 2001;Cathles et al., 2014), but using this range of shapes allows understanding the variability in radiation and albedo due to varying shapes.Moreover, this range of shapes covers most variability described over penitentes.

Radiative transfer
The 2-D representative penitente surfaces were subsequently used in an intrasurface radiative transfer analysis to calculate the shortwave radiation received/reflected along the surface or measured by the sensor.This was done by dividing the surface in 2.5 cm segments and calculating the view factors for each segment s.The view factors account for viewing obstructions and multiple reflections between segments   3. The blue/green shaded areas for the convex/concave shapes represent the variability in shapes due to the uniform distributions U (min, max) in Table 3.The dotted/dashed grey lines indicate the position of the penitente trough/tip, whereas the arrows illustrate the angle of direct incoming radiation during experiments A-D.Note the different y scales.
as they quantify the proportion of radiation coming from another segment s using the assumption of Lambertian surfaces that follow the cosine law of illumination.Based on the view factors the amount of incoming/outgoing radiation for each segment can be calculated by solving where S s in is the amount of incoming radiation on segment s, S s out is the amount of outgoing radiation on segment s, I d is the component of direct incoming radiation from the sun, I i is the component of indirect diffuse incoming radiation from the sky, α s →s is the albedo of segment s for the radiation coming from segment s , and F s →s is the view factor for radiation going from s to s.The calculations of the view factors were performed based on the adaptive integration approach of Walton (2002), which allows calculating for all segments s the amount of radiation coming from segments s based on the distance between segments, possible viewing obstructions, and the angle between the segment normals.
Equation (2) contains albedo terms α d→s , α i , and α s →s that are dependent on the source of the incoming radiation (d: direct sunlight; i: diffuse radiation; s : radiation from segment s ) and that account for the albedo dependence on the incidence angle of incoming radiation (Warren, 1982), where the incidence angle is the angle between the radiation rays and the normal of the surface segment.To include this dependence in the ISRT model, the parametrization of Gardner Table 3. Equations used for simulating the different penitente shape geometries based on the size parameters H and W (note: for the penitente trough below the sensor H = H and W = W ), where x is the horizontal coordinate centered around each penitente trough, y is the vertical coordinate and z is a sample of the uniform distribution U (min, max) for the convex and concave shapes.These equations relate to the penitente shapes shown in Fig. 4.

North face
South face Parameter and Sharp (2010) for broadband albedo was used to calculate α d→s , α i , and α s →s with where u is the incidence angle of incoming radiation from segment s on s, and α s mat is the material albedo of segment s (i.e., albedo for u = 0 ).Equation ( 3) is also applied for α i and α d by adopting an effective incidence angle of u = 50 • for pure diffuse radiation (Wiscombe and Warren, 1980), or by using the solar incidence angle u on segment s for the direct sunlight albedo.
One of the advantages of the radiative transfer calculations is that they allow determining the amount of outgoing radiation that effectively leaves the penitente troughs and/or reaches the sensor at height h or depth d: where S eff out is the upward radiation flux leaving the penitente trough through cross section W (Fig. 3) and S app out is the upward radiation flux from the surface that reaches the sensor that also follows the cosine law of illumination.
To solve Eqs. ( 1)-( 3) we calculated the view factors, incidence angles of incoming radiation, and shading for each segment (i.e., if shaded I d = 0).Moreover, we determined the fraction of direct and diffuse radiation, I d and I i , respectively (Table 1), based on the SBDART (Santa Barbara DISORT Atmospheric Radiative Transfer) model (Ricchiazzi et al., 1998) for a tropical, dry, clear sky atmosphere without aerosols at 5500 m elevation, where SBDART takes into account the radiation angles and the location to calculate I d and I i .As a result, the only unknown in Eqs.(1)-( 3) is the material albedo α s mat .However, by making assumptions on the spatial variability of the material albedo and using the apparent albedo measurements at h = 0 m, we can solve Equations ( 1)-( 3) and ( 5) and derive α s mat .For example, S s in , S s out and α s mat can be determined for each segment by assuming a uniform material albedo or a spatial albedo gradient within a penitente trough (e.g., the material albedo decreases when going from tip, typically made of firn or icy firn, to bottom, typically made of ice with standing or running water in mid-summer).Since field observations of the spatial variability in α s mat were lacking, both the uniform and gradient assumptions were tested as examples of the effect of spatial variability in α s .The assumption of a spatially uniform α s mat is perhaps not very realistic, but it allows a straightforward interpretation of the difference between the material albedo α mat , flat surface albedo α flat and effective albedo α eff (e.g., Warren et al., 1998).The albedo gradient, however, represents more realistic conditions since meltwater, debris, and dust tend to accumulate in the trough bottom (Cathles et al., 2014).

Material albedo, apparent albedo, and effective albedo
Finally, the ISRT model allows quantifying the difference between flat surface, effective and apparent albedo by comparing α flat and α eff with α app at different sensor heights/depths, as defined below.For this purpose, the flat surface albedo can be derived by introducing α s mat in Eq. ( 3) for the local solar conditions, whereas The flat surface albedo should be interpreted as the albedo of the surface material for a flat surface without roughness, the apparent albedo as the albedo measured by a sensor under particular geometric conditions, and the effective albedo as the ratio of radiation leaving/entering a penitente trough.Consequently the difference between flat surface albedo and effective albedo could be considered the macroscopic effect of the surface roughness on the surface albedo.The local solar conditions for each experiment (i.e., solar irradiance, sun geometrical conditions) were used to compare the α flat , α app , and α eff with the measured albedo and The Cryosphere, 8, 1069-1086, 2014 the effect of the position of the sensor and size/shape of the penitentes on them.Moreover, the diurnal variation in α app and α eff was calculated at h = 2 m (i.e., a typical height of an AWS radiation sensor) to verify if the differences between effective and apparent albedo vary throughout the day, which is important for the interpretation of diurnal albedo data (e.g., diurnal mean albedo).

Temporal albedo evolution
To analyze the evolution of broadband albedo over a penitent-covered surface an AWS was operated in the ablation zone of Glaciar Tapado (Fig. 1) during two ablation seasons (2009-2010; 2011-2012).Short-wave radiation was measured every 10 s and averaged to hourly values.Sensor specifications are summarized in Table 2.The radiation sensor was installed horizontally over a 10 • inclined surface with an aspect of 134 • (values derived from a 2010 GeoEye DEM, resampled to 20 m cell-size).The same correction as in Abermann et al. (2014) following Grenfell et al. (1994) has been applied to adjust albedo for the slightly sloping surface; however, deviations from the uncorrected data are small (mean difference between uncorrected and corrected 0.02) as radiation values are close to solar noon and the aspect of the surface is not very different from the north-south axis.Additionally, a 95 % confidence interval on the AWS albedo data was calculated by combining the measurement error (Table 2) with the error due to the use of the apparent albedo instead of the effective albedo when using AWS measurements: where ε AWS , ε sensor , and ε app−eff are respectively the standard errors for the AWS measurements, the sensor, and the use of apparent instead of effective albedo and where the 95 % confidence interval is ±2 times the ε AWS .
The apparent AWS albedo measurements were subsequently compared to satellite-derived albedo from MODIS (Moderate Resolution Imaging Spectroradiometer) and Landsat sensors.MODIS albedo (α MOD ) was derived from the MODIS daily snow product (MOD10A1) processed by the National Snow and Ice Data Center (Hall et al., 2007;Stroeve et al., 2006).MOD10A1 provides daily estimates of the snow cover, snow fraction and albedo of snow surfaces at 500 m spatial resolution based on the Terra overpass (i.e., around 10:30 Local Solar Time (LST) at the Equator) with viewing angle closest to nadir.Additionally, Landsat albedo (α L7 ) was derived from Landsat ETM+ (Enhanced Thematic Mapper Plus) images based on the mean and standard deviation of the nine pixels surrounding the AWS location.A total of 13 Landsat ETM+ images above the Tapado region were selected over both ablation seasons based on data availability of cloud-free Landsat scenes (both visually and from metadata).Landsat ETM+ data provides spectral radiance in seven bands at the top of the atmosphere (30 m spatial resolution, 16 day temporal resolution, acquisition around 10:00 LST at the Equator), which was converted to Landsat spectral reflectance at the surface and broadband albedo using the methodology of Klok et al. (2003) based on the anisotropic reflection factor of snow/ice in combination with the parameterization of Knap et al. (1999).

Radiation within the penitente trough
Figure 5 shows the effect of penitente geometry and sun geometry on the variation in outgoing radiation along a penitente surface and illustrates the importance of surface orientation, multiple reflections, and shading on the radiation distribution within a penitente trough.
Surface orientation plays an important role as segments orientated perpendicular to the incoming radiation will receive/reflect more incoming radiation than tilted surfaces (Vanonckelen et al., 2013;Veraverbeke et al., 2010).For example, it is clear that north-facing slopes that are oriented towards the solar incoming radiation, receive/reflect more radiation than south-facing slopes.The shape of the penitentes also has a strong influence on the radiation distribution within a penitente trough since it will determine the exact orientation of each segment s .This effect of shape can easily be perceived when looking at the spatial distribution of S out in experiment A, where the maximum radiation is received/reflected either at the penitente bottom for a convex shape, or at the tips for a concave shape depending on which segments s are more oriented towards the incoming radiation from the sun.
Multiple reflections have an important effect on the spatial distribution of incoming/outgoing radiation as they determine how much radiation is received from other segments.This effect can range up to 50 % of the incoming radiation for an α mat of 0.65 in Fig. 5, but will vary depending on α mat .Furthermore, for multiple reflections the radiation distribution is strongly influenced by the surface orientation as it determines how much radiation each segment "sees" from the other segments.
Shading also has a strong influence on the radiation distribution within a penitente trough as it produces I d = 0 in Eq. ( 1), resulting in a strong decrease in total incoming radiation along the penitente surface.Experiment D, for example, was performed with shading of the south-facing slope and part of the north-facing slope (Fig. 2), which has a dominant impact on the variability in S in and S out (Fig. 5).

Outgoing radiation within penitente trough
The spatial variability in incoming/outgoing radiation has a strong effect on the radiation measurements within the penitente trough as can be seen in Fig. 6 depths.It shows that the range of S out measurements by the sensor can be represented by the ISRT model for different geometries, except for experiment D where the decrease of S out with depth is overestimated due to the overestimation of the amount of shading within the trough.Additionally, Fig. 6 demonstrates that the location of the sensor has a strong effect on the measurements of outgoing radiation as it will determine the contribution of each segment to the radiation received by the sensor.For example, the segments at the bottom of the trough will contribute relatively more to the measured radiation when the sensor is positioned low within the trough.As a result, the S out measured by a sensor near the bottom is higher for convexly than for concavely shaped penitentes since both S s in and S s out of the bottom segments are higher (Fig. 5).

Apparent and effective albedo vs. sensor height
Table 4 shows the difference in material, flat surface, apparent and effective albedo for the penitente troughs for each of the experiments based on the ISRT model.Firstly, it illustrates that both size and shape of penitentes have a strong effect on the differences between material and flat surface albedo, on the one hand, and apparent and effective albedo, on the other hand.The size of the penitentes primarily determines the amount of light trapping in the trough as pen-itentes with a higher H /W ratio imply larger reductions in the α app and α eff in comparison with α flat .The albedo reduction can be up to 0.4 for the experiments in this study.The shape however also has an influence on the albedo reductions relative to a flat surface as concavely shaped penitentes result in lower albedo reductions than convexly shaped penitentes.This again can be explained by differences in segment orientation, as convex shapes are more open and consequently cause more light trapping.
Secondly, Table 4 illustrates the effect of using a sensor halfway between the tips (α app , when h = 0) to quantify the effective albedo of a penitente trough (α eff ).Although the differences between α app and α eff generally are small, they can reach errors of up to 0.07 depending on the radiation distribution within a penitente trough.These errors become particularly large when there is no uniform α mat within the trough.
Figure 7 demonstrates the effect of sensor height on albedo measurements by displaying the measured and modeled apparent albedo as a function of sensor height.The measurements show that during experiment A, when snow penitentes were ∼ 1 m high and the tips of the individual penitentes closely spaced, the apparent albedo rose from 0.32 at h = 0 m to 0.39 at h = 3.The solid lines represent the S out for the uniform material albedo whereas the blue/green shaded area is the variability in S out due to the uniform distributions U (min, max) in Table 3.
with sensor height.The large range of modeled changes in α app with height, however, indicates that the individual measurements during the experiments are difficult to interpret as the changes in α app with height are strongly dependent on the penitente shape and spatial distribution of material albedo, which was not assessed for the experiments.Nevertheless, the observed changes of α measured in Fig. 7 suggest that the penitente surface of experiment A corresponds more to cosine-shaped penitentes with an albedo gradient, whereas experiments B-D do not display this increase in α measured related to cosine-shaped penitente and/or an albedo gradient.Although data are lacking to verify this interpretation, this could be related to the flatter penitente tips and stronger albedo gradients over snow-ice surfaces in the early-summer experiment A compared to the pointy penitentes in B-D that consist only of ice later in summer (Fig. 2).Additionally, Fig. 7 reveals the difference between the albedo measured by a sensor at a specific height (i.e., apparent albedo) and the true albedo of a penitente field (i.e., effective albedo).It illustrates how the mean apparent albedo gradually converges to the effective albedo when the sensor is positioned higher above the penitente tips (i.e., sampling bias or α app − α eff < ±0.01).Below 2 m, however, systematic biases between the apparent and effective albedo occur.This is also clear in Fig. 8, which shows the sampling bias of the individual ISRT model runs and also contains the systematic underestimation/overestimation of α app for the lowest meter in A and D, respectively.This systematic underestimation/overestimation is the consequence of the location of Table 4. Overview of material albedo α mat , flat surface albedo α flat , apparent albedo α app at h = 0, effective albedo α eff and satellite albedo (α L7 and α MOD ) for each experiment and for the different shapes and assumptions on spatial variability in material albedo.The albedo gradients are based on the assumption that the penitente tips/bottoms have a material albedo which is 0.1 higher/lower than the mean material albedo and that there is a linear gradient between them.The satellite albedos also mention the closest acquisition date that corresponds to the four experiments.MODIS and Landsat acquisition times are around 10:30 and 10:00 LST at the Equator, respectively.the sensor in combination with the large viewing obstructions and heterogeneity in S out along the surface.This effect is illustrated in Fig. 9, which shows the spatial variability in S out along the penitente surface for experiment D in combination with the view factors of the sensor indicating how much each segment contributes to the radiation received by the sensor.
For h = 0 m, only the central penitente trough is seen by the sensor.When the sensor is positioned higher, the sensor receives less radiation from the central trough and more radiation from the neighboring penitente troughs.However, when the sensor is still low (e.g., h = 1 m), it specifically receives radiation from the segments of the neighboring troughs that have an above average S out .This explains why the apparent albedo for D at h = 1 m is overestimated, resulting in a sampling bias.Once the sensor is positioned higher (h > 2 m), this bias becomes less important, since the sensor measures more homogeneously over the neighboring valleys.Although Fig. 8 suggests that the mean difference between apparent and effective albedo is relatively small when the sensor is high enough (h > 2 m), the wide confidence interval of α app − α eff indicates that individual measurements of albedo still can contain sampling biases of up to ±0.05 (95% confidence interval for α app − α eff differences).This bias is specifically strong for narrow penitentes, where the effect of viewing obstructions is more pronounced.
Finally, Fig. 10 displays the effect of sun geometry on the sampling bias by presenting the modeled diurnal evolution of apparent albedo, effective albedo and the difference between both.It demonstrates that the difference in α app −α eff remains relatively constant during the day, implying that multiple observations throughout a day cannot balance the sampling bias for individual albedo measurements.Moreover, it shows that shading alone does not have a very strong influence on the effective albedo of the penitentes as no clear changes in α eff are observed when the penitentes are partly shaded or not.For example, in experiment A shading occurs before 11:00 LST and after 13:00 LST, whereas the effective albedo remains relatively constant between 09:00 and 16:00 LST.

Temporal albedo evolution
Figure 11 shows the temporal evolution of both satellite and ground-based broadband albedo for the 2009/2010 and 2011/2012 seasons together with the confidence intervals.The confidence interval is based on sampling biases of up to ±0.05 (95 % confidence interval), which corresponds to ε app−eff = 0.025, and combined with ε sensor = 0.029 in Eq. ( 8), results in an AWS albedo of ±0.08 (95 % confidence interval).Analysis of the AWS albedo time series reveals that both seasons differ significantly in their albedo evolution.The 2009/2010 season started with higher albedo values of more than 0.6 and then continuously decreased until mid-January, where a significant snowfall event raised the albedo to 0.8, after which metamorphosis and ice-exposure reduced it again to 0.2.In 2011/2012, values in early December were around 0.45.Throughout the season several fresh snow events temporarily raised albedo.From late January onwards the albedo showed values of about 0.25-0.30.The difference in albedo time series between the seasons is clearly attributable to a very distinct evolution of penitentes (Fig. 12).Whereas in 2009/2010 penitentes started to form in December and did not get higher than 1 m before February, in 2011/2012 there were already small penitentes during the first visit in late November and by February they were higher than 1.5 m.
Comparison of the ground-based albedo observations with Landsat-and MODIS-derived albedo show that both satellite-derived albedo products capture the albedo evolution including snowfall events (mean difference of −0.02 and root mean square difference of 0.08 (Landsat: 13 images) and mean difference of −0.12 and root mean square difference of 0.15 (MODIS: 35 images), respectively).Moreover Fig. 11 shows that for 12 of the 13 Landsat images the 95 % q q q q A q q q q α eff triangle α eff convex α eff concave α eff cos q q q q D Figure 7.The measured and modeled (i.e., mean of 75 samples) albedo changes in function of height above the penitente surface.The dots/squares represent the effective albedo α eff of the penitente surface that is not height dependent.The lines illustrate the apparent albedo α app in function of the sensor height, where solid lines are used for a uniform material albedo and the dashed lines for a material albedo gradient within the trough.
confidence intervals overlap.Nevertheless, the comparison of AWS and satellite albedo also illustrates the shortcomings of the satellite imagery with clear biases for the MODIS albedo (e.g., underestimation in 2009/2010 with deviations of up to 0.22) that complicate the monitoring of albedo evolutions among seasons.

ISRT model
Radiative transfer models such as the ISRT model provide a useful instrument to assess the effect of surface topography on the surface albedo (Cathles et al., 2011(Cathles et al., , 2014;;For-tuniak, 2007;Pfeffer and Bretherton, 1987).Although more detailed high-resolution measurements and wider footprints are necessary to actually reproduce the albedo measurements over penitente surfaces (e.g., the actual height dependence of α app in experiment A), the range of assumptions on shape and material albedo used in the ISRT model allows one to understand the variability due to shape and material albedo.Simultaneously, improved radiative transfer models that simulate the radiation over 3-D surfaces and that account for the anisotropy factor of light within snow/ice (Dumont et al., 2010) and its spectral variability will be required to fully understand the effect of roughness on surface albedo, because in the current setup snow/ice was assumed to be a Lambertian reflector for broadband shortwave radiation without taking the spectral variability in material albedo and incoming radiation into account.Although these assumptions will have an effect on the results, the 2-D representation in the ISRT model already provides a quantitive indicator of the uncertainties, because the ISRT model provides results that closely correspond to the range of observations (e.g., height dependence in Sect.4.3 or measured diurnal variability in Abermann et al. 2014).Nevertheless, we strongly encourage repeating analogue experiments using more complex radiative transfer models that account for spectral dependence and more detailed high-resolution measurements to improve our understanding of the effect of penitentes on the surface albedo and the energy balance of glaciers.

Effective albedo of a penitent
One of the advantages of the ISRT model is that it allows quantifying the difference between flat surface albedo and the effective albedo of penitente surfaces.For the penitentes in this study, for example, albedo reductions of up to 0.4 of a rough surface relative to a flat surface are modeled based on the ISRT model.The magnitude of these reductions primarily depends on the penitente opening angle (i.e., more open penitentes with a lower H /W have lower reductions).This effect is in accordance with the albedo reductions obtained for crevasses by Pfeffer and Bretherton (1987), who found very similar values when using similar opening angles, and for sastrugis by Warren et al. (1998), however, they limited their study to H /W < 1.The obtained albedo reductions differ from the values reported by Corripio and Purves (2005), who found smaller albedo reductions (i.e., 8 % for penitentes with H /W = 1.69) based on the comparison of two locations on one glacier (i.e., one with and one without penitentes) and the assumption that both locations have similar material albedos.The IRST model, however, proves that actually assessing the spatial variability of the material albedo in combination with the penitente shape is crucial to determine the actual albedo reductions over penitentes, as both have a large effect on the amount of light trapping within the trough.This importance of the penitente's shape was also shown by Cathles et al. (2011Cathles et al. ( , 2014)).
The ISRT model also indicates that shading within a penitente trough has no strong effect on the effective albedo as the α eff remains relatively constant over time periods where shading and no-shading occurs.The ISRT model shows only significant increases in effective albedo for high solar zenith angles, but this is in accordance with the diurnal changes in albedo observed by Abermann et al. (2014) under similar geographical conditions.
The large reductions in effective albedo over penitente surfaces stress the importance of understanding the variability in albedo due to surface roughness.Moreover, it is important to include surface roughness in the current parameterizations of surface albedo (e.g., Gardner and Sharp, 2010;Kuipers Munneke et al., 2011;van Angelen et al., 2012), certainly when they are used over surfaces with prominent macroscopic surface roughness.This is because over these surfaces the assumption that albedo reductions due to surface roughness occur only at high solar zenith angles, which have little impact on the overall surface energy balance (Gardner and Sharp, 2010;Warren et al., 1998), is no longer true.For example, this can be seen in experiment A where large albedo reductions also occur for small solar zenith angles with high incoming radiation.

Apparent albedo vs. effective albedo
The modeled differences between α app and α eff highlight the importance of understanding the effect of using albedo measurements under particular conditions to determine the effective albedo, because viewing obstructions and spatial variability in incoming and outgoing radiation may affect the radiance received by the sensor (Pirazzini, 2004).These effects are particularly strong when the sensor is low above the surface (e.g., h < 2 m), because viewing obstructions will then play a larger role.Nevertheless, although the viewing obstructions become less important when the sensor is higher above the surface, our results indicate that even at 4 m above the tips individual measurements can still show differences with the "true" or effective albedo of ±0.04.

Apparent albedo vs. remote sensing albedo
Despite the possible sampling bias in albedo measurements over rough surfaces, the correspondence between AWS albedo and MODIS-derived albedos (mean difference of  0.12) are similar to values obtained by Dumont et al. (2012) over mountain glaciers with lower roughness (mean difference of 0.11).On the other hand, less agreement between AWS and MODIS albedo is obtained than in Box et al. ( 2012) (mean underestimation of −0.02).However, they focused on monthly averages over more homogenous surfaces on the Greenland ice sheet and are therefore likely less influenced by surface heterogeneity.The Landsat albedo shows a better correspondence with the AWS albedo.These accuracies are in the order of magnitude of the Landsat albedo accuracies derived by Klok et al. (2003), who found an overestimation of 0.03, but with a lower root mean square difference (0.07).These biases between Landsat-derived and AWS albedo could be caused by the possible sampling biases due to albedo measurements from one location in combination with other possible explanations such as spatial heterogeneity (i.e., the experiment footprint is smaller than 9 m × 30 m) in dust and debris cover or the assumption of a flat surface when using the anisotropic reflection factor in the method of Klok et al. (2003).This flat surface is obviously not the case for penitente fields.This can also be seen in Table 4, where the albedos of the experiments are compared with the Landsat-derived albedo of the closest date (note: the different illumination conditions due to different date/hour result in changes in solar zenith angle of ∼ 10 • ) and where different biases occur when the penitentes further develop.The higher accuracy for the Landsat albedo compared to the MODIS albedo likely is related to the improved spatial resolution of the Landsat pixels (30 m vs. 500 m).In this context and given the spatial extent of the glacier, the MODIS-derived albedo is more related to the combination of glacier and surrounding-surface albedo.This spatial heterogeneity also explains the lower than daily temporal resolution of the MODIS albedo, although daily albedo values should be retrieved given the daily MODIS overpass.In reality, however, this daily temporal resolution is not obtained for the MOD10A1 product as the glacier is often misclassified as land/cloud.These misclassifications (i.e., only 15 % of the observations is classified as snow/ice, vs. 78 % land and 7 % cloud) show that the commonly reported snow vs. cloud misclassification (e.g., Dozier et al., 2008) does not play an important role in this area with low cloud cover (Gascoin et al., 2013), but that the snow vs. land misclassifications are more important.This could be explained by the subpixel fraction of the Glaciar Tapado within a MODIS pixel, causing the omission of the glacier's snow/ice cover in the α [ ] q q q q q q q q q q q q q q q q q Dec 2009 Jan 2010 Measured albedo MOD L7 α [ ] q q q q q q q q q q q q q q q q q q q q q q q Dec 2011 Jan 2012 MOD10A1 product.This misclassification effect due to surrounding land will moreover be enhanced as the pixel footprint of the MODIS sensor is often much larger than 500 m due to off-nadir viewing conditions (Dozier et al., 2008).As such, the advantage of the MOD10A1 product (i.e., daily temporal resolution) is relatively small compared to the spatial resolution of Landsat to retrieve effective albedo over small mountain glaciers with large spatial variability.

Implications for interpretation of albedo measurements
The large uncertainty in sampling bias due to use of albedo measurements over rough surfaces stresses the importance of a well-designed data collection setup.In this framework Pirazzini (2004) expresses the importance of minimizing the sampling bias over sastrugi fields by collecting a large number of samples at various places or over a long time span (Carroll and Fitch, 1981) or by increasing the height of the sensor above the surface, so that the irradiance measured will also be more representative of the effective reflectance of the surface (Warren et al., 1998).For penitentes, however, the ISRT models shows that the albedo uncertainty for individual measurements still is ±0.04 at 4 m above the surface.Moreover, we see that the uncertainty remains relatively constant throughout the day, indicating that daily means will contain similar sampling biases as individual measurements.Consequently, the only beneficial approaches to minimize the sampling bias of surface albedo over such rough surfaces is to increase the height of the sensor even more until the viewing obstructions are minimized or to use a large number of samples at various places.However, the practical difficulties to apply both approaches over a rough terrain should be taken into consideration.The benefit of the latter approach based on a large number of samples at various places is also confirmed by the IRST model as the mean of the representative samples, which could be considered the albedo of different locations, shows an uncertainty on the sampling bias below ±0.01.

Conclusions
The results of the ISRT model combined with the field data show that surface roughness due to the development of penitentes can produce large effective albedo reductions relative to the albedo of a flat surface.The magnitude of these reductions depends primarily on the penitente opening angle, but also the shape and spatial variability of the material albedo play a major role.Including surface roughness in the current parameterizations of surface albedo (e.g., Gardner and Sharp, 2010;van Angelen et al., 2012) is therefore essential if we want to correctly represent albedo in energy balance models over complex, rough terrain.
Up to now, the representation of penitentes in energy balance models has typically been done based on albedo measurements derived from shortwave radiation sensors or remote sensing data (e.g., Corripio and Purves, 2005;Pellicciotti et al., 2008;Winkler et al., 2009).The ISRT model, however, shows that individual albedo measurements or diurnal mean albedo data may contain sampling biases of up to ±0.05 depending on the location of the sensor and the roughness of the surface.Consequently, the only effective approach to reduce the sampling bias of surface albedo over rough penitente surfaces is to increase the height of the sensor even more until the viewing obstructions are minimized or to use a large number of samples at various places.This is important if we want to assess the accuracy of remote sensing measurements over heterogeneous glaciers.
The comparison of the ground-based effective albedo observations with Landsat-and MODIS-derived albedo showed that both satellite-derived albedo products capture the albedo evolution, but also illustrates the problems related to temporal resolution and spatial heterogeneity over heterogeneous glaciers and the difficulty to interpret this correspondence given the possible sampling bias.

Figure 1 .
Figure 1.Glaciar Tapado with the location of the AWS where albedo experiments were carried out.The inset shows the location in South America close to the Chilean-Argentinean border.

Figure 2 .
Figure 2. Illustration of the four experiments over different penitente surface topographies during the 2012/2013 ablation season.Details for each experiment are listed in Table1.

Figure 4 .
Figure 4.The different penitente shapes used in the ISRT model based on the equations in Table3.The blue/green shaded areas for the convex/concave shapes represent the variability in shapes due to the uniform distributions U (min, max) in Table3.The dotted/dashed grey lines indicate the position of the penitente trough/tip, whereas the arrows illustrate the angle of direct incoming radiation during experiments A-D.Note the different y scales.

Figure 5 .
Figure 5. Variation in outgoing radiation with and without multiple reflections (S out and S out , respectively) along a penitente surface for a constant material albedo (α mat = 0.65) for each experiment (A-D) and for the different shapes.The dotted/dashed grey lines indicate the position of the penitente trough/tip, whereas the arrows illustrate the angle of direct incoming radiation during experiments A-D.

Figure 6 .
Figure6.Measured and modeled outgoing radiation by the sensor placed at different depths within the penitente trough.The measured outgoing radiation has been normalized to have uniform incoming radiation.The solid lines represent the S out for the uniform material albedo whereas the blue/green shaded area is the variability in S out due to the uniform distributions U (min, max) in Table3.

αFigure 8 .
Figure8.Mean difference between apparent and effective albedo (lines) and the 95 % confidence interval (color shaded area) based on 75 samples per experiment for the different penitente shapes and for uniform material albedo.

Figure 9 .
Figure 9. Illustration of the modeled penitente geometry and the effect of height on the radiation received by the sensor: (a) simulated 2-D penitente surface of experiment D with a convex shape; (b) S out along the surface for a uniform material albedo; and (c) view factors (i.e., proportion of radiation measured by the sensor coming from segment s ) of the sensor along the surface for different sensor heights (h).

Fig. 10 .Figure 10 .
Fig. 10.Diurnal modeled evolution of the apparent albedo at h = 2 m, effective albedo and sampling bias (i.e., α app −α eff ) of two sample 2D penitente surfaces with convex and concave shape for the experiments A-D for solar zenith angles <75 • and for a constant α mat

Figure 11 .
Figure 11.Seasonal evolution of broadband albedo on Glaciar Tapado for the ablation seasons 2009/2010 (a) and 2011/2012 (b).The solid lines are arithmetic means of hourly values (slope corrected) between 10:00 and 11:00 LST for every day at the AWS, whereas the grey shaded area is the 95 % confidence interval for the AWS data.The blue dots represent the MODIS (MOD: 10:30 LST Terra overpass at the Equator) and the red squares the Landsat (10:00 LST overpass at the Equator) derived albedos with the 95 % confidence interval for the nine pixels surrounding the AWS.

Figure 12 .
Figure 12.Compilation of photographs from the AWS during the seasons 2009/2010 (left side) and the respective nearest available point in time in 2011/2012 (right side).Note the different penitente heights from since the beginning and their evolution throughout the season (Fotos: CEAZA glaciology group).

www.the-cryosphere.net/8/1069/2014/ The Cryosphere, 8, 1069-1086, 2014Table 1 .
Details of the individual experiments.θ is the solar zenith angle and φ the sun's azimuth angle during the experiment, H is the penitente height (tip to trough) of the trough below the sensor, W the distance between the nearest penitente tips.The ratio H /W is dimensionless.F d /F i are the fractions of direct (F d ) and diffuse (F i ) radiation in percent.

Table 2 .
Sensor specifications and accuracies.The same type of sensors have been deployed on both AWSs.The instrument heights give minimum and maximum values.