3D surface properties of glacier penitentes over an ablation

Penitentes are a common feature of snow and ice surfaces in the semi-arid Andes where very low 10 humidity, in conjunction with persistently cold temperatures and sustained high solar radiation favour their 11 development during the ablation season. As penitentes occur in arid, low-latitude basins where cryospheric water 12 resources are relatively important to local water supply, and atmospheric water vapor is very low, there is potential 13 value in understanding how penitentes might influence the runoff and atmospheric humidity. 14 The complex surface morphology of penitentes makes it difficult to measure the mass loss occurring within them 15 because the (i) spatial distribution of surface lowering within a penitente field is very heterogeneous, and (ii) steep 16 walls and sharp edges of the penitentes limit the line of sight view for surveying from fixed positions and (iii) 17 penitentes themselves limit access for manual measurements. In this study, we solved these measurement problems 18 by using a Microsoft Xbox Kinect sensor to generate small-scale digital surface models (DSMs) of small sample 19 areas of snow and ice penitentes on Tapado Glacier in Chile (30°08’S; 69°55’W) between November 2013 and 20 January 2014. The surfaces produced by the complete processing chain were within the error of standard terrestrial 21 laser scanning techniques. However, in our study insufficient overlap between scanned sections that were mosaicked 22 to cover the studied sites can result in three-dimensional positional errors of up to 0.3 m. 23 Mean surface lowering of the scanned areas was comparable to that derived from point sampling of penitentes at a 24 minimum density of 5 m over a 5 m transverse profile. Over time the penitentes become fewer, wider, deeper and 25 the distribution of slope angles becomes more skewed to steep faces. These morphological changes cannot be 26 captured by the interval sampling by manual point measurements. Roughness was computed on the 3D surfaces by 27 applying previously published geometrical formulae; one for a 3D surface and one for single profiles sampled from 28 the surface. For each method a range of ways of defining the representative height required by these formulae was 29 used, and the calculations were done both with and without using a zero displacement height offset to account for 30 the likelihood of skimming air flow over the closely spaced penitentes. The computed roughness values are in the 31 order of 0.01-0.10 m during the early part of the ablation season increasing to 0.10-0.50 m after the end of 32 December, in line with the roughest values previously published for glacier ice. Both the 3D surface and profile 33 methods of computing roughness are strongly dependent on wind direction. However, the two methods contradict 34 each other in that the maximum roughness computed for the 3D surface coincides with airflow across the penitente 35 lineation while maximum roughness computed for sampled profiles coincides with airflow along the penitente 36 lineation. These findings highlight the importance of determining directional roughness and wind direction for 37 strongly aligned surface features and also suggest more work is required to determine appropriate geometrical 38 roughness formulae for linearized features. 39 1 The Cryosphere Discuss., doi:10.5194/tc-2015-207, 2016 Manuscript under review for journal The Cryosphere Published: 19 January 2016 c © Author(s) 2016. CC-BY 3.0 License.

Mean surface lowering of the scanned areas was comparable to that derived from point sampling of penitentes at a 24 minimum density of 5 m -1 over a 5 m transverse profile. Over time the penitentes become fewer, wider, deeper and 25 the distribution of slope angles becomes more skewed to steep faces. These morphological changes cannot be 26 captured by the interval sampling by manual point measurements. Roughness was computed on the 3D surfaces by 27 applying previously published geometrical formulae; one for a 3D surface and one for single profiles sampled from 28 the surface. For each method a range of ways of defining the representative height required by these formulae was 29 used, and the calculations were done both with and without using a zero displacement height offset to account for 30 the likelihood of skimming air flow over the closely spaced penitentes. The computed roughness values are in the 31 order of 0.01-0.10 m during the early part of the ablation season increasing to 0.10-0.50 m after the end of 32 December, in line with the roughest values previously published for glacier ice. Both the 3D surface and profile 33 methods of computing roughness are strongly dependent on wind direction. However, the two methods contradict 34 each other in that the maximum roughness computed for the 3D surface coincides with airflow across the penitente 35 lineation while maximum roughness computed for sampled profiles coincides with airflow along the penitente 36 lineation. These findings highlight the importance of determining directional roughness and wind direction for 37 strongly aligned surface features and also suggest more work is required to determine appropriate geometrical 38 roughness formulae for linearized features. 39

Introduction 40
Penitentes are spikes of snow or ice ranging from a few centimetres up to several metres in height that can form 41 during the ablation season on snowfields and glaciers under the right conditions. The conditions required for 42 penitentes to form are dew point below 0°C, persistently low air temperatures and sustained strong solar insolation 43 (Lliboutry, 1954). These conditions are frequently met at high elevation, low-latitude glaciers and snowfields, such 44 as in the subtropical Andes (e.g. Hastenrath  Observations show that penitente geometry is aligned with the arc of the sun across the sky and tilted toward the sun 47 at local noon, highlighting the importance of solar radiation in penitente formation (Lliboutry, 1954;Hastenrath and 48 Koci, 1981; Bergeron et al., 2006). Indeed, the alignment and restricted latitudinal range of penitentes (within 55° of 49 the equator on horizontal surfaces) can be explained by solar-to-surface geometry alone (Cathles et al., 2014). The 50 process of penitente growth involves geometric focusing of incident solar radiation by surface irregularities that 51 causes depressions to receive more radiation than surrounding peaks (Amstutz, 1958;Corripio and Purves, 2005;52 Lhermitte et al., 2014; Claudin et al., 2015). Consequently, energy receipts and ablation are initially enhanced in the 53 hollow due to multiple reflection of irradiance, and the surface irregularity becomes amplified. However for 54 substantial penitente growth it is crucial that, at the tips of penitentes, ablation occurs by sublimation and the 55 snow/ice temperature remains below the melting point, while in the troughs between penitentes, melting can occur 56 once a more humid microclimate is established within the hollow (Lliboutry, 1954;Drewry, 1970;Claudin et al., 57 2015). Once the snow/ice in the hollows has reached the melting point, the spatial differentiation of ablation 58 processes serves to further amplify the penitente relief as melting only requires approximately an eighth of the 59 energy of sublimation to remove the same amount of ice. 60 The impact of penitentes on the surface energy balance and ablation of snow and ice is of interest in arid mountains 61 catchments, where penitentes are widespread and meltwater can be a substantial contribution to local hydrological 62 resources (Kaser et al., 2010). Previous studies have shown that penitentes alter the surface energy balance of snow 63 and ice surfaces by reducing effective albedo by up to 40% compared to flat surfaces (Warren et al, 1998;Corripio 64 and MacDonell, 2015). Thus, the presence of penitentes is expected to alter the rate of mass loss and meltwater 67 production of snow and icefields during the ablation season, and, on the basis of the radiative balance it has been 68 postulated that they will accelerate the snow and ice mass loss rates (Cathles et al., 2014). However, the 69 development of penitentes on the surface will also alter the roughness properties in both space and time, but this, as 70 well as its impact on the resultant turbulent fluxes is not quantified. The wind direction-dependence of surface 71 roughness over linearized surface features has been previously observed in wind profile measurements over snow 72 sastrugi, for which the derived aerodynamic roughness length varied from 1-70 mm over 120° range of wind 73 direction (Jackson and Carol, 1978). While penitentes are a relatively rare form of linearized surface feature in many 74 glacierized environments, in contrast linear crevasses are widespread, and although the impact of wind direction on 75 roughness and the resultant turbulent heat fluxes is generally not treated in glaciology, penitentes offer a unique test 76 bed for investigating the significance of linearized features on effective surface roughness for various wind 77 directions. 78 In general, the physical roughness of snow and ice surfaces are particularly prone to varying in space and time (e.g. 79 Smeets et al., 1999;Brock et al., 2006;Fassnacht et al., 2009), it is desirable to be able to replace relatively 80 logistically and technologically challenging methods of determining roughness parameters from atmospheric profile 81 or eddy covariance measurements, with methods based on more readily measurable surface terrain properties (e.g. 82 Kondo and Yamazawa, 1986;Munro, 1989;Andreas, 2011), or properties such as radar backscatter that can be 83 derived from spaceborne instruments (e.g. Blumberg and Greeley, 1993). The most comprehensive surface of 84 methods to determine apparent aerodynamic properties from surface morphometry was carried out by Grimmond 85 and Oke (1999) who tested several methods in urban environments, which are among the roughest surface 86 conditions encountered in boundary layer atmospheric studies. The morphometric estimates of roughness properties 87 were compared with those from aerodynamic methods from numerous field and laboratory studies. Many of the 88 aerodynamic studies were found to be flawed, and the study demonstrates that, despite the considerable effort in 89 obtaining such measurements, their reliability in complex and rough terrain is contested as the computations rely 90 upon theory that is developed for flat homogenous terrain, and in general the aerodynamic results show a similar 91 amount of spread as the various geometrical methods tested. Although, Grimmond and Oke (1999)  Xbox Kinect sensor as a close-range mobile distance ranger to produce a series of small-scale digital surface models 106 (DSMs). These surface models are used to perform (i) the first detailed examination of the geometry of penitentes 107 and how they change over the course of the core ablation season; (ii) an examination of the geometrical roughness 108 properties of penitentes and (iii) compare the volume changes computed from differencing the DSMs with the 109 volume changes estimated from manual measurements of surface lowering within a penitente field. These 110 measurements enable evaluation of how accurately simplified penitente surfaces used in theoretical modelling 111 represent the true surfaces found in nature, improved parameterization of surface roughness in energy balance 112 models applied to glacier and snowfields with penitentes, and the performance of energy balance models over 113 penitente surfaces to be evaluated against mas loss derived from the measured surface changes. 114

Description of fieldsite 116
Tapado Glacier (30°08'S; 69°55'W) lies in the upper Elqui Valley of the semi-arid Andes of the Coquimbo Region 117 of Chile (Figure 1). This glacier is relatively easily accessible and previous research indicates that the glacier surface 118 develops penitentes every summer (Sinclair and MacDonell, 2015). Two separate study areas were analysed. Firstly, 119 a test site was established at a patch of snow penitentes within a dry stream bed at 4243 m a.s.l. in the glacier 120 foreland (Figure 1b). This site was used to (i) trial instrumental setups in order to optimize the field operation of the 121 Kinect sensor, and (ii) compare the performance of the Kinect sensor against a Terrestrial Laser Scanning (TLS)  122 system. This location was chosen due to the logistical difficulties of transporting the TLS to the glacier. 123 Subsequently, two study plots were established at an elevation of 4774 m a.s.l. within the glacier ablation zone. 124 These surfaces at these sites were measured repeatedly using the Xbox Kinect (see section 2.3) during the core 125 ablation season between the end of November 2013 and the beginning of January 2014. An automatic weather 126 station on a free-standing tripod was installed beside the two plots to provide meteorological context for the 127 measurements. 128 The location and layout of the two sites is shown in Figure 1. Site A (5 m by 2 m) was measured four times, on 25 129 November, 11 December, 20 December and 3 January. Site B (2 m by 2 m, Figure 1c) was only measured on the last 130 three dates (Figure 1c). The corners of the study sites were marked with 2 m lengths of plastic plumbing piping 131 hammered vertically into the snow, or drilled into the ice. In order to locate the study sites in space and to provide a 132 common reference frame for each survey date, marker stake positions were measured using a Trimble 5700 133 differential GPS with Zephyr antenna on the 25th November, with a base station in the glacier foreland. On each 134 visit to the glacier, when possible, the stakes were hammered further into the snow and the resultant lowering of the 135 stake top was noted. The maximum standard deviations of the GPS stake positions were < 1.0 cm, 1.1 cm and 136 1.7 cm in easting, northing and elevation respectively, with combined XYZ standard deviation < 2.0 cm for all 137 stakes (Supplement A). Error on the manual measurements of height offsets of the marker stakes on subsequent 138 survey dates is conservatively estimated to be 2.0 cm. This results in total positional errors of the ground control 139 points at each scan date of between 2.3 and 2.7 cm depending on the stake. 140

Terrestrial laser scanning 141
Surface scans of snow penitentes at the test site were undertaken with both a terrestrial laser scanner (TLS) and the 142 Kinect sensor in order to compare the surface scans produced by the well-established TLS method and the relatively 143 new Kinect sensor application. The TLS system used was an Optech ILRIS-LR scanner, which is a long-range 144 terrestrial laser scanner especially suitable for surveying snow and ice surfaces thanks to a shorter wavelength laser and humidity) and trimmed with ILRIS Parser software, aligned with Polyworks IMAlign software into a common 153 local coordinate system and georeferenced with differential GPS measurements using Polyworks IMInspect 154 software. The alignment error of the point clouds as estimated by software is 0.36-0.87 cm and comparison with 155 ground control points gives an error of 5.65 cm. Unfortunately, the scans of snow penitentes could not be carried out 156 with both the TLS and Kinect on the same day, so direct comparison of the TLS and Kinect scans is instead 157 performed on a reference boulder lying on the ground beside the test site, whose surface is assumed unchanged 158 between different scan dates. The TLS scan of the snow penitentes is presented as an example of the nature of the 159 DSM that can be obtained within a penitente field using a TLS ( Figure 2). 160

Kinect scans of surface change 161
The Kinect sensor emits a repeated pattern of structured infra-red (IR) beams, and records the pattern distortion with 162 an onboard IR camera. The depth of field calculation is performed via a proprietary algorithm and a distance map is 163 the raw data output. Using the standard calibration the static raw depth field resolution of the Kinect is 1 mm and the 164 Kinect-measured distance at the center of the field of view is within 1% of the real distance (Mankoff et al., 2013), 165 implying an error of < 1.0 cm at the distance range of the penitente scans. 166 For its original gaming usage, the Kinect is in a fixed position and proprietary software uses feature tracking to track 167 the movements of players moving within the field of view of the Kinect. However, the inverse of this workflow can 168 also be applied wherein the Kinect sensor is moved interactively around a static surface or 3D body, using the same 169 feature tracking to compute the position of the sensor relative to the object and thereby allowing a point cloud 170 reconstruction of the object to be constructed. In this work we apply the second work flow sampling Kinect data 171 using the ReconstructMe™ 2.0 software package. In common with alternative reconstruction packages that are 172 compatible with the Kinect, ReconstructMe™ performs bilateral filtering on the output depth map frame and 173 converts the pixel version of each depth map frame to 3D coordinate maps of vertices and normals. An iterative 174 closest point (ICP) alignment algorithm is then applied frame by frame at three scales to repeatedly rotate and 175 translate the depth field to determine camera position and an aligned surface, giving weighted preference to portions 176 of the surface that are perpendicular to the line of sight. This software has the advantage of producing surface 177 meshes in real-time, so that the operator can visibly check the scan quality and coverage at the time of capture, but 178 the disadvantage that the raw point cloud is not saved and if the real-time tracking is lost a new scan sample must be 179 started. 180 The Xbox Kinect was connected via a 5m powered USB extension cord to an MSI GE60 gaming laptop, powered 181 using a 240V 600W inverter connected to the 12V 160Ah battery of the automatic weather station on the glacier. 182 Scans were carried out by two people; one moving the Kinect across the penitente field and the other monitoring the 183 quality of the surface being generated. The return IR signal of the Kinect is swamped by natural radiation in bright 184 conditions, and this is especially true over bright, rough snow and ice surfaces, which reflect the shortwave 185 radiation, and absorb or scatter much of the longwave radiation signal. To solve this, scanning was carried out at 186 twilight or just after nightfall. Sudden movements caused by the operator slipping or the snow compacting underfoot 187 can result in the ReconstructMe software losing its tracking of common reference points used to generate the 188 continuous surface mesh. Consequently, each study site was scanned in small sections and three to thirteen separate 189 surface meshes were used to cover the area of each study site. 190

Mesh processing 191
Freely available Meshlab software was used to initially align the surface meshes covering each study site using a 192 pairwise alignment procedure. The full mesh processing procedure is presented in Supplement B, and briefly 193 described here. Small surface components, unreferenced and duplicated vertices were removed from the meshes 194 using inbuilt filters. The Meshlab alignment algorithm was applied to objectively optimize the alignment and 195 compute the alignment error. This alignment procedure uses an ICP algorithm to iteratively align the component 196 meshes and distribute the alignment errors evenly across the mosaicked surface mesh. Alignment solutions 197 consistently had mean distributed error < 4 mm (Supplement B). The aligned meshes were flattened into a single 198 layer, remeshed using a Poisson filter and finally resampled to reduce the point density by setting a minimum vertex 199 spacing of 2.5mm. 200 The surface mesh for each scan date was georeferenced using the known coordinates of the base of the marker 201 stakes at the time of each scan because the upper portions of the stakes are often poorly represented in the scans due 202 to the fact that ReconstructMe™ does not handle symmetrical objects well. It proved difficult in some cases to 203 locate the surfaces in space such that the locations of all marker stakes were consistent with the ground control 204 points. This is most likely an artifact of a combination of (i) reduced mesh quality at the margins of the component 205 scans and (ii) insufficient overlap between some scan sections producing distortion within the mesh alignment. The 206 mismatch evident in the georeferencing step (Table 1) is much larger than the mesh alignment error (Supplement B). 207 To eliminate the marker stakes and any data gaps near the margins of the study areas, each surface mesh was sub-208 sampled within the staked area. The sub-sampled area for site A is a 2.0 by 3.5 m horizontal area (7.00 m 2 ), and site 209 B is a 1.5 x 1. upward-facing triangles were computed column-wise from these projected areas using the height coordinate of the 217 triangle centroid as the height dimension for each column. These were summed and volumes for overhanging 218 triangles, calculated in the same way as the up-ward facing volumes were subtracted to derive a total volume 219 between the reference surface and the scanned penitente surface. Successive volumes were subtracted to obtain the 220 volume change over each measurement interval. 221

Calculations of geometric surface roughness 222
The aerodynamic roughness length (z 0 ) is the distance above the surface at which a logarithmic windspeed profile 223 under neutral conditions would be extrapolated down through the surface layer and reach zero. Over taller roughness 224 elements the level of action of momentum transfer between the airflow and the surface roughness elements is 225 displaced upwards by a distance, termed the zero-plane displacement (z d ). Above particularly rough surfaces, a 226 roughness sub-layer is formed in the lowest part of the surface layer within which surface roughness elements create 227 a complex 3D flow that is almost chaotic. Where roughness elements are widely spaced, the separated flow over 228 obstacles reattaches to the surface before the subsequent obstacle is reached. More closely packed roughness 229 elements experience a wake interference regime, and in the most densely packed arrays of roughness elements 230 skimming flow occurs (Grimmond and Oke, 1999). At the top of the roughness sublayer individual wakes caused by 231 surface obstacles are smeared out and the flow is independent of horizontal position, and thus, observations at this 232 level represent the integrated surface rather than individual surface obstacles. This level is known as the blending 233 height (z r ). All these properties are dependent on the size and arrangement of surface roughness elements. 234 There are a number of formulations for deriving z 0 from geometrical measurements. For example, the simplest 235 approach is to take the standard deviation of the surface elevations as a measure of roughness (Thomsen et al., 236 2015). In this work, the surface meshes were analysed for roughness on the basis of a widely-used relationship 237 established by Lettau (1969), initially developed for isolated, regular obstacles distributed over a plane: 238 where h is the height of the obstacles, s is the upwind silhouette area of each obstacle and S is the specific area 240 occupied by each roughness element obstacle, also referred to as its lot area. The roughness values computed using 241 Equation 1 over 3D snow surfaces has been shown to vary widely depending on the methods of surface interpolation 242 used (Fassnacht et al., 2014), due to the influence on interpolation method on the unit surface area occupied by each 243 roughness element. However in this work the high resolution meshes used can be expected to adequately capture the 244 surface properties as no extrapolation or interpolation procedure is needed. Isolated roughness elements of regular 245 geometry distributed over a horizontal plane are a poor analogy for the irregular surface topography of a penitente 246 field, and the applicability of this formulation over penitentes has not been established. Nevertheless, we apply the 247 analysis as an illustration of the nature of the results generated from such an approach over penitentes and hope that 248 future aerodynamic roughness lengths can be compared to these geometrically derived ones. Macdonald and others 249 (1998) state that for irregular obstacles h can be replaced by average obstacle height, s with the sum of all the 250 upwind silhouette areas, and S with the total area covered by the obstacles. While the upwind silhouette area, and 251 indeed surface area in any direction, is relatively easily defined for each surface mesh area using trigonometry, it is 252 difficult to define individual roughness elements and their representative heights, due to the lack of an apparent base 253 level. Here we first detrend the surfaces to remove any general surface slope at the site, then compute the roughness 254 for the detrended 3D meshes assuming that the roughness elements cover the whole surface area (i.e S = plot area), 255 and for four possible representations of average obstacle height (h) as follows: (i) the maximum range of the 256 detrended mesh; (ii) twice the standard deviation of the detrended surface mesh; (iii) mean mesh height above the 257 mesh minimum; and (iv) median mesh height above the minimum. element density exceeds 20-30%, as is expected for penitente fields (Macdonald et al., 1998). High density 260 roughness elements means that they interfere with the airflow around each other, and the zero wind velocity level is 261 displaced upwards, and effective roughness is a result of the roughness elements above this zero velocity 262 displacement plane. The zero displacement height in this sense, gives an indication of the penetration depth of 263 effective turbulent mixing into the penitente field. Accordingly, we additionally present sample calculations of 264 three-dimensional roughness on the detrended surface meshes using three possible realizations of z d , as z d is also 265 unknown in the case of the penitente fields being sampled. In the first case, z d is taken to be h, in the second 2/3 h, 266 which is a widely used standard in forests and other complex terrain applications (Brutseart, 1975), and in the third 267 1/3 h for comparison, both computed for the four realizations of h used as before. Equation 1, (for irregular 268 obstacles) is then applied to the roughness elements remaining above the plane of the general surface slope offset by 269 a distance z d above the minimum height of the surface mesh. The representative height h for this portion of the mesh 270 exceeding the plane is taken to be the mean area-weighted height of all triangles above this plane, s is the summed 271 frontal area of all mesh triangles above z d that face into the chosen wind direction and S is the total horizontal area of 272 the surface components above z d . 273 Munro (1989Munro ( , 1990) modified the formula of Lettau (1969) to be applied to a single irregular surface cross-section 274 of length X, sampled perpendicular to the wind direction. This modified formulation is easier to work with on a 275 glacier where the roughness elements are irregular, closely spaced, and generally poor approximations of objects 276 distributed over a plane. Instead of having to define an obstacle height above the plane, h is replaced with an 277 effective height h* expressed as twice the standard deviation from the standardized mean profile height; s is replaced 278 with h*X/2f, in which f is the number of profile sections that are above the mean elevation; and S is replaced with 279 (X/f) 2 . This approach approximates the surface elevation profile as rectangular elements of equal size, and has been 280 shown to give results within 12% of the silhouette area determined by integrating between true topographic minima 281 (Munro, 1989). Importantly, roughness values derived this way over snow, slush and ice surfaces show reasonable 282 agreement with roughness values derived from wind profiles (Brock et al., 2006). To investigate the nature of the 283 roughness computed this way for north-south and east-west impinging wind directions, cross profiles longer than 284 1.5 m at 0.1m intervals orientated E-W and N-S were extracted from each scanned surface. Cross-sections were 285 detrended to remove the influence of any general surface slope at the site, and roughness was computed following 286 the modifications of Munro for each detrended surface profile. Mean profile roughness for these two wind directions 287 are presented for each sampled surface. 288

Manual measurements of surface change 289
Traditional stake measurements of glacier surface lowering made at a single point are unreliable within the 290 inhomogeneous surface of a penitente field, as multiple measurements are required to characterize the complex 291 surface. One alternative is to measure surface lowering at intervals along a profile perpendicular to the main axis of 292 alignment of the penitentes. Such a reference was installed along the 5 m-long eastern margin of site A, between two 293 longer corner stakes drilled 3 m into the ice using a Kovacs hand drill. The distance between a levelled string and 294 the glacier surface was measured using a standard tape measure at 0.2 m intervals on 23 November. Subsequent 295 measurements, on the 12 and 21 December and on 4 January, were made at 0.1 m intervals. All measurements were 296 recorded to the nearest centimetre, and the error on each measurement is conservatively estimated to be 2.0 cm, 297 which is assumed to capture the error associated with the horizontal position of the measurements along the 298 reference frame and the vertical measurements of the distance to the surface beneath.

Evaluation of the quality and suitability of penitente scans by TLS and Kinect 301
At the test site, the snow penitentes were well-developed and between 0.5 and 1.0 m in height (Figure 1b). TLS 302 scans were made of these penitentes to illustrate the capabilities of this more conventional scanning system in 303 capturing the penitente surfaces. TLS scans were taken from five different vantage points positioned above the 304 penitentes. The penitente surface produced by the TLS had surface slope ranging between -30 and 90 degrees, 305 indicating that overhanging surfaces within the penitente field are captured, however only 58% of the total surveyed 306 horizontal area could be scanned as the deepest parts of the troughs were obscured from the view of TLS by the 307 surrounding penitentes (Figure 2a). By comparison, the hand-held, mobile nature of the Kinect means that 100% of 308 the surface of the penitente field can be captured as the field of view can be adjusted into limitless close-range 309 positions. The long range of the TLS makes it easier to cover large areas in comparison to the close range Kinect 310 sensor, but as only penitente tips are scanned the utility of this larger areal coverage is limited. 311 The Kinect scan of the reference boulder produced from three mosaicked meshes was aligned to that produced from 312 TLS point clouds. The TLS scan was incomplete, with parts of the top and overhanging surfaces of the boulder 313 missing due to being obscured from the TLS survey positions, while the Kinect scan achieved complete coverage of 314 the boulder. The difference between the two aligned meshes where overlapping data existed was always < 2 cm 315 (Figure 2b), which is well within the error of the georeferenced TLS surface model. Larger differences in Figure 2b positions are all < 2.0 cm. However, in this study the three-dimensional georeferencing error is large (Table 1)  322 compared to the other sources and can be taken as a reasonable value for the error of the total process chain. Errors 323 given on the seasonal mass, volume and surface changes are based on summing the squares of the mean elevation 324 difference between the marker stakes and ground control points (GPCs) at each site on the first and last survey dates. 325

Meteorological conditions 326
During the study period one significant snowfall event occurred on the 8 th December 2013, when the sonic ranger 327 recorded an increase of surface height of 0.09 m over the course of the day, and temperature and incoming longwave 328 radiation increase progressively ( The morphometry of the sampled penitentes changed visibly over the measured intervals (Figures 3 and 4). The 356 strong east-west preferential orientation predicted from theory developed early and was maintained throughout study 357 period. The expression of this alignment is more convoluted in the stages of development studied here than the 358 parallel rows of penitentes used in model representations ( From the onset of measurements the surface aspect distribution is strongly dominated by north and south facing 369 components and this becomes more pronounced in the latter measurements and the preferred orientation rotates 370 slightly over the course of the season (Figure 4 e & f). 371

Surface roughness assessments 372
Given that aerodynamic measurements to determine the most suitable representative height and zero displacement 373 level for penitentes are thus far unavailable, the approach taken here was to do an exploratory study and compute 374 geometric surface roughness values using various ways of expressing h and z d . As a consequence the results are 375 purely illustrative and while patterns can be drawn from them that have meaning for understanding the nature of the 376 computation, the applicability of these values in turbulent exchange calculations remains to be established. The 377 representative height, h, used in the calculations increases over time in all cases, and is bounded by the case taking h 378 to be the range of the detrended surfaces (maximum) and the case taking h as twice the standard deviation of the 379 detrended surface ( Figure 5). For clarity, the other two case values are not included in the plots shown here. 380 Differences within a single method between the two sites can reach as much as 0. the way in which the representative height is expressed, the time of year and the wind direction (Figure 7). However, 389 given the close spacing of the penitentes it seems appropriate to also explore what the calculated z 0 would be like 390 when applying a zero displacement height offset, although again, in the absence of validation data these numbers 391 can be only indicative of the pattern of roughness computed by these methods. Introducing the zero displacement 392 height reduces the maximum calculated roughness by about half, and also reduces the variability between different 393 representative heights (Figure 7), as a smaller h value translates into a smaller z d so that the calculation is performed 394 on a larger portion of the mesh. 395 Surface roughness assessments on the basis of calculations following Munro's modification for single profile 396 measurements were applied to cross profiles longer than 1.5 m yielding 20 (6) profiles orientated N-S and 33 (7) E-397 W at site A (B). Surface amplitude increases over time, and the amplitude of the N-S running cross profiles is 398 generally larger than the E-W running cross profiles, as illustrated in the example of site B (Figure 8). The profile-399 computed roughness length increases monotonically over time at site B, but shows a reduction over the first period 400 at site A, associated with snowfall during this period. Both the range and relative increase in roughness over time is 401 larger for the N-S running profiles. The computed roughness at both sites is 4.3 to 6.8 times larger for airflow 402 impinging on the penitente field in an E-W direction than for airflow in the N-S direction. This is contrary to the 403 results computed on the full 3D mesh surface, but is understandable because this formulation relies on the amplitude 404 of the surface, which is generally larger in the N-S orientated cross profiles than the E-W running cross profiles. 405 Prevailing wind direction differs only slightly in each period with an increasing northwesterly component in the 406 second two periods compared to the first. This may be related to the occurrence of snow during the first period, 407 which can be expected to alter the thermally driven valley wind systems. Over the whole study period wind direction 408 is predominantly from the south-westerly sector, but swings through southerly to easterly thereby encompassing 409 both extreme wind angles used in the roughness calculations here (Figure 9). This indicates that the effective 410 roughness can be expected to differ significantly depending on the wind direction. 411

Manual measurements of reference cross-profile 412
Using data sampled at 0. The underlying ice surface identified by the snow probing, does not appear to be influencing the structure of the 437 snow penitentes developing in the current season. However, it is difficult to draw a firm conclusion based on 438 measurements at only 0.2 m spacing, particularly as, while the surface of the penitentes was still snow on the 3 439 January, in several instances the surface had lowered below the level of the ice interface indicated by the initial 440 probing. 441

Methods of measuring change of rough glacier surface elements 443
The test site for scanning penitentes with a TLS was chosen as it provided the most optimal viewing angles possible 444 from scanning positions, as the penitentes lay in a river bed and scanning positions could be established on the 445 surrounding river banks to look down into the penitente field. Nevertheless, the terrestrial laser scanning could only 446 capture the tips rather than the whole surface of the penitentes and, as ablation is at its maximum in the troughs, TLS 447 data is therefore not able to determine the true volume change ongoing in penitentes. The coverage would be 448 increased if a higher viewing angle could be achieved, but the steep, dense nature of penitente fields makes it 449 difficult to imagine where sufficient suitable locations can be found surrounding glaciers or snowfields with 450 penitentes. In contrast, the mobile Kinect sensor can be moved across the complex relief of the penitente field to 451 make a complete surface model. Although it is in principle possible to capture a large area with the ReconstructMe 452 software used here, and it offers the advantage of providing real time feedback on the mesh coverage, it proved 453 difficult to capture the study sites in a single scan given (i) the reduced signal range of the sensor over snow and ice 454 (Mankoff et al., 2013), and (ii) the difficulty of moving around the penitente field. As a result, partial scans were 455 obtained, with the disadvantage that subsequently combining these introduces a substantial degree of additional error 456 associated with alignment if the component scans were not of high quality at the margins, or did not overlap 457 adjacent scan areas sufficiently. A combination of these two techniques might allow the extrapolation of small-scale 458 geometry changes and volume loss determined from a Kinect surface scan to be extrapolated usefully to the glacier 459 or snowfield scale. 460 Despite not visually capturing the complex surface properties of the penitentes, manual measurements of surface 461 height change in a penitente field along a profile cross-cutting the penitentes are robust for determining mean 462 surface lowering rates, and show good agreement to the volume changes computed from differencing the digital 463 surface models scanned in detail using a Kinect. Thus, the detailed surface geometry need not be known in order to 464 reasonably calculate the total volume loss over time within penitente fields. Comparison of the manual sampling at 465 different intervals suggest that five samples per meter is adequate to characterize penitentes. Over the 39 days of the 466 study, the mass loss calculated from 26 points spaced at 0.2 m intervals along a 5 m profile crosscutting the 467 penitentes differed from that calculated from volume change computed on surface meshes consisting of over 1.3 468 million points and covering an area of 7 m 2 by only 28 kg m -2 . Although this difference was within the error of the 469 two measurement types, the seasonal difference, assuming that this difference applies to a whole ablation season of 470 120 days would be 86 kg m -2 , and applied to the whole glacier (3.6 km 2 ) would amount to an underestimate of mass 471 loss over an ablation season of 0.3 gigatonnes. As a side note, the probing of snowdepth carried out as part of this 472 study highlights the difficulty in identifying the underlying ice surface, or summer ablation surface, in this way 473 within a penitente field, suggesting that a single location must be sampled very densely to obtain a characteristic 474 snowdepth in this way. 475

Penitente morphometry and change in time 476
The manual measurements at 0.2 m intervals are adequate to determine the mean surface lowering within a penitente 477 field, giving confidence to this type of simplified measurement on seasonal timescales. However, the interval 478 measurements cannot capture the surface morphometry, or how it changes in time. 479 At all times the penitente surface represents a much larger total surface area than the equivalent non-penitente 480 surface. Over time the surface relief, and slope angle, increases as the penitentes deepen, unless a snowfall event 481 occurs to partially fill the troughs, which also reduced the mean surface slope. The control of solar radiation on 482 penitente morphology means that the vast majority of the surface consistently dips steeply to the north and south at 483 all stages of development. This means that the angle of incidence of direct solar radiation is reduced, decreasing 484 both penitente field compared to a flat surface, the total absorbed shortwave is a third higher in the modeled penitentes. 489 At Tapado Glacier, penitentes are initially overhanging to the north, and the southfacing sides are convex compared 490 to the northfacing overhanging faces. Over the season the penitentes become more upright as the noon solar angle 491 gets higher. Idealized modelling based on measurements at Tapado Glacier, shows that concave and convex slopes, 492 as well as penitente size have been shown to impact the apparent albedo as measured by ground and satellite sensors 493 (Lhermitte, et el., 2014), and there may be some value in assessing the impact of these morphometry changes on 494 albedo over time. For the idealized penitente surface at 33°S during summer solstice case, modeled increase in net 495 shortwave radiation over penitentes is not compensated by modelled changes in net longwave radiation, meaning 496 that the excess energy receipts must be compensated by either turbulent energy fluxes or consumption of energy by 497 melting (Corripio and Purves, 2005). 498 In the context of the numerical theory of Claudin and others (2015), progressive widening of the penitente spacing, 499 as observed at both site A and B, is indicative of changes in the atmospheric level at which water vapor content is 500 unaffected by the vapor flux from the penitente surface. Simultaneous field or laboratory measurements of penitente 501 spacing evolution and vapor fluxes above the surface would be required to solidly confirm this, but the field 502 measurements provided here can be used as an indication of the level to which vapor flux from the surface is 503 influencing the boundary layer vapor content. 504

Surface roughness 505
In this work a single, simple geometric relationship (Lettau 1969) was investigated because a profile-based version 506 of this formulation has previously been tested against aerodynamic measurements over glacier surfaces (Munro, 507 1989(Munro, 507 , 1990Brock et al., 2006). Certainly other relationships could be explored in the context of linearized glacier 508 features, but given the wide spread of values produced in previous comparisons such an analysis might be of limited 509 value in the absence of simultaneous aerodynamical investigations (Grimmond and Oke, 1999). Furthermore, the 510 results of Grimmond and Oke (1999) indicate that for the cities sampled, the Lettau method gives z 0 values that are 511 in the middle of the range of all the methods. The analysis of geometric computations of roughness properties in 512 Grimmond and Oke (1999) highlight the importance of correctly determining z d , and limited sensitivity analyses 513 show the computed z d and z 0 to be strongly dependent on the dimensions of the obstacles. Lettau's (1969)  point of a logarithmic wind profile to be lower than formulations that include z d in their computation of z 0 . In this 518 work however we computed z d in a separate preceding step to explore the impact of z d on the computed z 0. 519 As penitentes fields present very densely packed roughness elements, the frontal area of the surface tends to be large 520 compared to the ground area, and the limits of the ratio of frontal to planar area found in this study implies that 521 skimming flow is almost always occurring over penitente fields, such that turbulent airflow in the overlying 522 atmosphere does not penetrate to the full depth of the penitente fields. This is in agreement with the theory of 523 formation and growth of penitentes in which the development and preservation of a humid microclimate within the 524 penitente hollows is required to facilitate differential ablation between the trough and tip of the penitente. As the 525 spacing between the penitentes also increases over the ablation season the features become less densely packed over 526 time, although the available data are insufficient to determine if the spacing increases sufficiently by the end of the 527 season to comply with the applicable limits of the roughness calculation used here. ice (Smeets et al., 1999;Obleitner, 2000). The topographic analysis clearly shows that in the absence of intervening 534 snowfall events, this roughness increase is related to the deepening of the penitentes over time and an increase of the 535 surface amplitude. The patterns of the computed roughness properties is consistent between the two neighbouring 536 sites, but individual values can differ, suggesting that local relief varies substantially and sampling a larger area 537 would be beneficial in order to capture mean properties. 538 The strong alignment of penitentes means that roughness calculated on the 3D surface meshes is higher for wind 539 impinging in a north-south direction as the large faces of the penitentes form the frontal area in this case. In contrast, 540 if roughness is computed for individual profiles extracted from the mesh to mimic manual transect measurements in 541 the field, roughness is between 3 and 6 times larger for air flow along the penitente lineation (E-W) than it is across 542 the lineation (N-S). While clearly highlighting that the surface roughness of the strongly aligned penitente fields is 543 dependent on wind direction, this contradiction poses a conundrum as neither approach has been specifically 544 evaluated against independent surface roughness derived from atmospheric profile measurements over penitentes. 545 Consequently, although surface roughness calculations on the basis of profile geometry have been evaluated against 546 aerodynamic roughness over rough ice surfaces, the available data is insufficient to distinguish which pattern is 547 more appropriate for calculating turbulent fluxes over penitentes. It principle it sounds reasonable to expect airflow 548 across the penitente lineation to maximize turbulence as the penitentes present a large surface area to the wind, yet, 549 if skimming flow is established, with the result that only the tips of the penitentes are determining the structure of 550 the turbulence then roughness in this direction would be strongly reduced, and perhaps even be less than for air flow 551 along the penitente lineation, for which the smaller frontal area reduces the likelihood of skimming flow. Further 552 investigation of this in order to quantify the impact of penitentes on turbulent fluxes for various airflow patterns 553 requires measurement of turbulent fluxes using eddy covariance or atmospheric profile methods, which would 554 demonstrate the nature of the directional roughness and establish the impact of penitentes on turbulent energy fluxes 555 for different wind directions. Such measurements would be best implemented in a manner which can sample all 556 wind directions equally, and eddy covariance systems for which analysis is limited to a sector of airflow centred 557 around the prevailing airflow source, might not be able to capture the nature of the directional dependence correctly. 558 In this study we did not explicitly compute the blending height as available formulae are dependent upon z 0 and z d . 559 Estimates penitentes would have to be carried out at some height above the surface to capture mean surface properties rather 564 than the effects of individual roughness elements. The mathematical model of Claudin and others (2015) indicates 565 that the level at which the vapour flux does not is constant in horizontal space, and therefore is the product of mean 566 surface properties, is related to the spacing of the penitentes. Taking this to be representative of the blending height 567 would imply that a formulation for the blending height might be possible on the spacing of penitentes alone, and that 568 this in turn might contain useful data for understanding the structure and efficiency of turbulence above penitentes. 569 However, exploring these ideas requires information from meteorological measurements as well as the geometrical 570 information offered in this paper. 571

Conclusion 572
Surface scanning technology and software is an area of rapid development, and a number of potentially superior 573 alternative set-ups and data capture sensors and software is now available. This study demonstrates that the 574 Microsoft Kinect sensor can work successfully at close range over rough snow and ice surfaces under low light 575 conditions, and generate useful data for assessing the geometry of complex terrain and surface roughness properties. 576 The data collected offers the first detailed study of how the geometry of penitentes evolves through time, 577 highlighting the rate of change of surface properties over an ablation season that can serve as a guideline for 578 parameterizing surface properties required for energy and mass balance modelling of penitente surfaces. The 579 measurements confirm that even relatively crude manual measurements of penitente surface lowering are adequate 580 for quantifying the seasonal mass loss. 581 Aerodynamical roughness properties and related metrics over very rough surfaces remain poorly quantified and both 582 geometric and meteorological determinations of these values show a wide spread; consequently it remains unclear 583 what the best methods to use are or what values modellers would be best to use (Grimmond and Oke, 1999). In this 584 context penitentes and further study of them offers a useful opportunity as (a) their morphometric evolution over 585 time allows various geometries to be evaluated by instrumenting and scanning a single site, and (b) they offer a 586 bridge between wind tunnel and urban field experimentation of turbulence and roughness over extreme terrain. 587 Although validity of surface roughness calculations based on surface geometry remains to be established for 588 penitentes, this study highlights that (i) skimming flow is expected to persist over penitentes field, but is more likely 589 under wind directions perpendicular to the penitente alignment; (ii) z d is certainly greater than zero, and while the 590 depth of penetration of surface layer turbulence into a penitente field is not clearly established it is likely to evolve 591 with the developing penitentes, and values of z d ~2/3h give results that are theoretically reasonable in the framework 592 outlined by Grimmond and Oke (1999); (iii) the two methods of geometric computation of surface roughness 593 applied here give conflicting results as to whether the effective surface roughness of penitentes is greater for airflow 594 along or across the penitente lineation and (iv) more complete understanding of the impact of penitentes on the 595 turbulent structure, its evolution in time, and its directional dependency, would require atmospheric measurements 596 with no directional bias. 597 Potential future applications and analyses of the surfaces generated in this study include (i) using surface properties 598 and roughness values as a guide for input into surface energy balance models; (ii) assessing the performance of 599 models against the measured volume loss over time and (iii) evaluating how well simplified representations of 600 penitente surfaces used in small scale radiation models and turbulence models capture the real-world complexity. 601 Such studies would help establish the nature of the likely micro-climatic distribution of the surface energy balance 602 within a real penitente field, and as a result the impact of penitentes on runoff and exchange of water vapour with 603 the atmosphere.           Munro (1989Munro ( , 1999.