A multi-season investigation of glacier surface roughness lengths through in situ and remote observation

The roughness length values for momentum, temperature, and water vapour are key inputs to the bulk aerodynamic method for estimating turbulent heat flux. Measurements of site-specific roughness length are rare for glacier surfaces, and substantial uncertainty remains in the values and ratios commonly assumed when parameterising turbulence. Over three melt seasons, eddy covariance observations were implemented to derive the momentum and scalar roughness lengths at several locations on two mid-latitude mountain glaciers. In addition, two techniques were developed in this study for the remote estimation of momentum roughness length, utilising lidar-derived digital elevation models with a 1× 1 m resolution. Seasonal mean momentum roughness length values derived from eddy covariance observations at each location ranged from 0.7 to 4.5 mm for ice surfaces and 0.5 to 2.4 mm for snow surfaces. From one season to the next, mean momentum roughness length values over ice remained relatively consistent at a given location (0–1 mm difference between seasonal mean values), while within a season, temporal variability in momentum roughness length over melting snow was found to be substantial (> an order of magnitude). The two remote techniques were able to differentiate between ice and snow cover and return momentum roughness lengths that were within 1–2 mm ( an order of magnitude) of the in situ eddy covariance values. Changes in wind direction affected the magnitude of the momentum roughness length due to the anisotropic nature of features on a melting glacier surface. Persistence in downslope wind direction on the glacier surfaces, however, reduced the influence of this variability. Scalar roughness length values showed considerable variation (up to 2.5 orders of magnitude) between locations and seasons and no evidence of a constant ratio with momentum roughness length or each other. Of the tested estimation methods, the Andreas (1987) surface renewal model returned scalar roughness lengths closest to those derived from eddy covariance observations. Combining this scalar method with the remote techniques developed here for estimating momentum roughness length may facilitate the distributed parameterisation of turbulent heat flux over glacier surfaces without in situ measurements.


Introduction
The turbulent fluxes of sensible and latent heat can form a major component of the surface energy balance (SEB) of a glacier, and substantially influence its rate of surface melt (Hock and Holmgren, 1996;Anderson et al., 2010;Fitzpatrick et al., 2017). With a lack of direct measurement on glaciers, the bulk aerodynamic method is commonly used to parameterise the turbulent fluxes, requiring input of roughness length values for momentum ( 0 ), temperature ( 0 ), and water vapour 5 ( 0 ). Observations of roughness length are rare on glacier surfaces, however. The majority of SEB studies use values and ratios from previous research on similar surface types (e.g. Gillet and Cullen, 2011;Giesen et al., 2014), or treat roughness lengths as model tuning parameters (e.g. Braun and Hock, 2004;Sicart et al., 2005), rather than obtaining site-specific measurements. This approach introduces uncertainty into turbulent flux estimation, as the transferability of roughness lengths between locations and seasons is unknown. Furthermore, parameterisation of the turbulent heat fluxes has been 10 shown in previous studies to be highly sensitive to the implemented roughness lengths (up to a doubling of the calculated flux for one order of magnitude increase in 0 ), and to dominate over stability corrections as a source of uncertainty (Munro, 1989;Braithwaite, 1995;Brock et al., 2000;Fitzpatrick et al., 2017). The importance of accurate roughness length selection, as identified in these studies, highlights the need for further research on the spatial and temporal variability of their values on glacier surfaces, and on the methods used in their estimation . 15 The roughness length values are defined as the lower limits of integration in the bulk-gradient or 'K' theory parameterisation of the turbulent fluxes (Stull, 1988). 0 can be thought of as the height above the surface at which wind speed, extrapolated downwards along an assumed logarithmic profile, will reach its surface value. Similarly, 0 and 0 can be considered to be the heights at which temperature and specific humidity reach their surface values, respectively. 0 accounts for the effects of form drag on the near-surface wind profile due to the interaction of airflow with features on the surface. In many glacier 20 studies and climate models (e.g. Van As, 2011;Fausto et al., 2016), 0 values of 1 mm and 0.1.mm are used for ice and snow surfaces, respectively, and are often assumed constant with time. Where measurements have been obtained on glacier surfaces, however, a large range of 0 values have been recorded, with several orders of magnitude of variation between different glaciers and seasons (e.g. Van den Broeke et al., 2005;Brock et al., 2010). In addition, existing values for 0 on glaciers (observed through mast-based vertical wind profile measurements, or estimated from eddy covariance (EC) 25 observations or microtopography surveys) only provided values for an individual location or turbulent footprint. Implementing these single values in a glacier-wide distributed model or in a point model at another location on the glacier may not account for the potential variability in surface roughness that may exist across a glacier surface (e.g. Smith et al., 2016). The turbulent footprint referred to above is the source region for the turbulent fluxes received at a given location. It represents the upwind area that influences and contributes to the observed fluxes, and hence, the surface properties that 30 modulate turbulence generation. Broadly speaking, the turbulent footprint for fluxes measured at a given height will extend upwind by a distance of roughly 100 times the measurement height (Burba, 2013). Efforts have been made in previous boundary-layer studies over different land surfaces to determine momentum roughness length values for large areas, including over forestry, scrubland, and outwash plains (e.g. Nield et al., 2013;Li et al., 2017).
A range of remote sensing techniques have been implemented in such studies, including the use of light detection and ranging (LiDAR) systems. Paul-Limoges et al. (2013) used digital elevation models (DEMs), obtained from airborne LiDAR, to estimate 0 values over a harvested forest surface ( 0 = 0.13 m), and found good agreement with corresponding 5 EC-derived values ( 0 = 0.12 m). Similar studies on mountain glaciers are extremely rare. Smith et al. (2016) used terrestrial-based structure-from-motion photogrammetry and laser surveying to generate a distributed map of 0 estimates for a glacier. Meteorological-based evaluation of the returned 0 estimates was not carried out, however.
The scalar roughness lengths ( 0 and 0 ) are commonly estimated in SEB studies using a fixed ratio with 0 , and are generally assumed to be equal to or one to two orders of magnitude smaller than the momentum roughness length (e.g. Hock 10 and Holmgren, 2005;Sicart et al., 2005;Hoffman et al., 2008). Molecular diffusion controls the rate of scalar transfer with a surface, and having a smaller spatial scale than the form drag processes driving momentum transfer, it is likely that the scalar roughness lengths would be smaller (Beljaars and Holtslag, 1991). The persistence of this ratio with time is uncertain, however. Surface renewal methods have been implemented in some studies (e.g. Andreas, 1987;Smeets and van den Broeke, 2008), where variation in this ratio is described as a function of the roughness Reynolds number * . Changes in 15 mean air temperature and relative humidity have also been proposed as drivers of scalar roughness length variation (e.g. Calanca, 2001;Park et al., 2010).
The initial goal of this study is to obtain in situ values of the momentum and scalar roughness lengths from multiple locations over several seasons. EC-observed data will be implemented into the bulk aerodynamic method to derive these values. The temporal variability of roughness lengths on a glacier will be examined, and the transferability of values between 20 location and years will be assessed. Commonly assumed values and ratios from the literature will be compared with the obtained data, and predictive relationships for the scalar roughness lengths will be tested. The second goal of this study is to develop remote methods for estimating momentum roughness lengths for a glacier surface, which would facilitate SEB modelling for glaciers without in situ observations, and distributed modelling for glaciers with point measurements only.
Digital elevation models will be obtained for each study location and will be used to provide surface height data for the two 25 roughness methods developed in the study. Turbulence footprint modelling will be employed in one of these methods to identify the region of the glacier surface influencing the EC-derived roughness length values. The estimates from both remote methods will then be compared with those from corresponding in situ observations.

Field Campaign
Observations were carried out over three melt seasons (2014-2016) on two glaciers in the Purcell Mountains of British Columbia, Canada (Fig. 1). Nordic Glacier (51°26' N, 117°42'W) is a small (~5 km 2 ), north-facing glacier, ranging in elevation from 2,000 m to 2,900 m above sea level (a.s.l.), approximately. An automatic weather station (AWS) was 5 installed in the ablation zone of the glacier through July and August 2014 (NG14). Conrad Glacier (50°49' N, 116°55'W) is located 87 km to the southeast of Nordic, with an area of ~15 km 2 , and an elevation range of 1,800 m to 3,200 m a.s.l., approximately. A total of four AWS deployments were executed on Conrad during 2015 and 2016; two stations in the ablation zone from July to September 2015 (CG15-1 and CG15-2), and one in both the ablation (CG16-1) and accumulation (CG16-2) zones from June to August 2016 (Table 1). An exposed ice surface was present during observations at NG14, 10 CG15-1, CG15-2, and for most of the observation period at CG16-1, while a snow surface was present throughout at CG16-2, and for the first 10 days at CG16-1. A transitional snow surface was present for the first four days at NG14, with partial snow cover diminishing to a fully bare ice surface.

AWS
The AWS developed for this project (see Fitzpatrick et al., 2017) was equipped with an array of meteorological and 15 glaciological sensors to observe the complete SEB, with additional sensors added to the stations each year (Table 2). Open and closed path eddy covariance (OPEC and CPEC) systems were used in this project to observe the turbulent heat fluxes, with both forms installed on the same station, in some cases (CG15-1 and CG16-1). Both systems were comprised of a 3D sonic anemometer, and an infrared gas analyser; the OPEC analyser has a sample space that is open to passive air flow, while the CPEC analyser has a closed sample space into which air is drawn using a pump. Implementing these methods 20 together helped minimise gaps in the turbulence dataset (OPEC analysers are susceptible to errors during precipitation), and enabled a comparison of their values and performance in a glacial environment. The EC data was recorded in raw 20Hz format, with observations from the remaining sensors stored in one-minute averages.
The meteorological sensors were housed on a four-legged quadpod, which provided a stable platform (verified by an inclinometer sensor) that lowered as the ice melted, and maintained a constant height of the sensors above the surface. EC 25 measurements were carried out at a constant height (~2 m at each station) to avoid substantially varying the turbulence footprint area and to reduce the risk of elevating the sensor above the turbulence coupled with the surface (Burba, 2013;Aubinet, 2008). The installation site for each station was selected based on the criteria of a relatively uniform upwind footprint and slope angle, so as to minimise the corrections required in the EC (and radiation) data processing. The EC systems were installed on the upslope side of each station, so as to be the first point of contact with the prevailing wind 30 (downslope), and to help minimise flow distortion. Time lapse cameras at each location were used to observe the surface and atmospheric conditions over a season, and to monitor station behaviour. Over the three melt seasons, the stations performed well, operating continuously over each study period. The solar power systems for the stations had been designed to have sufficient battery storage for approximately a week of operation without sufficient recharge (due to persistent overcast conditions or covering of the solar panels by snow/ice.). If battery voltages dropped below a critical level, the system was designed to restrict power supply to the higher consuming sensors (e.g. CPEC system) to ensure continued operation of the bulk of the instruments, and to allow the batteries to recharge. This occurred at only one station, CG16-2 in the accumulation 5 zone, after consecutive periods of snowfall and persistent low cloud, resulting in four intermittent gaps in the CPEC dataset (28% of total observation time).

LiDAR
Airborne LiDAR was employed to obtain high resolution topographic data over each of the study locations, using a Riegl 580 laser scanner and dedicated Applanix PosAV 910 Inertial Measurement Unit. In general, flights were performed over 10 Nordic and Conrad glaciers twice per year (Table 3), close to the end of the winter and summer seasons (April and September), as part of an ongoing mass balance survey of the study glaciers (B. Pelto, unpublished data). By analysing the altimetry data from these times of the year, it was hoped that the variation in surface roughness due to the transition from a snow-covered to bare ice surface could be captured. In addition, the repeat mapping of each location from one year to the next would help identify the persistence in surface roughness. In 2014, April flights were not performed over the glaciers (a 15 July flight was performed over Nordic), while in 2015, the September flight over Conrad captured usable data for the accumulation zone, only.

Eddy Covariance Data
Prior to calculating observed values for the turbulent heat fluxes and roughness lengths, the raw (20 Hz) EC data were 20 passed through a series of preprocessing steps using the EddyPro data package (LI-COR, 2016). These steps are described in detail in Fitzpatrick et al. (2017), but a summary of the main techniques is provided below. A planar fit coordinate rotation method (Wilczak et al., 2001) was applied to all of the sonic anemometer data to account for misalignment of the axis of the sensor with the component of the mean air flow. For the OPEC water vapour measurements, the Webb-Pearman-Leuning correction (Webb et al., 1980) was used to correct for the density effects of air temperature fluctuations, while 25 readings from periods affected by precipitation on the analyser windows were removed. These corrections were not required for the CPEC water vapour data. The turbulence data were averaged over 30-minute blocks, and the calculated fluxes were filtered using quality tests for steady state and developed turbulent conditions, following Mauder and Foken (2004).

LiDAR Data
The trajectories of each LiDAR flight had been previously post processed using a network of permanent GPS base stations in British Columbia. The positional uncertainties of the flight trajectories were typically better than 5 cm, with the total uncertainty in the processed LiDAR point clouds better than ±10 cm, while the average point density for the LiDAR surveys over the ice-covered terrain was 1-2 laser shots per m 2 (B. Pelto, unpublished data). LAStools (Isenburg, 2006) was utilised 5 to classify the LiDAR data into ground and non-ground laser returns. The ground returns were subsequently gridded into DEMs with a 1 m 2 grid cell, the grid lines aligned with true north and east.

In Situ Roughness Length Values
Roughness length values were calculated by implementing EC data into the bulk method, with separate values calculated for OPEC and CPEC systems when both sensors were used at the same station: 10 where is the von Kármán constant (0.4), is the sensor height, and , , , * , * , and * are the 30 minute ECobserved values for mean wind speed, air temperature, specific humidity, friction velocity, and the surface layer scales for 15 temperature and specific humidity, respectively (Conway and Cullen, 2013).
� � and ℎ � � are the vertically integrated stability functions for momentum and heat (Beljaars and Holtslag, 1991;Dyer, 1974), where is the Monin-Obukhov length. Glacier surface specific humidity is calculated from atmospheric pressure , and the surface vapour pressure ( ) which is assumed to be at saturation at the glacier surface temperature ( = 0.622 ⁄ ). To minimise potential errors and to obtain roughness lengths representative of the conditions at each site, an extensive series of filters 20 were applied to the 30-minute values (see Fitzpatrick et al., 2017, for full details). These filters included a 90° wind direction window centred on the main axis of the EC sensor (to minimise the influence of flow distortion due to the station structure), minimum values for wind speed (> 3 m s -1 ) and * (> 0.1 m s -1 ), minimum differences between measurement and surface height values of air temperature (> 1°C) and vapour pressure (> 66 Pa) (Calanca, 2001;Conway and Cullen, 2013), a minimum scalar roughness length value of 1 x 10 -7 m based on the mean free path length of molecules , a 25 precipitation filter, and a test for stationarity of the turbulence (Foken, 2008). Only roughness length values calculated during near-neutral stability conditions (−0.1 < < 0.2) were retained, to minimise the uncertainty associated with the stability functions applied in Eq. (1-3) during non-neutral conditions (Smeets and van den Broeke, 2008;Conway and Cullen, 2013).

Remote Momentum Roughness Length Estimation 10
The set of 1 x 1 m grid cell DEMs obtained for the study glaciers from the LiDAR data were utilised to remotely estimate momentum roughness length values. Estimates were determined at the location of each station using the DEMs from the same year the station was in place, and compared with the EC-derived 0 _ values. September DEMs were used to estimate roughness length values for bare ice surfaces, and April DEMs for snow-covered surfaces (both the April and September DEMs at CG16-2 in the accumulation zone represent a snow-covered surface). The DEM for Nordic Glacier in July 2014 15 was used to estimate roughness lengths for the transitional snow-ice surface at NG14. The estimation of 0 was also repeated on DEMs from periods without a station present at that location to allow for an examination of the temporal variation of roughness properties at each site over the three years. Two methods were developed in this study, referred to as the (i) block and (ii) profile methods. Both methods assume that a DEM with a 1 x 1 m grid cell can adequately resolve the scale of the surface features that have the primary influence on roughness length. Where airflow encounters a dense 20 distribution of roughness elements (as can be present on an ablating glacier surface), the flow is likely to experience wakeinterference or skimming (Wieringa, 1993), reducing the relative influence of smaller scale roughness features on 0 (Smeets et al., 1999), and increasing the influence of elements that are potentially resolvable at the DEM scale.
Both methods draw on the empirical theory of Lettau (1969) for the estimation of 0 from microtopography measurements: where ℎ * is the average effective height of the roughness elements above the surface, is the average crosswind silhouette or face area of the roughness elements encountered by oncoming air flow, is the lot area, equal to the total area of the site divided by the number of roughness elements on its surface, and the value 0.5 represents an average drag coefficient. The original application of the above theory assumes that the surface is composed of regularly spaced roughness elements of similar size and shape, an assumption that may not always hold for a glacier surface. 30 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-232 Manuscript under review for journal The Cryosphere Discussion started: 12 November 2018 c Author(s) 2018. CC BY 4.0 License.

Block Estimation
The first method developed in this study to estimate 0 aimed to account for the variation in shape and distribution of roughness elements on a glacier surface. First, the form drag generated by the features on an individual portion or block of the surface was estimated, before combining the influence of each portion over a footprint to determine the momentum roughness length value for a given downwind location. Similar methods were proposed and evaluated by Kondo and 5 Yamazawa (1986) for estimating 0 over irregular surfaces. To account for the often dense distribution of roughness elements on a melting glacier surface, and the effects of this distribution on airflow, the block method developed here also considers the relative height differences and potential sheltering influence of neighbouring features on the surface.
As the method would be evaluated using roughness measurements derived from the EC systems, it was applied to subareas of each DEM that contained the potential turbulent footprint for a given station. Each subarea was 2,000 x 2,000 m in 10 dimension, and centred on the grid cell containing the station site. For each grid cell in the subarea, a one-cell-thick border was selected around the cell of interest, creating a 3 x 3 m block of cells (Fig. 2), representing a roughness element and its surrounding area of influence. A localised drag value ( _ ) was estimated for each block, by utilising Eq. (6), and building on the methods of Smith et al. (2016). The heights of the cells in the block were detrended for the mean slope of the glacier in the region of the station, as it was assumed that mean airflow was parallel to this plane. The height values within 15 the block were normalised, and the mean height of all the cells above the zero plane was assigned to ℎ * . A value for was calculated for each cardinal wind direction, as follows. The heights of the first line of cells in the block perpendicular to the oncoming wind (ℎ 1 ) set the base levels for the silhouette area, and the maximum height of the cells in each row set the upper level. The sum of the silhouette areas of each row was then assigned to the value for that block and wind direction: where is the number of rows. The area of the block was assigned to the value for . _ values were then calculated for each of the four cardinal wind directions for each grid cell; the block in Fig. 2 shifting by one cell each step: A range of border thicknesses around each grid cell, from one to five cells (3 x 3 to 11 x 11 m block area), was also implemented to test the performance sensitivity to this choice. 25 To estimate a momentum roughness length value at the location of a station, the effective influence of the _ values over the entire footprint must be determined. The flux footprint of the turbulence observed at each station was estimated using the model of Kljun et al. (2015). This model involves a two-dimensional parameterisation of a more complex, backward Lagrangian particle dispersion model (the LPDM-B model in Kljun et al., 2002). In the above study, the parameterisation was developed and evaluated for a wide range of boundary layer conditions and surface types, and was 30 shown to agree with the footprint estimates of the more complex model. To estimate the footprints for the glacier stations in this study, EC-observed values for mean wind speed and direction, 0 _ , , * , and the standard deviation of lateral wind velocity were implemented into the parameterisation. Flux footprint maps were generated from the model, with a 1 x 1 m grid cell and total area of 2,000 x 2,000 m, centred on the station location, to match the selected DEM subareas. Each grid cell was assigned a flux footprint value ( ), representing its normalised contribution to the turbulent flux observed at the station. Maps were generated for every 30-minute period in the EC data, from which an average seasonal footprint for the station was determined. For stations with two EC systems, separate footprint maps were generated for each to investigate sensitivity to the observation method. 5 The seasonal flux footprint map for a given station (or EC system) was overlaid over the corresponding _ values for the wind direction of interest. The _ value for each grid cell was then weighted by its flux footprint contribution, and summed over the subarea to obtain 0 _ : where is the number of grid cells in the subarea. This process was then repeated for the DEMs available from each season. 10 Standard error propagation methods were used to calculate the uncertainty in 0 _ by considering the uncertainties in the LiDAR height data (< ±0.1 m) and the normalised mean square error in the values from the footprint model (0.48; Kljun et al., 2015).
The primary application of a remote technique to estimate momentum roughness lengths would be to obtain values for where in situ observations are not available, and therefore, where the turbulent flux footprint for a given site is unknown. 0 _ 15 values were first calculated with EC-derived footprints, as above, to evaluate the effectiveness of the local form drag estimation (Eq. 8). To test the performance of the block method in situations when EC data is not available, the observed turbulent footprints were then replaced with a series of assumed footprint areas at each site and applied to the corresponding _ values to calculate 0 _ .

Profile Estimation 20
The second method developed in this study takes a profile-based approach to estimating momentum roughness lengths, and aims to identify the length scales relevant to form drag over that profile, rather than using the element by element approach of the previous technique. As with the block method, the first step was to detrend the surface height values for the mean slope of the glacier. Beginning with roughness estimation for the downslope (southerly) wind direction, a profile of grid cells was selected from a given DEM along the glacier slope; 600 m in length, one grid cell wide, and centred on the location of a 25 station. A linear trend was fitted to this profile to identify the slope, and the trend was then removed from the original height data ( Fig. 3a-b). This step was repeated for 50 parallel profiles on either side of the central 'station' profile (101 profiles, in total). The next step was to determine the scale of the features relevant to form drag, that is, the features that act as obstacles to air flow, and to remove large scale surface features or waves which air flow may follow rather than be impeded by. The power spectrum was calculated for the detrended profile, and a cut-off wavelength was visually identified between large and 30 small scale features. In Fig. 3c, an example of the mean power spectrum over 101 detrended profiles is shown in log-log for CG16-1 in September 2016. In this case, a cut-off wavelength of 0 = 35 m was visually identified from the plot as differentiating between wavelengths with large or small power or amplitude. With the cut-off wavelength identified, a fast Fourier transform (FFT) high pass filter was applied to the detrended profile to remove the large wavelengths (Fig. 3b), and to obtain a filtered profile. The filtering was performed in the wavenumber ( ) domain with the following steps: (i) FFT was applied to the detrended profile h(y) in Fig. 3b to get H(k); (ii) H(k) was modified by setting its values to zero for < 2π/ 0 ; (iii) an inverse FFT was applied to the modified H(k) to get the filtered profile h(y) in Fig. 3d. Finally, Lettau's method (Eq. 5 6) was applied to the filtered profile to estimate a value for momentum roughness length. was calculated as the width of the profile ( = 1 m) multiplied by the length of the fetch ( ) upwind of the station. A range of values for were applied from 0 to 2 0 in 1 m increments. The height of the grid cells along a given fetch was assigned to an array from ℎ 0 to ℎ , where is the number of grid cells in the fetch, and the standard deviation of the height array along was assigned to ℎ * . A value for was obtained from the sum of the height differences between adjoining grid cells: 10 and substituted into Eq. (6), with and ℎ * , to estimate a momentum roughness length value for a given fetch ( 0 _ ℎ ). The mean of the 0 _ ℎ values from = 0 to 2 0 was then assigned to the momentum roughness length for the station grid To examine roughness length variability in the vicinity of the station grid cell, and to determine the uncertainty in the 15 presented results, the above process was repeated for all grid cells in the 101 x 101 m area upwind of the station (i.e. 50 m either side of the station profile). The profile method was also applied over a range of angles in addition to the prevailing downslope, southerly direction, to examine the effects of changing wind direction on momentum roughness length (Fig. 4).
To do so, the x-y grid matrix of a patch of grid cells (101 m wide and 351 m long, containing the station site) was multiplied by a rotation matrix (in 5° increments between 90° and 270°). The height values from the DEM grid cells were then bi-20 linearly interpolated to the rotated grid to derive new rotated height values. A value for 0 _ was then calculated as above for profiles in line with the long axis of the patch, for each 5° increment in direction.
The sensitivity of the profile method to the use of a DEM with a finer (1 x 0.1 m) or coarser resolution (3 x 3 m) than the original 1 x 1 m DEM was tested. As a 1 x 0.1 m DEM could not be derived from the LiDAR data, a synthetic test surface was created using data from microtopography profile measurements obtained at CG16-1 at the end of the melt season. Four 25 surface height profiles, 2 m in length and with 0.1 m resolution, were obtained at distances of 10 m, 50 m, 100 m, and 150 m upwind of the station (Fig. 5a). The profiles were taken perpendicular to the prevailing wind direction (downslope), and measured using a 2 m snow probe, horizontally laid on the surface and allowed to partially melt in place. The long axis of the probe was set as the zero plane, and the height of the surface was measured relative to this level at 0.1 m spacings.
Height variability parallel to the downslope direction was expected to be smaller than in the parallel direction which 30 crosscuts supraglacial channels on the surface. Therefore, in the absence of microtopography measurements in this direction, the profile from the cross-slope direction with the smallest variance i.e. the 10 m upwind profile (Fig. 5b), was used to represent the slope-parallel variance. This 2 m profile was demeaned at a 1 m interval and lined up in a repeated sequence to obtain an extended (600 m long) synthetic microtopography profile. The final test profile was constructed by adding this extended synthetic profile to the detrended profile in the downslope wind direction from the 1 x 1 m DEM. The same synthetic profile was added to the detrended profiles from each side of the station, at 1 m distance apart, yielding the synthetic 1 x 0.1 m DEM. The 3 x 3 m DEM was created by applying a 2-D smoothing of the original 1 x 1 m DEM, using a 3-point running mean in both x (Easting) and y (Northing) directions. The profile method was then applied to both the 1 x 0.1 5 m and 3 x 3 m DEMs for the 600 x 101 m area upwind (slope-parallel) of the station, using the same steps as outlined previously. The same threshold wavelength, 0 = 35 m, was used to filter the profiles. Figure 5c displays examples of filtered profiles, h(y), as derived from the three DEM resolutions.

EC-Derived Roughness Lengths 10
The geometric means of the roughness length values calculated from each EC dataset are presented in Table 4, with separate 0 _ values for periods with snow and ice surfaces. Each of the observed 30-minute roughness length datasets were found not to have a normal distribution (using one-sample Kolmogorov-Smirnov tests), but one that was approximately lognormal. For presenting mean EC-derived values in the remainder of this study, geometric means are used to avoid excessively weighting the larger roughness values (Andreas et al., 2010). Stable atmospheric conditions persisted over the 15 glaciers for much of each season, limiting the number of suitable 30-minute periods for roughness calculation after application of the filters discussed in Sect. 2.5 (number of available measurements presented in Table 4). Across all test sites, 0 _ had a mean of 2.3 mm and 1 mm for ice and snow, respectively, while the scalar roughness lengths had mean values of 0.05 mm for 0 _ and 0.11 mm for 0 _ . Where OPEC and CPEC systems were used on the same station, the OPEC system returned slightly larger mean 0 _ values (2.8 mm and 1.4 mm, respectively). Mann-Whitney U tests applied to the 20 30-minute roughness values from CG15-1 rejected the null hypothesis that the 0 _ values from the OPEC and CPEC systems had the same distribution (p < 0.01), but the hypothesis could not be rejected for the scalar values (p > 0.5).
The ice 0 _ values were within the expected range for moderately rough glacier ice (Brock et al., 2006). Where measurements were repeated in the same area a year apart (CPEC observations on CG15-2 and CG16-1), persistence in the mean ice roughness length values was noted (0.86±7.4 mm, and 0.74±6.4 mm, respectively), with a failure to reject the 25 hypothesis of equal distributions (p = 0.16). Within a season, substantial variability was noted in the 30-minute 0 _ values for each ice surface (Fig. 6a), but with no evident trend in 0 _ due to changes in surface roughness over time. Mean momentum roughness lengths for snow were also within previously observed values on glacier surfaces, with a particularly large mean value observed at CG16-2 in the accumulation zone (2.4±16 mm). Extensive variability was also present in the 30-minute 0 _ values for CG16-2 (Fig. 6b), with a general increasing trend in roughness over the season. Across all 30 stations and seasons, substantial variability was noted in the mean scalar roughness lengths, with 0 _ , in particular, showing a range of two and a half orders of magnitude. 0 _ exhibited less variability (~one order of magnitude), with similar mean values observed for CG15-2 and CG16-1 (0.03±0.28 mm and 0.05±0.29 mm), and a failure to reject the null hypothesis of equal distributions (p = 0.11).
The ratios of the 30-minute EC-determined scalar roughness lengths to 0 _ were expressed as a function of * using the data from all stations and seasons (Fig. 7). These values were compared with the surface renewal models of Andreas (1987) 5 andvan den Broeke (2008). The seasonal mean ratios and * were also compared with these models. In general, the roughness ratios were shown to decrease with increasing * , with substantial scatter in the 30-minute values. The

Block Method
_ maps were generated from LiDAR-derived DEMs using the block estimation method (Fig. 8a-b) for all available years and seasons, and for each of the four cardinal wind directions. Substantial variation in _ was observed across each glacier surface, ranging from 10 -4 m for snow-covered grid cells to 10 -0.5 m for large crevasses. Figure 9 displays the seasonal turbulent flux footprint maps generated using the model of Kljun et al. (2015) for each EC sensor deployment. In 15 general, the fluxes were sourced from regions to the south of each station, in line with the prevailing downslope winds at each site. Over 80% of flux contribution came from an area within 200 m upwind of each station, with concentrated peak source regions 15-20 m upwind, on average. The flux footprints of each EC dataset were merged with the corresponding _ maps (Fig. 8c), producing a series of 0 _ values for each site. As stated, wind direction was predominately from the south during each station deployment, so the roughness estimates for this wind direction (Table 5)   As previously stated, the sensitivity of roughness length estimation to the selected block size was tested by varying the border thickness around the grid cell of interest. Overall, increasing the block area was found to lead to an increase in estimated roughness length for a given footprint, with a border thickness of 1 cell (3 x 3 m block area) returning roughness lengths closest to the EC-derived values at all stations (e.g. CG16-1 ice 0 _ = 1.6, 1.9, 2.0, 2.2, and 2.4 mm for an increasing border thickness range of 1-5 cells). 5

Profile Method
The detrending and filtering of the surface height data, as shown in Fig. 3, was performed for downslope profiles at each station site using the DEMs for all available years and seasons. The same approximate value for the cut-off wavelength ( 0 ≈ 35 m) was identified at each station site. 0 _ values were then estimated for each station location, and for each grid cell in a 101 x 101 m upwind area (Fig. 10a), from all corresponding DEMs. Table 5 presents the 0 _ values for each station 10 and LiDAR flight. Mean 0 _ values for ice and snow surfaces, over all sites and seasons, were 4.3 mm and 1.1 mm, respectively. Where repeated over the same location, the 0 _ values displayed substantial differences from one year to the next over ice surfaces (up to 5 mm), in contrast to the noted 0 _ persistence. was reduced by a factor 10 (from dm to cm scale) and 0 _ recalculated. The resulting roughness length values were reduced and matched more closely the original 0 _ from the 1 x 1 m DEM, however, still yielding up to 10% larger values than original (Fig. 10b).

In Situ vs Remote Methods
The estimates from both DEM-based roughness methods were compared with the EC-derived values ( Fig. 11 and Table 6). 30 In cases where LiDAR data were not available from the same year a station was in place, the averages of the roughness estimates from the two other years were utilised for the comparison. Overall, estimates from both DEM-based roughness methods provided values for ice and snow surfaces in line with previous observations on glacier surfaces (Brock et al., 2010), and were generally within 2 mm and 0.2mm (< half order of magnitude) of the corresponding 0 _ observations over ice and snow, respectively. Over ice surfaces, the 0 _ values were slightly smaller than the corresponding 0 _ values (mean values of 3.1 mm and 4.3 mm, respectively), and tended to align more closely with the 0 _ estimates (mean 5 of 2.3 mm). For the snow surface at CG16-2 in the accumulation zone, the mean roughness lengths from both DEM methods (0.4 mm) substantially underestimated the 0 _ value (2.4 mm). Potential causes for this deviation will be discussed in Sect. 4.1.2. For the transitional snow/ice surface present at NG14 during the first four days of observations, the 0 _ and 0 _ values from the July 2014 flight (4.5 mm and 6.8 mm) aligned more closely with the mean 0 _ value for ice over the season (4.5±28.8 mm) than with the 0 _ value obtained during the four day period (0.5±3.0 mm). The mean 0 _ 10 value for this period, however, was based on a very limited number of EC observations after filtering (n=16) with substantial scatter.

Wind Direction and Momentum Roughness Length
The 30-minute EC data and the rotated 0 _ values were used to examine the influence of wind direction on the effective roughness length at each location. It should be restated at this point that the 0 _ values had been filtered to remove values 15 when wind direction was beyond ±45° of the main axis of the EC sensor to minimise the influence of flow distortion due to the station structure. Therefore, only a limited direction window is available in the 0 _ data over which to examine this dependence. For the ice surface of CG16-1, 0 _ values were observed to increase and become more scattered as the wind direction veered towards the southwest, a pattern that was also detected in the 0 _ (Fig. 12a). Similar behaviour was noted at the same location in 2015 (CG15-2), with greater variation in 0 _ with wind direction (Fig. 12b). 20 The rotated 0 _ values were also used to examine a wider angle of wind direction than was possible with the 0 _ data.

Ice Surfaces
Variation in both the 0 _ and DEM-based roughness length values was noted across test sites with a melting glacier ice surface (e.g. 4.5 mm and 0.7 mm for mean 0 _ at NG14 and CG16-1, respectively). An assumed 0 value for ice (e.g. 5 1mm), applied uniformly to all locations in this study, would have substantially misrepresented the surface roughness characteristics, and the resulting turbulent flux parameterisations. In the case of NG14, implementing the commonly assumed 0 value for ice of 1 mm in the bulk parameterisation of turbulent heat fluxes, rather than the mean observed value of 4.5 mm, would result in a ~20% reduction in the mean estimated fluxes. Furthermore, stations throughout the study were installed in secure regions of the glaciers with relatively smooth and uniform surfaces, and away from crevasse fields and 10 glacier margins where the surface drag on airflow would be higher (Fig. 8). Therefore, the true range of roughness length values over the entire surface of the study glaciers would be greater than that represented by the values estimated for the Over the study period, the mean momentum roughness length estimates for ice at each site showed little temporal variance 15 from one year to the next. This persistence in seasonal ice roughness values may allow for the use of 0 estimates from preexisting EC or DEM campaigns at a site of interest. The period of validity of these estimates may vary, however, depending on the surfaces processes of each glacier. Within a single melt season, there was substantial scatter observed in the 30minute 0 _ values (Fig. 6a). Changes in momentum roughness length due to the evolution of the ice surface through the season were not evident in the 0 _ values, however. Previous glacier roughness studies (e.g. Sicart et al., 2014) have also 20 noted persistence in 0 despite extensive ice melt. Smith et al. (2016) noted that this persistence was most evident over ice surfaces with defined melt features, such as supraglacial channels, similar to the ice surfaces of this study. While estimated using EC-observed data, the 0 _ calculations are still derived from the bulk aerodynamic method (Eq. 1). Extensive filtering was applied to 0 _ values, in particular, to avoid uncertainty in the bulk method due to non-neutral stability conditions. However, using a filter that allows values from near-neutral conditions (−0.1 < < 0.2) rather than strictly 25 neutral, only ( = 0), may introduce some uncertainty and variability to the 0 _ estimates. Furthermore, previous studies have suggested that some assumptions of the bulk method, namely, constant momentum and heat flux values with height, may not be valid during katabatic conditions with shallow wind maximums, which can develop frequently over sloped glaciers (e.g. Denby and Smeets, 2000).

Snow Surfaces
Large differences in 0 between sites were also noted in this study for snow-covered surfaces. The annual persistence in roughness values observed over ice was also present in the snow surface values, with similar 0 _ values returned for the same time each year when repeated at the same location. Where both in situ and remote values over snow surfaces were available, agreement between 0 _ and the DEM-obtained roughness values varied substantially. In the case of CG16-2, 5 which had a snow-covered surface throughout, the relatively large mean 0 _ value (2.4±16 mm) was substantially greater than 0 _ and 0 _ (both 0.4 mm). This difference may be due to the temporal variance in roughness of a snow surface within a melt season (as observed in Fig. 6b), and the difference in observation time between the EC and LiDAR data. Images from the time lapse camera installed at CG16-2 ( Fig. 14a-b) illustrate the variety in roughness conditions of the snow surface at that site. Two periods were selected with visually apparent roughness differences and an adequate number of 10 30-minute 0 _ observations; a moderately smooth, melting snow surface (June 30 th -July 3 rd ; 78 observations), and a rough, sun-cupped surface (Aug. 19 th -21 st ; 38 observations). Examining the 0 _ values, an order of magnitude difference was noted between the mean values for the moderately smooth (1.0±4.2 mm) and rough (9.6±21.7 mm) snow surfaces. In view of this short-term variability in snow roughness, the 0 _ and 0 _ values, derived from LiDAR flights in April and September, cannot be considered comparable to the 0 _ values from the summer. Imagery taken in the same location 15 as the station site a few days after the April LiDAR flight (Fig. 14c) show a very smooth snow surface. With fresh snowfall in late August and September, a similar surface was likely present during the second flight, resulting in the small DEMbased values returned. Relatively large 0 _ and 0 _ values were obtained for NG14 during the April LiDAR flights, possibly in response to a rough snow surface. Comparable in situ imagery of the site was not available for these periods, however. The effect of the size of the roughness elements on a melting snow surface is discussed further in Sect. 4.2. 20

Wind Direction
Evidence of roughness length dependence on wind direction was observed in the 30-minute EC data at some locations, and in the rotated 0 _ values, also. The strongest dependence on wind direction in the 0 _ values was noted for the ablation zone of Conrad Glacier, at the location of CG15-2 and CG16-1. Elongated roughness features, including meltwater channels, were present on the surface during these observations, with the orientation of their long axes pointing in a 25 southeast to northwest direction (Fig. 12c). As the wind veered to the southwest, airflow became perpendicular to the faces of these features, likely resulting in increased form drag, which produced the larger roughness lengths observed. The rotated 0 _ values for the three stations in Conrad's ablation zone revealed an increase as wind direction approached a crossglacier orientation. At NG14, the pronounced increase in roughness over the ice surface at 240° was likely due to a crevasse field to the west of the station. This feature was not evident in the April values, suggesting snow cover had smoothed the 30 surface in that region. Dependence of momentum roughness length on wind direction has been observed in several other glacier studies (e.g. Munro, 1989;Brock et al., 2006;Smith, 20014). Over all seasons and locations in this study, wind direction was found to be within 45° of the mean slope angle for approximately 93% of the time. This persistent, katabatic downslope wind is a common feature in glacial boundary layers, and as a result, will substantially reduce the influence of surface roughness anisotropy on the variation in the effective roughness lengths and mean generated turbulence.

Performance of DEM-based z0v Estimation
The methods developed here for remotely estimating 0 were found to returning roughness length values within 1-2 mm 5 (<< an order of magnitude) of those determined from in situ EC measurements, and were shown to respond to changes in surface cover from snow to ice. Using a DEM with a 1 x 1 m grid cell appears to resolve the length scales influencing 0 on the ice surfaces of this study. With a dense distribution of roughness elements (Fig. 12c), the previously mentioned effects of wake-interference and skimming of the airflow over the ablating ice may have reduced the influence of the smaller roughness elements on 0 , as noted in previous studies (e.g. Wieringa, 1993;Smeets et al., 1999). During the April flights 10 over Conrad Glacier, the DEM methods returned roughness values in line with previous observations over smooth, fresh or compacted snow surfaces (e.g. Brock et al., 2006). Over rough, undulating snow surfaces, larger-scale features will have the dominant influence on roughness length (Fassnacht et al., 2009), and are potentially resolvable in the utilised DEM, as may have been the case with the April values for NG14. Over smoother surfaces, however, it is likely that the roughness elements influencing 0 are not resolvable with a 1 x 1 m DEM, making the usefulness of these methods over a melting snow surface 15 uncertain (in addition to the temporal variation discussed in Sect. 4.1.2).
The profile method developed here has been shown to return values in line with in situ estimates of the momentum roughness lengths without the need for the assumptions employed by the block method. The estimation of a similar cut-off wavelength at each station ( 0 = 35 m) indicates a common length scale, above which, the surface features are no longer an influence on roughness length. The 0 _ values did show a tendency towards overestimation, relative to the 0 _ values. 20 In addition, the persistence between seasons in roughness length, noted in the 0 _ and 0 _ values, was less evident in the 0 _ values, suggesting that the profile method is sensitive to changes in small scale features which may not have a substantial influence on the observed ( 0 _ ) roughness values. The profile method also displayed sensitivity to the choice of DEM resolution, arising from substantial differences in the estimate of (Eq. 10) for different resolutions (>50% difference between 1 x 1 m and 1 x 0.1 m resolutions). 25 The block estimation method returned roughness length values that were smaller than those from the profile method, and more in line with mean 0 _ , in general. The technique used in the 0 _ method to calculate across overlapping block areas (as shown in Fig. 2 and Eq. 7) was developed in an effort to account for the shadowing of elements from airflow by upwind features. Rather than assuming that each feature above the mean surface has an additive influence on roughness length, as done in the 0 _ method (below the cut-off wavelength) and other profile-based methods (e.g. Munro, 1989;30 Arnold and Rees, 2003), the relative height differences and potential sheltering influence of neighbouring features in the block are considered. On glacier surfaces, where elongated roughness features such as melt channels are common, the block approach may also help account for the channelisation of air flow and the shadowing of the roughness element by the upwind continuation of the feature, which in turn, may reduce the effective roughness length. The response of the block method to this effect can be seen when the _ estimates for the southerly (downslope) wind direction are compared with those for the westerly (cross-slope) wind direction (Fig. 15). Drag values estimated for the meltwater channels on the surface are lower when air flow is close to parallel to these features, and higher when air flow is perpendicular to the channels. This 5 effect may have led to the smaller 0 _ values, relative to the 0 _ values. When implemented with an assumed turbulent footprint, the block method returned roughness length values in line with those calculated using a footprint model or from EC data (Table 5), indicating the potential for its use where turbulence observations are unavailable.
To apply the block approach, a number of additional assumptions were required, however. The choice of block size corresponds to an assumption on the size of the dominant roughness elements influencing 0 on the glacier surface, and an 10 assumption on the range of a feature's shadowing effect. The downwind shadowing generated by a feature will likely vary with wind speed, and this variation is not accounted for here. The optimal block size may vary between locations and wind regimes, and require tuning for application to other surfaces. Over the range of surfaces in this study, however, a 3 x 3 m block (applied to a 1 x 1 m DEM) was shown to be optimal, and to respond to changes in surface roughness due to snow and ice cover. As a test of robustness, 0 _ values were also estimated for a region of forest captured in the LiDAR data. This 15 forest was located on a valley floor to the east of Conrad Glacier, and consisted of tall (~20 m), coniferous trees. The 0 _ value for a 200 x 200 m subarea within this forest was 1.28 m. This value is in line with existing 0 measurements over coniferous forest (Wieringa, 1993). While the 0 _ method (including the LiDAR data utilised) is not configured nor intended for use over forestry, this test indicates that its configuration (including selected block size) is responsive to a wide range of roughness element sizes, beyond the scale of those encountered on the glacial surfaces of this study. 20

Scalar Roughness Relationships
Whilst displaying similar mean values over the entire dataset (0.05 mm for 0 _ and 0.11 mm for 0 _ ), the scalar roughness lengths differed substantially from each other when examined on a site-by-site basis. There was no evidence of a consistent ratio between 0 _ and 0 _ , with their seasonal means ranging above and below each other by up to an order of magnitude. Between the momentum and scalar roughness lengths, seasonal 0 _ displayed a more consistent relationship 25 with 0 _ , being approximately one and a half orders of magnitude smaller than 0 _ in most cases. This relation did not hold for NG14 and CG16-2, however, and between 0 _ and 0 _ , there was no persistent ratio. Calanca (2001) observed 0 to be a function of the temperature gradient between the air and a melting ice surface, while Park et al. (2010) found a relation between relative humidity at 2 m height and 0 . In this study, variation in the scalar roughness lengths was compared with fluctuations in air temperature gradient and relative humidity, but no dependent relationship was evident. The 30 surface renewal model of Andreas (1987), where the ratio of momentum to scalar roughness was expressed as a function of * , showed relatively good performance, particularly for seasonal values of 0 . If momentum roughness length values have been obtained for a given surface (through remote or in situ methods), this model appears to be the best available method for estimating the scalar values.

Conclusions
Over three melt seasons, in situ and remote methods were implemented to determine the momentum and scalar roughness lengths on the surface of two glaciers in the Purcell Mountains of British Columbia, Canada. EC sensors were employed to 5 obtain continuous in situ measurements throughout each melt season, while LiDAR-derived DEMs were utilised in the development of two remote estimation techniques. Seasonal mean momentum roughness length values, estimated from eddy covariance observations at each location, ranged from 0.7-4.5 mm for ice surfaces, and 0.5-2.4 mm for snow surfaces. For representative turbulent flux modelling, this study suggests that site-specific 0 values are necessary, particularly in the case of distributed glacier models. From year-to-year, 0 values were noted to remain relatively consistent at a given location 10 (<0.2 mm difference between seasonal mean values). Within a melt season, continuous EC observations and camera imagery noted greater temporal variation in roughness for snow surfaces than for ice. These findings indicate that site-specific 0 values on an ice surface may be valid to implement over multiple melt seasons, while over snow surfaces, the utilised roughness values require intraseasonal updating. Wind direction was also noted to affect 0 variability where elongated features such as melt channels dominated the surface topography. Persistence in wind direction on sloped glacier surfaces, 15 however, reduces the influence of this variability.
Observations of the scalar roughness lengths differed substantially from the corresponding momentum values, showing considerable variation between location and season, and little agreement with fixed ratios commonly assumed with 0 . In general, the Andreas (1987) surface renewal method showed agreement with the observed ratios between EC-derived scalar and momentum roughness lengths, and would seem to be the appropriate method to implement where continuous EC 20 observations are not available, but site-specific 0 values have been established.
The DEM-based methods described in this study were shown to perform well over most surfaces, differentiating between ice and snow cover, and returning momentum roughness values that were within 1-2 mm (<< an order of magnitude) of ECderived values for the corresponding footprints. The features dominating glacier surface roughness, particularly over ice, were found to be of a scale that was resolvable using a 1 x 1 m DEM, and persistent enough for a DEM-based roughness 25 estimate to be usable over an extend period of time. This may allow for the potential upscaling of these methods with high resolution satellite imagery, greatly expanding the number of glaciers for which roughness length estimates could be obtained. Over melting snow surfaces, the validity time of a retrieved DEM is reduced due to the discussed temporal variability in roughness, and as a result, the estimated roughness lengths may quickly become unrepresentative. In addition, the roughness features observed to develop on melting snow in this study may not be resolvable using a 1 x 1 m DEM, and 30 further testing over snow, with simultaneous in situ and remote observations, would be useful. Author contributions. N. Fitzpatrick prepared and executed the on-glacier field campaign, conducted the data analysis, and proposed the roughness estimation methods presented in the chapter, as well as writing the text. V. Radić helped develop these methods and the associated sensitivity testing, supervised the analysis, and provided financial and field support. B.
Menounos instigated the LiDAR imaging flights, and provided the resulting digital elevation models to the study, in addition to logistical and minor financial support to the field campaign. 5 Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Funding support for this study was provided through the Natural Sciences and Engineering Research Council (NSERC) of Canada (Discovery grants to V. Radic and B. Menounos). The eddy covariance system and 10 meteorological equipment were supported by NSERC Research Tools and Instruments grant and Canada Foundation for Innovation grant (to V. Radic). The authors wish to thank Ben Pelto for his assistance in the field and with the provision of the DEMs, Steve Conger and Tannis Dakin for logistical support, Zoran Nesic for assistance with station design and sensor calibration, and Jörn Unger for his work on the quadpod construction.            Table 4) extends beyond the y-axis range.         and  Table 6. Comparison of momentum roughness length values (in mm) for each station, as observed from the EC systems ( _ ), and as estimated using the DEM-based methods ( _ and _ ). *For years where LiDAR data was not available from the same year a station was in place, the averages of the roughness estimates from the two other years were utilised for evaluation.