Articles | Volume 15, issue 12
Research article
09 Dec 2021
Research article |  | 09 Dec 2021

Aerodynamic roughness length of crevassed tidewater glaciers from UAV mapping

Armin Dachauer, Richard Hann, and Andrew J. Hodson

The aerodynamic roughness length (z0) is an important parameter in the bulk approach for calculating turbulent fluxes and their contribution to ice melt. However, z0 estimates for heavily crevassed tidewater glaciers are rare or only generalised. This study used uncrewed aerial vehicles (UAVs) to map inaccessible tidewater glacier front areas. The high-resolution images were utilised in a structure-from-motion photogrammetry approach to build digital elevation models (DEMs). These DEMs were applied to five models (split across transect and raster methods) to estimate z0 values of the mapped area. The results point out that the range of z0 values across a crevassed glacier is large, by up to 3 orders of magnitude. The division of the mapped area into sub-grids (50 m× 50 m), each producing one z0 value, accounts for the high spatial variability in z0 across the glacier. The z0 estimates from the transect method are in general greater (up to 1 order of magnitude) than the raster method estimates. Furthermore, wind direction (values parallel to the ice flow direction are greater than perpendicular values) and the chosen sub-grid size turned out to have a large impact on the z0 values, again presenting a range of up to 1 order of magnitude each. On average, z0 values between 0.08 and 0.88 m for a down-glacier wind direction were found. The UAV approach proved to be an ideal tool to provide distributed z0 estimates of crevassed glaciers, which can be incorporated by models to improve the prediction of turbulent heat fluxes and ice melt rates.

1 Introduction

The aerodynamic roughness of a glacier influences the turbulent heat exchange between the glacier surface and the atmosphere (Rees and Arnold2006). Both sensible and latent heat fluxes lead to this heat exchange on the surface and therefore have a large impact on the meltwater production and the surface energy balance of glaciers (Hock2005). The bulk approach for the calculation of those turbulent fluxes is very popular due to its low requirements for data collection. It only requires basic atmospheric measurements (e.g. wind speed, temperature), as well as the aerodynamic roughness length of the surface (Chambers et al.2020). The aerodynamic roughness length, also called z0, is a length scale that represents the height above the surface at which the wind speed drops to zero (Chappell and Heritage2007). It is often considered a surface characteristic and describes the loss of wind momentum that can be attributed to surface roughness (Smith2014).

Fitzpatrick et al. (2019)Obleitner (2000)Smeets et al. (1999)Smith et al. (2016)Smith et al. (2016)

Table 1Published aerodynamic roughness length values for crevassed glacier areas.

Download Print Version | Download XLSX

Recently, a series of studies (e.g. Irvine-Fynn et al.2014; Miles et al.2017; Smith et al.2016) have determined z0 values of glacier surfaces based on the bulk approach while using digital elevation models (DEMs) as the basis for calculations. Terrestrial light detection and ranging (lidar) systems (used for instance in the studies of Smith et al.2016, Nicholson et al.2016, and Nield et al.2013a) constitute a powerful way to effectively produce such DEMs. However, they are very expensive (Uysal et al.2015) and limited in the area they cover (Irvine-Fynn et al.2014). Thus, uncrewed aerial vehicles (UAVs), often also called unmanned aerial vehicles or drones, provide a cheap alternative to overcome these limitations (Uysal et al.2015) because they are more flexible in their use and less limited by local topography as they provide a bird's-eye perspective. In recent years, UAVs have presented new opportunities for detailed mapping of the earth surface and have become more and more popular in the field of glaciology (Bhardwaj et al.2016) and in Svalbard (Hann et al.2021). The main advantage of UAVs is the possibility of collecting high-temporal- and high-spatial-resolution data at low costs (Casella and Franzini2016) and to overcome the gap between sparse field observations and coarse-resolution space-borne remote sensing data (Bhardwaj et al.2016).

Several studies already investigated the z0 values of non-crevassed (Irvine-Fynn et al.2014; Nield et al.2013a), debris-covered (Quincey et al.2017), sparsely crevassed (Smith et al.2016) or rough (van Tiggelen et al.2021) glacier ice surfaces using different approaches. However, still little is known about the effect of heavily crevassed glacier surfaces on the turbulent heat exchange between glaciers and the atmosphere. A broader-scale, heterogeneous surface topography of obstacles (i.e. crevasses) makes the definition of z0 values challenging (Quincey et al.2017). Typical values of z0 on glacier ice range from less than 0.0001 m for smooth ice to 0.02–0.1 m for rough glacier ice (Brock et al.2006; Smeets et al.1999; van Tiggelen et al.2021), but published z0 values for large and deep crevasses are rare. Table 1 gives a closer overview of such published aerodynamic roughness values and shows z0 values up to 0.5 m for large crevasses (Fitzpatrick et al.2019). Furthermore, most surface energy balance models only consider one single z0 value for the whole glacier (commonly 0.001 m; Smith et al.2020), regardless of any spatial and temporal variability (Quincey et al.2017). The aerodynamic roughness length is a key parameter for the calculation of turbulent fluxes (Chambers et al.2020) since a change in z0 by an order of magnitude can double the estimated turbulent fluxes (Brock et al.2006; Munro1989). The uncertainty in z0 values therefore presents a serious challenge for the calculation of surface ice melt (Smith et al.2016), and its accurate parameterisation is crucial, especially for complex ice surfaces.

The objective of this study is to assess the application of UAVs for capturing spatially variable z0 values of heavily crevassed ice surfaces. A photogrammetry method is used to build DEMs of the mapped glaciers from the aerial images, which are then utilised to estimate the aerodynamic roughness length of crevassed tidewater glaciers in Svalbard. The main advantage of this approach is that increasingly available UAV technology can be used to estimate turbulent fluxes that are usually difficult to measure in the field. Furthermore, the chosen DEM approach allows glaciers to be divided into areas of different aerodynamic roughness length values, leading to a better spatial representation of the turbulent fluxes and therefore surface ice melt on glaciers.

2 Data and methods

The following section describes how a DEM was generated from aerial imagery of crevassed glaciers in Svalbard. Images were obtained using off-the-shelf UAVs. In addition, different methods to calculate the aerodynamic roughness length from the DEMs are introduced.

Figure 1Field sites: map of Svalbard with marked locations of the investigated tidewater glaciers (red dots). Additionally, Sentinel-2 satellite images (ESA2020) of the glaciers Nordenskiöldbreen (a), Tunabreen (b), Fridtjovbreen (c) and Heuglinbreen (d) taken in the according week of fieldwork provide a closer look. The lines mark the mapped front area for each glacier in 2019 (blue) and 2020 (red).

Table 2Overview of the tidewater glaciers visited and the date of their UAV survey. Additionally, the size of the mapped area and DEM resolution, as well as the average height and length of roughness elements for a up-/down-glacier (cross-glacier) wind direction are listed.

Download Print Version | Download XLSX

2.1 Field sites

Four heavily crevassed tidewater glacier termini in central Spitzbergen (see Fig. 1 and Table 2) were visited during three field campaigns. Nordenskiöldbreen was visited in summer 2019 and 2020. In spring 2020 further fieldwork was conducted on Fridtjovbreen, Heuglinbreen and Tunabreen. Fridtjovbreen is a single tidewater glacier of about 13 km length, flowing southwards and terminating in Van Mijenfjorden on the western side of Spitzbergen island. Here precipitation and temperature are relatively high with an annual temperature of about 4 C and precipitation of up to 1000 mm (Hanssen-Bauer et al.2019). Both Nordenskiöldbreen and Tunabreen are outlet glaciers (flowing from northeast to southwest) draining the large Lomonosovfonna ice cap, where the precipitation usually is lower than on the west coast (Hagen et al.1993). While Nordenskiöldbreen is a roughly 15 km long, wide tidewater glacier terminating in Billefjorden, Tunabreen is narrower, with a length of about 20 km and terminating in Tempelfjorden. Additionally, both Tunabreen and Fridtjovbreen are known to have experienced a surge event. While the surge on Fridtjovbreen already happened during the 1990s (Murray et al.2012), Tunabreen surged more recently in the years 2003 to 2005 and had another advance of the glacier front about 10 years later (Ericson et al.2019). Heuglinbreen is a tidewater glacier flowing southwards and terminating into Mohnbukta, a bay on the east coast of Spitzbergen. This region is known to be particularly cold (annual temperature of about 10 C) and dry (annual precipitation up to 700 mm) (Hanssen-Bauer et al.2019).

Figure 2The originally mapped DEM (a) of Fridtjovbreen was rotated and cropped (black frame) before being used for the model calculations. Additionally, the elevation values of the sea ice area were removed, and the whole DEM was divided into sub-grids of 50 m× 50 m (b).


2.2 Field data collection

UAV-based aerial imagery was collected during each fieldwork campaign with off-the-shelf UAVs (a DJI Mavic 2 Enterprise and a DJI Phantom 4 Pro). In order to have a sufficient resolution of the crevasse fields, a target DEM resolution of about 0.25 m px−1 was chosen. To achieve this DEM resolution, the UAV imagery aimed for a ground sampling distance (GSD) of at least 0.1 m px−1, a forward overlap of 90 % and a side overlap of 80 %. During fieldwork the UAVs were planned to operate at an altitude of 200 ma.g.l. (above ground level), taking nadir-viewing pictures. Flights were conducted with pre-programmed waypoints and a run separation of 70 m. Since information about the glacier surface elevation was unknown to the UAV, the true altitude above the glacier surface was typically less than 200 m. As a result, the GSD ranged between 0.04 and 0.07 m px−1, and the subsequent DEM resolution ranged between 0.17 and 0.28 m px−1 (see Table 2).

2.3 DEM generation and preparation

The high-resolution images from the UAVs were processed with a structure-from-motion (SfM) multi-view stereo (MVS) photogrammetry method using Agisoft Metashape version 1.6.2 (Agisoft LLC2020). Building DEMs in Agisoft is a three-step process: image alignment, construction of a dense cloud and DEM generation (Verhoeven2011). The software runs on a fully automatic workflow. However, the whole processing comes along with many parameter settings that can be selected to improve the DEM quality depending on the input data and the output purpose. To determine an ideal set of parameters for our approach, many different combinations of settings were tested, leading to the optimised final settings (see Dachauer2020). Due to the inaccessibility of the glacier surface, no ground control points (GCPs) could be placed on the mapped area for georeferencing. No alternative georeferencing platforms such as real-time kinematic correction (Chudley et al.2019) were available. While this is a recommended procedure for future applications of this technique, we point out that computation of z0 requires quantification of relative topographic differences, and so the impact of this shortcoming is minor.

Before starting the model calculations of the aerodynamic roughness length, the DEMs obtained from the SfM processing were rotated in such a way that the glacier flow direction corresponds to the column alignment of the DEM, with the front at the lower part (see Fig. 2). This allowed for the estimation of z0 values for the following four wind directions: down-glacier, up-glacier and cross-glacier from both sides. Furthermore, the DEMs were cropped to a shape that contains the most crevassed zones. This area is specific to each glacier and varies significantly due to the differences in glacier size, time available for mapping, number of UAV batteries and limitations due to GPS interference. The sea ice and water area in front of the glaciers was removed since it makes no contribution to the turbulent heat exchange of the glacier.

Figure 3Graph (a) shows a random transect of two roughness elements on Fridtjovbreen (blue line). The transect is linearly detrended either over the whole profile (green) or in 10 m intervals (purple). Graph (b) illustrates the linearly detrended surface data applied at the whole transect (green line) or applied at 10 m intervals (purple line).


This study assumes that the mean airflow is blowing parallel to the slope of the glacier. This means that the glacier slope has no effect on the aerodynamic roughness (Fitzpatrick et al.2019). Furthermore, it justifies the assumption that the aerodynamic roughness is in the first place influenced by the macro-structure of the surface (crevasses, large obstacles, etc.) and not only by small-scale surface roughness on the crevasse obstacles. In other words, only looking at the small-scale surface roughness elements would lead to a wrong roughness parameterisation since they might be located on the inner side of a large-scale roughness obstacle not exposed to the whole mean airflow. Accordingly, the chosen grid size must be large enough to include the macro-structure of the surface because small-scale roughness elements alone do not represent the real topographic expression. Linear detrending over long baselines manages to represent areas of high curvature (Smith2014) and is therefore appropriate for this purpose. Figure 3 presents the impact and importance of the chosen sub-grid size (i.e. length of detrended transect) on the modelled surface roughness. The blue line in Fig. 3a shows a random transect of two roughness elements of 30 m width on Fridtjovbreen. Two linear detrending methods were applied to the surface data of this illustration. First, the green line detrended the whole transect. Second, the transect was detrended in 10 m intervals (purple line). The detrended data (Fig. 3b) show that the 10 m grid size only captured the small-scale surface roughness of the obstacles (purple line). In contrast, the linearly detrended transect length of 60 m (green line) managed to represent the two large roughness elements. Therefore, a sub-grid size of at least the width of an average obstacle should be chosen to account for the macro-structure surface roughness of the crevassed glaciers. On average, the mapped tidewater glaciers had obstacle widths between 30 and 50 m (calculated according to the transect method of Munro (1989), see Table 2). Thus, the final DEMs were subdivided into rectangular sub-grids of 50 m× 50 m (grid in Fig. 2), each estimating one z0 value.

2.4 Models for aerodynamic roughness length estimation

The aerodynamic roughness length was calculated for each DEM sub-grid and all four wind directions with five models after 2D linear detrending (see Fig. 3). The most common models in glaciology for calculating z0 using the bulk method are based on the work of Lettau (1969), who developed the following equation for the bulk aerodynamic roughness length:

(1) z 0 = c d h s S ,

where h is the effective obstacle height (m), s is the silhouette or frontal area of the obstacle (m2), and S is the horizontal ground area (m2). The value cd=0.5, first proposed by Lettau (1969), corresponds to the average drag coefficient of a characteristic roughness element.

Table 3Overview of parameters from Eq. (1) of Lettau (1969) used for the z0 calculation of each model.

Download Print Version | Download XLSX

The definition of the parameters of Eq. (1) is more complex in glacial environments since the individual roughness elements are non-uniform and vary substantially in height, size and density (Chambers et al.2020). Therefore, the original equation has been adjusted and further developed in several studies, leading to the five different models used in this study (Table 3). The five models can be subdivided into two groups according to how they determine and measure roughness elements. One group counts the number of roughness elements in a transect (hereafter called the “transect method”), while the other group is based on a raster approach (hereafter called the “raster method”) using every DEM cell value for the calculation of z0. The choice of the models is justified by comparing results of the two methods and to give further insights into the impact of variable parameterisation on z0 estimates for each method.

Two models using the transect method were included in this study, only differing in their definition of the effective obstacle height h (m). Each row of the sub-grid was detrended and treated as a separate transect of length X (m), whereof the transition frequency f from below to above the mean elevation was recorded (often referred to as “zero-up-crossing” in literature) for the calculation of the transect z0 value. Thus, the final z0 value for each sub-grid was then calculated by averaging the individual transect z0 values within the sub-grid. First, the Lettau model calculates h by taking the average vertical extent of the detrended roughness elements, as described by Lettau (1969). Second, the Munro model simplifies Eq. (1) of Lettau (1969) such that the height of roughness element h is calculated by taking twice the standard deviation of elevations along the detrended transect, as described by Munro (1989). In contrast to the study of Munro (1989), which used wind-perpendicular surface transects for the roughness calculation, we used wind-parallel profiles for both models of the transect method. This is because if a crevasse is aligned perpendicular to the prevailing wind direction, a wind-perpendicular transect is not able to detect the crevasse, yielding a relatively low z0 value (for further explanation, see Smith et al.2016). Thus, such an adaptation is essential for heterogeneous, anisotropic and naturally streamlined roughness elements as those investigated in this study since wind systems are influenced by the large-scale catchment topography and therefore often flow up or down the glacier (Quincey et al.2017). The transect method presents a simple approximation of the roughness elements across a profile and, in contrast to the raster method, assumes all roughness elements to be of equal height, uniformly distributed, isotropic and not affected by any sheltering effects.

The raster method is also based on Eq. (1) of Lettau (1969), but the elevation differences between two adjacent cells define the surface roughness in the end. In the raster method all sub-grids were detrended row-wise, and areas below the detrended plane were neglected, assuming that they would be effectively sheltered. Three models following the raster method were included in this study. First, the Smith model was based on the “DEM-based” approach described by Smith et al. (2016). The effective obstacle height h was calculated as the mean elevation of all the cells of each row above the zero plane. The second model of the raster method (hereafter called Chambers model) was based on the “DEM method” described by Chambers et al. (2020). Its only difference to the Smith model is again the definition of the effective obstacle height h, which is twice the standard deviation of elevations above the detrended plane. The third model of the raster method (hereafter called Fitzpatrick model) was based on the “block estimation” described by Fitzpatrick et al. (2019). While the two previous raster methods calculated some of the parameters row-wise, this method followed a moving window approach (see Table 3). The obstacle height h corresponded to the mean of all the detrended elevation values above the zero plane within the window. In this study, a window size of 30 m was chosen, which corresponded to the lower estimation of an average roughness element on the investigated glaciers (see Table 2). For all three models, s and h were calculated individually for each row of the sub-grid (moving window for the Fitzpatrick model). Accordingly, the ground area S was assigned to the area of the sub-grid (moving window), and the final z0 value for each sub-grid was then calculated by taking the mean of its row (moving window) values.

Figure 4Variability in z0 values for Nordenskiöldbreen (2019) depending on the calculation models – Smith (a), Chambers (b), Fitzpatrick (c), Munro (d) and Lettau (e) – for a sub-grid size of 50 m× 50 m and a down-glacier wind direction. Graph (f) shows the DEM with the underlaid hillshade layer.


3 Results

The DEMs obtained from the UAV-based imagery and processed with the SfM-MVS method illustrate that the crevasses of the mapped glaciers are in general aligned perpendicular to the glacier flow direction. The crevasses closer to the front are often deeper and larger in terms of spacing (width of crevasse opening) compared to crevasses located upstream of the glacier (see Fig. 4). The DEMs were then used to calculate the aerodynamic roughness lengths with the two transect method models and the three raster method models. The results (e.g. Figs. 4, 6 and 7) show that the spatial variability in z0 values across the mapped area is up to 3 orders of magnitude. In general, the larger a roughness element is, the larger its aerodynamic roughness length will be. The largest z0 values (decimetre to metre scale) were found close to the glacier front, where crevasses are big and steep, using the transect methods and for winds blowing parallel to the flow direction of the glacier. The lowest values (millimetre to centimetre scale) are estimated with the raster methods for smooth, crevasse-free ice and for cross-glacier wind directions.

3.1 Model results

Figure 4 shows the estimated aerodynamic roughness length values for each of the five applied models for the down-glacier wind direction on Nordenskiöldbreen (2019). The results reveal that all models agree on the relative spatial z0 patterns across the glacier. Accordingly, for the flatter and less crevassed part (e.g. to the right of the mapped area on Nordenskiöldbreen), all models show lower sub-grid z0 values (red) compared to the heavily crevassed part close to the glacier front (yellow). To further investigate the relative agreement between the models a statistical correlation test on the z0 data of Nordenskiöldbreen (averaged over all four wind directions) was conducted. The correlation test showed that all models are strongly correlated to each other, leading to R2 values of 0.877 and higher (Dachauer2020). A similar correlation was also observed for the other glaciers. While this indicates that the z0 values are correlated, it does not provide any conclusion about the quality of the individual models.

Table 4Overview of mean z0 values (m) for each glacier and model either for the down-glacier or cross-glacier (left-to-right) wind direction.

Download Print Version | Download XLSX

Figure 5Variability in z0 values for Fridtjovbreen depending on the wind direction for a sub-grid size of 50 m calculated with the Smith model. Winds blowing across the glacier either from the left-to-right (a) or right-to-left (c) produce smaller z0 values than down-glacier (b) or up-glacier (d) wind systems.


Figure 4 further illustrates that the absolute values of the Munro and Lettau models show greater roughness values on the same sub-grid than the other three models. In more detail, the Lettau model estimates are generally greater than those of the Munro model and Chamber z0 values greater than those of the Smith model (see also Table 4). The three models of the raster method provide sub-grid z0 values of about 0.001 m for slightly crevassed areas and up to 1 m for heavily crevassed areas. The same sub-grids calculated with the two transect methods produced z0 values which are locally up to 1 order of magnitude larger. Table 4 illustrates the down-glacier and cross-glacier (left-to-right) mean z0 values for all glaciers and models. Table 4 also shows that average z0 values across the glacier vary almost up to half an order of magnitude between the models. In summary, despite the clear relative agreement between the models, the estimated magnitude of the z0 values varies substantially between the models – especially between the raster and transect method models. This might be explained by the fact that the transect method does not account for sheltering of an obstacle (Smith et al.2016). The raster method on the other hand assumes the areas below the detrended plane to be effectively sheltered. This plane indicates how far the effective turbulent mixing advances into the crevasses, and the corresponding z0 values are expected to be lower if the sheltering effect is considered (Nicholson et al.2016).

3.2 Wind direction variability

Figure 5 illustrates the map of Fridtjovbreen with z0 values obtained with the Smith method for all four wind directions. The results show that z0 values are higher for wind directions that face the crevasses perpendicularly (i.e. up- and down-glacier) and lower for wind that blows parallel to the elongated crevasse features (i.e. cross-glacier). The two cross-glacier (up- and down-glacier) wind directions lead to very similar z0 values since they are both calculated on the same transect but from opposite wind directions. Since crevasses are mainly oriented perpendicular to the glacier flow direction, the mapped areas show a strongly anisotropic pattern of the glacier surface. This wind dependency effect is visible on all five applied models and is independent of the roughness element size. However, Fig. 5 illustrates that larger roughness elements, which can be found close to the glacier front for instance, present a stronger wind dependency because they vary more strongly with changing wind directions (from decimetre to metre scale) compared to areas that are less crevassed like on the upper part of Fridtjovbreen (similar z0 values in millimetre scale for all wind directions).

Figure 6Boxplot visualisation of sub-grid z0 values for all four wind directions and each applied model determined on Fridtjovbreen. The wind direction is either down-glacier (orange), up-glacier (green), cross-glacier from right to left (blue; rtl) or cross-glacier from left to right (red; ltr). Whiskers are visualising the variability outside the upper and lower quartiles up to 1.5 times the interquartile range.


Figure 7Variability in z0 values for the glaciers Tunabreen (a) and Heuglinbreen (b) calculated with the Smith model for a down-glacier wind direction and a sub-grid size of 50 m.


Additionally, Fig. 6 shows the boxplot graph of z0 estimates illustrating the wind direction dependency of z0 values on Fridtjovbreen. The boxplot z0 medians of all models vary in a range of 0.012 to 0.037 m for crosswind directions and 0.058 to 0.2 m for up- and down-glacier wind systems. The mean and median values between up- and down-glacier wind systems (left-to-right and right-to-left crosswinds, respectively) never differ by more than 10 % and only rarely more than 5 % independently of the chosen model or glacier. In summary, the wind direction has a large impact on the resulting aerodynamic roughness length values. Its effect on average or median z0 values slightly exceeds the model variability. Locally, however, both mean and median z0 estimates can vary by about 1 order of magnitude with changing wind direction.

3.3 Glacier variability

All glaciers showed a similar range of z0 values on the decimetre to metre scale for heavily crevassed areas and millimetre to centimetre scale for less crevassed areas (see Figs. 4, 5 and 7). However, the results of Table 4 show that the mean z0 values for the Smith model and down-glacier wind direction are somewhat larger on Tunabreen (0.223 m) and the extract of Nordenskiöldbreen mapped in 2020 (0.231 m) compared to the other glaciers. The same mean z0 values on Nordenskiöldbreen measured in 2019 are lower, at only 0.129 m. Heuglinbreen (0.15 m) and especially Fridtjovbreen (0.082 m) show lower mean z0 values than the glaciers mentioned above. Nevertheless, all the patterns found in this study (e.g. between wind direction or sub-grid size and the z0 values) are independent of the glacier and vary (if even) only in magnitude rather than relative patterns in between the glaciers.

On a side note, it is interesting to study the results of Nordenskiöldbreen 2019 and 2020 in more detail. The comparison provides insights into the inter-annual temporal variability in z0 values for 2 consecutive years. The mapped area in summer 2020 (1 d of fieldwork) was smaller compared to the field campaign approximately 1 year earlier in 2019 (3 d of fieldwork). Therefore, only the overlapping area has been used for the temporal variability investigation. Furthermore, for this particular comparison the DEM extract of Nordenskiöldbreen 2019 (0.17 m px−1) was resampled to the resolution of the Nordenskiöldbreen 2020 DEM (0.22 m px−1). The z0 values on Nordenskiöldbreen were very similar for the 2 consecutive years. The mean z0 value for the Smith model and down-glacier wind direction in 2019 was, at 0.25 m, only slightly greater than the corresponding value of 0.23 m at the same area 1 year later. The observations are in line with other studies (i.e. Fitzpatrick et al.2019), which also did not observe a large difference in z0 estimates for the same location over 2 consecutive years. Relatively, the mean z0 values never deviated more than 20 % and mostly less than 10 % between the 2 years. In summary, the differences in z0 estimates across the models exceed the inter-annual temporal z0 variability by far, independently of the model calculation.

Figure 8Scale-dependent z0 values for the Smith model and down-glacier wind direction applied on Fridtjovbreen for sub-grid size resolutions of 5 m (a), 10 m (b), 30 m (c), 40 m (d), 50 m (e), 60 m (f), 70 m (g) and 100 m (h) with underlaid hillshade layer for orientation.


3.4 Sub-grid size dependency

The mapped glacier areas were divided into sub-grids with a grid size of 50 m× 50 m to account for the spatial variability in z0 on the glacier. However, to investigate the grid size dependency of the sub-grids on z0 values, a small case study on the glacier Fridtjovbreen was conducted. The Smith model with a down-glacier wind direction was used to calculate aerodynamic roughness lengths for sub-grid sizes between 5 and 150 m. Figure 8 illustrates the scale dependency of z0 values to the selected sub-grid size. The results show that a grid size of 5 to 10 m mostly produces z0 values on a scale of centimetres. Between sub-grid sizes of 10 and 50 m a higher grid size results in higher z0 values. From 50 m onwards, the z0 values are mostly at the decimetre scale and do not change substantially. This behaviour was also observed with the other models and for all wind directions. For small sub-grid sizes, the z0 values are likely representing microtopography rather than the macro-scale surface roughness of the crevasses.

Figure 9Scale dependency of mean z0 values for the five applied models and a down-glacier wind direction on Fridtjovbreen with chosen sub-grid sizes from 5 to 150 m. The black line indicates the chosen grid size of 50 m which incorporates the upper boundary of the average obstacle sizes between 30 and 50 m found on all glaciers (see Table 2).


Figure 9 summarises the mean z0 results of the mentioned case study for the down-glacier wind direction. For all the models, the mean z0 values increased by at least half an order of magnitude between a sub-grid size of 5 and 150 m. All models showed a similar pattern with strongly increasing z0 values for small grid sizes (5 to 30 m) and only slightly increasing estimates afterwards. Grid sizes of 70 m or more lost their strong scale dependency effect, leading to stable z0 estimates. The same effect is visible in Fig. 8, where the chosen grid size of 50 m represents similar z0 values (colours) as the higher grid sizes for the same location while still providing a considerable grid resolution.

4 Discussion

4.1 Model inputs for aerodynamic roughness length estimation

4.1.1 Validation of digital elevation model accuracy

The obtained DEM resolution with the SfM-MVS method (about 0.25 m) was accurate enough to capture the large crevasse structures. Given the advantages of a UAV compared to other devices (e.g. low cost, applicable to inaccessible areas; Hann et al.2021), this study shows that UAVs provide a reliable and effective way of data gathering for aerodynamic roughness length estimation on glaciers. Nevertheless, the depth of the crevasses must be seen as a minimum depth, and the crevasses might penetrate further into the glacier than actually measured. This is due to snow bridges or the lack of reflected light from the deep crevasses. The latter prevents the SfM-MVS methods from correctly constructing the deeper parts of the crevasse. Additionally, the lowest points of the crevasses are very narrow and may not be captured accurately (Ryan et al.2015). However, in an aerodynamic context those narrow crevasses are not likely to have a significant influence on the heat exchange since they lie below the penetration depth of effective turbulent mixing (Nicholson et al.2016). Equation (1) of Lettau (1969) and the transect methods do not define any penetration depth limit. The raster method, however, assumes that effective roughness only depends on the roughness elements above the detrended plane level, which indicates how far the effective turbulent mixing advances into the crevasses.

In the scope of this study the absence of GCPs was also considered. GCPs can significantly increase the georeferencing accuracy of the DEMs (Chudley et al.2019). However, it was practically impossible to place GCPs on the crevassed glacier surfaces due to safety reasons. Therefore, the georeferencing information was provided only from the on-board GPS and the measurements of camera orientation (James et al.2017). In other words, the positional accuracy of the DEM is limited by the internal GPS system of the UAV, which has a relatively low accuracy (Federman et al.2017). However, the objective of this study was not to obtain a high-accuracy DEM but rather a precise DEM combined with a detrending approach for the investigation of the effect of relative distances. Therefore, the given hover accuracy of ± 1.5 m horizontally and ± 0.5 m vertically for both UAVs (DJI2017, 2019) was considered sufficient (Dachauer2020).

Nevertheless, a comparison with Sentinel-2 satellite data (ESA2020) revealed horizontal positioning errors (data not shown). The deviation was classified as a systematic offset (for more details see Dachauer2020). The detected small horizontal distortions of about 1.7 % (horizontal length deviation in percent of DEM compared to the Sentinel-2 satellite image) mainly occurred on the ends of the DEMs. This is a typical feature appearing when only using nadir imagery and is related to self-calibration because the reconstruction software is not able to derive the accurate radial lens distortion leading to a systematic “doming” DEM deformation (James and Robson2014). However, a systematic error is of low significance for this study since only the relative accuracy of the roughness elements (i.e. distortion) influences the estimation of the aerodynamic roughness lengths. Furthermore, the influence of a small distortion of a few percent or several metres across the mapped area on the resulting z0 values is minor compared to other parameters such as wind direction, model calculation or scale dependency. For example, a distortion of 2 % led to a change in z0 of about 4 % (Dachauer2020). This means that the obtained DEMs in this study are a reliable data source to estimate the aerodynamic roughness length.

4.1.2 Scale dependency

Many studies have investigated and encountered the dependency of z0 estimates on the size of the sub-grid or the transect length (e.g. Fitzpatrick et al.2019; Miles et al.2017; Rees and Arnold2006) and reported that larger sub-grids or longer transects cause z0 values to increase (Chambers et al.2020). Also the case study conducted on Fridtjovbreen data revealed that the z0 values increase with larger grid size independently of the chosen glacier, model or wind direction. This can be explained by the fact that glaciers often have heterogeneous roughness elements (Quincey et al.2017). In general, it needs to be considered that the selection of an appropriate grid size comes with a large potential uncertainty. To find the most meaningful grid size a comparison with independent methods like aerodynamic wind profiles is recommended (Chambers et al.2020). Since this option was not available in our study, the validity of the chosen grid size has been evaluated theoretically. The grid size in our models corresponds to the ground area S, whose definition requires the grid size to be the size of one individual roughness element (Lettau1969), which is 30 to 50 m (see Table 2). According to Smith (2014), the definition of “roughness” is related to the grid scale that separates the grain roughness (representing the texture of a roughness element) and the form roughness (corresponding to the form drag of the roughness element itself). The grid size at which the transition from grain to form roughness occurs provides a useful reference point and can be found as a kink in the trend line of a figure plotting z0 values against grid sizes (Shepard et al.2001; Smith2014). Accordingly, given the large roughness elements investigated in our study, the “grain” roughness was assumed to belong to the texture on the crevasses. Thus, the transition from grain to form roughness is again located somewhere between a sub-grid size of 30 to 50 m (see Fig. 9). Therefore, in this study the grid size of 50 m was chosen since it is the smallest resolution possible to still include the size of an average obstacle. Typical z0 estimates for smooth glacier ice, which for instance can be found on the upper part of Fridtjovbreen, have a length of about 1 mm (Brock et al.2006). The choice of a 50 m grid size can further be justified since Fig. 8 shows that grid sizes below 30 m do not provide high enough values and therefore do not agree with literature values.

Figure 10Anisotropy ratio values for the two models of Smith (blue) and Munro (green) calculated for sub-grid z0 values on Fridtjovbreen. The dark green colour corresponds to an overlapping of both models. A positive ratio towards 1 means that parallel winds (up- and down-glacier) are dominant over perpendicular winds (cross-glacier).


4.2 Estimated aerodynamic roughness length

4.2.1 Wind direction dependency

The wind direction has a large impact on the magnitude of the z0 values on glaciers since they contain many anisotropic roughness elements (Chambers et al.2020). Our results confirmed this statement and further revealed that larger roughness elements in general present a stronger wind dependency (see Fig. 5). The wind dependency effect was additionally investigated with the calculation of the anisotropy ratio Ω (see Smith et al.2006, for further explanation):

(2) Ω = z 0 - z 0 z 0 + z 0 ,

where subscripts and  denote parallel and perpendicular to the ice flow direction, respectively. Figure 10 illustrates a histogram of the sub-grid's Ω values on Fridtjovbreen for the Smith (blue) and Munro (green) model. The results show that the sub-grid z0 values are strongly anisotropic all across the glacier, whereas wind directions parallel to the flow direction are dominant (since values are mostly positive). Both the raster (here Smith model) and the transect methods (Munro) show a similar pattern with the highest frequency at the range of 0.4 to 0.7 (Smith) and 0.4 to 0.8 (Munro). Although Tunabreen has still mostly positive ratio values, it is the glacier with the least anisotropic behaviour of all investigated glaciers. The importance of wind direction can be observed in many studies (e.g. Fitzpatrick et al.2019; Munro1989; Smith et al.2016) and is found to be the strongest on ablation zones where elongated features like meltwater channels and crevasses are frequent. Thus, wind directions that face these features perpendicularly lead to higher z0 values due to an increased form drag (Fitzpatrick et al.2019). Winds are in general very likely to blow in the direction of the mean slope angle or in a down-slope wind direction (e.g. Fitzpatrick et al.2019; Karner et al.2013). This is due to katabatically forced down-slope winds, which are common over the glacier terminus (Munro1989). Nevertheless, Esau and Repina (2012) found the katabatically forced wind systems to be of less significance for polar tidewater glaciers. This highlights the influence and importance of wind direction on effective z0 values. Furthermore, the wind direction dependency reveals the temporal variability in the aerodynamic roughness length due to changing wind directions from daily up to seasonal timescales.

4.2.2 Variability across the glaciers

Since almost no values of z0 for heavily crevassed glaciers are available in the literature, the validation of the results of this study is difficult. The roughness elements investigated by previous studies (see Table 1) were smaller than the crevasses of this study. Therefore, it is expected that the z0 values obtained here should be larger. This was the case, although there was a significant spread of z0 values across the estimation models. The mean results in Table 4, and there in particular the raster method results, fall mostly within the same order of magnitude. The range suggested by Fitzpatrick et al. (2019) fits most of the results from this study across different methods and glaciers. Outside the field of glaciology, the heavily crevassed glaciers might most effectively be compared with villages since buildings can have similar obstacle density and height. Accordingly, the z0 values for villages are about 0.2–0.4 m (Abbas et al.2021), which lie within the range of estimated roughness values in this study.

In general, Tunabreen and Nordenskiöldbreen (especially the part mapped in 2020) show higher z0 values than Fridtjovbreen and Heuglinbreen. This is in good agreement with the average length and height of roughness elements estimated for each glacier in Table 2. The crevasses on Tunabreen and Nordenskiöldbreen are generally deeper and steeper than the crevasses on the other two glaciers. This might be explained by different dynamical behaviour, leading to higher z0 values. In general, a faster flowing glacier leads to more crevasses in the terminus area of the glacier and a larger area of crevassing (Błaszczyk et al.2009). Additionally, Tunabreen was recently surging, which usually leads to very chaotically aligned crevasses (Mansell et al.2012). Thus, these results are in line with the observation that the anisotropy ratio value Ω is lowest on Tunabreen as perfectly aligned roughness elements increase the anisotropy effect. However, comparing average z0 values among different glaciers is challenging because the mean z0 values depend substantially on the mapped input area and therefore the size of included roughness elements.

4.3 Model performance assessment

If the roughness elements on a plot are too densely packed, then the objects begin to interfere aerodynamically with each other (Rounce et al.2015). They form a plateau-like new surface at their tops (Lettau1969) leading to a skimming flow. If this roughness density (frontal area divided by ground area) increases as far as is necessary to induce skimming flow, then the z0 values decrease (Grimmond and Oke1999). Accordingly, the results from the transect methods which are not considering any sheltering effects are likely to overestimate the roughness, especially for densely packed obstacles (Nicholson et al.2016). Several studies (e.g. Nicholson et al.2016; Rounce et al.2015; Smith et al.2016) used a roughness density threshold of 20 % to 30 % (as stated by Macdonald et al.1998) for Eq. (1) of Lettau (1969) to still be valid. In our study, all glaciers were tested with the Smith model for their roughness density. Apart from some single heavily crevassed sub-grids close to the glacier front which exceed the threshold of 30 %, all glaciers are below the given threshold with mean roughness density values between 0.1 and 0.15 for each glacier, indicating no skimming flow over the obstacles was likely. Therefore, this study assumes that Eq. (1) of Lettau (1969) also holds for heavily crevassed tidewater glaciers with respect to the sheltering effect.

The Lettau formula is based on empirical experiments of systematically placed bushel baskets for roughness simulation. Transferring the simple relation onto heterogeneous and complex surfaces as found on the heavily crevassed tidewater glaciers might lead to some uncertainties. This is because the Lettau relation for the calculation of z0 strongly simplifies the complex surface roughness and its obstacle size, shape and density. A simplified representation of the surface roughness might fail to capture the complete range of aerodynamic processes. The Smith and Chambers models, as well as Munro and Lettau models, only differ in the definition of h. An appropriate definition of h is crucial since the obstacle height is the most important control parameter over the output of z0 (Nield et al.2013b). Thus, the lack of a clear obstacle definition presents the main problem of the bulk method (Rounce et al.2015), especially in crevassed glacier areas, where an apparent base level is missing (Nicholson et al.2016). To find out which model and which definition of h might perform best, a comparison with alternative measurement methods (e.g. wind profiles) would be necessary, as is done in several studies (Fitzpatrick et al.2019; Quincey et al.2017; Smeets et al.1999). For this study, however, the area of interest was inaccessible and so this validation option was impossible. In any case, it clearly can be seen that the definition of h has a lower impact on the results than the choice of the basic method (i.e. raster/transect), which particularly affects the values of the silhouette area s. The definition of the ground area S which corresponds to the chosen grid size and the according profile length is a simple approximation of the real fetch footprint. However, the simplification allows for the estimation of z0 values for the four wind directions and additionally provides uniform parameterisation throughout all glaciers and models. The widely adopted drag coefficient of Lettau (1969) of cd=0.5 corresponds to an average form drag effect on roughness elements. Its rationale is widely discussed since it does not necessarily hold for complex and heterogeneous locations where drag may be large due to the irregular density and shape of roughness elements. This might lead to a higher effective drag and therefore underestimated z0 values (Quincey et al.2017).

Turbulent fluxes already contribute a large fraction of surface ice melt (Fausto et al.2016) and are supposed to have an even greater contribution to surface energy balance models under global warming (Smith et al.2020). Thus, a high accuracy of estimated z0 values is desirable. The validation of the model estimates remains a big challenge due to the lack of reference values, and it should be questioned whether the modelled z0 values are accurate enough for energy balance models. In general, an increase in z0 by 1 order of magnitude will more than double the value of turbulent fluxes (Brock et al.2000). Therefore, whilst grid-size and model choice can contribute up to 1 order of magnitude uncertainty in z0 values, it is far outweighed by the 3 orders of magnitude range in z0 observed in the present study. Spatial variability in z0 values caused by intense crevassing near the margin of tidewater glaciers therefore greatly exceeds the uncertainty introduced by modelling choices and can be constrained with sufficient accuracy through the methods pioneered here. Furthermore, this large spatial variability, which is commonly not considered in models, affects the partition of mass loss between frontal and surface ablation. Additionally, the relevance of further narrowing down the modelled range in the future depends a lot on the field of application. For large-scale, satellite-based investigations an average value between all models (e.g. 0.1 m) for heavily crevassed glacier areas might be a sufficiently accurate approximation. Small-scale investigations on individual glaciers on the other hand likely benefit from more accurate z0 estimates. It is here in particular, where our study shows that UAVs are the ideal platform for investigating aerodynamic roughness length.

5 Conclusions

The heavily crevassed terminus areas of the tidewater glaciers Fridtjovbreen, Heuglinbreen, Nordenskiöldbreen and Tunabreen were mapped with UAVs to build DEMs that revealed crevasse shape information. To take into account the spatial variability across the glacier, the DEM was divided into sub-grids of 50 m× 50 m, which was assumed to be large enough to include an average roughness element while still being small enough to account for the roughness variability across the glacier. Five different models (belonging to either the transect or the raster method) were applied to each DEM sub-grid to calculate the aerodynamic roughness length. The z0 estimates from the transect method were in general greater (up to 1 order of magnitude) than the raster method estimations. Wind direction and sub-grid size had a large impact on the z0 estimates, again producing a range of up to 1 order of magnitude for each parameter. Winds blowing parallel to the ice flow direction produced larger z0 values than cross-glacier winds. The chosen sub-grid size presents a large uncertainty in aerodynamic roughness length estimation. The resulting z0 values are strongly scale dependent such that a larger grid size leads to greater z0 values. If all parameters (i.e. model, wind direction, grid size) are included, the spread of the resulting z0 estimates is large, ranging from below 1 mm (snow-covered, smooth glacier surface) up to several decimetres (heavily crevassed ice) or locally even more. Averaged z0 values for down-glacier wind directions varied from 0.08 to 0.88 m depending on glacier and model. Nevertheless, all models managed to detect the same spatial variability across the glacier. The UAV approach allows several z0 values for each mapped glacier area to be derived, which is crucial for heavily crevassed glaciers since one value would be a poor representation of the real roughness across such a complex topography. Therefore, models can now incorporate distributed z0 estimates easily following UAV deployment, potentially leading to a better representation of turbulent heat fluxes and prediction of surface ice melt rates.

Spatial and temporal variability in crevassing and a dependence on wind direction were found to extend the range of z0 values on tidewater glaciers. Variability caused by sub-grid size and model calculation assumptions reveal uncertainties which should be addressed by future investigations. Some degree of uncertainty also comes with the unsatisfactory georeferencing of the DEM in crevassed areas because the inaccessible topography imposes practical limitations, especially on the use of GCPs. Furthermore, future work should seek a scale-independent method for z0 calculation and also assess model performance using meteorological measurements (e.g. wind profiles) or computational fluid dynamics simulations. In a next step, a combination of high-resolution DEMs from UAVs for reference z0 values and satellite-based crevasse density estimates might prove valuable for future research.

Code availability

The codes used in this paper can be downloaded from (Dachauer2021).

Data availability

The data used in this paper can be downloaded from (Dachauer et al.2021). Additional data can be obtained from the authors without conditions.

Author contributions

AD was responsible for the formal analysis, the investigation and visualisation of the paper, and the original draft preparation. RH accomplished the project administration and data curation. Furthermore, he supervised the work and managed the resources together with AJH. All authors contributed to the conceptualisation and the review writing and have read and agreed to the published version of the manuscript.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We much appreciate the help of Evan Stewart Miles by sharing his codes and knowledge.

Financial support

Fieldwork was funded by the Arctic Field Grant (RiS ID: 11148), provided by the Svalbard Science Forum (SSF) which is coordinated by the Research Council of Norway (RCN). The work is partly sponsored by the RCN through the Centre of Excellence funding scheme, project number 223254, AMOS, and project CLIMAGAS with project number 294764.

Review statement

This paper was edited by Brice Noël and reviewed by Evan Miles and Maurice Van Tiggelen.


Abbas, M. R., Bin Rasib, A. W., Abbas, T. R., Ahmad, B. B., and Dutsenwai, H. S.: Assessment of Aerodynamic Roughness Length Using Remotely Sensed Land Cover Features and MODIS, IOP C. Ser. Earth Env., 722, 012015,, 2021. a

Agisoft LLC: Agisoft Metashape, available at:, last access: 27 July 2020. a

Bhardwaj, A., Sam, L., Akanksha, Martín-Torres, F. J., and Kumar, R.: UAVs as remote sensing platform in glaciology: Present applications and future prospects, Remote Sen. Environ., 175, 196–204,, 2016. a, b

Błaszczyk, M., Jania, J. A., and Hagen, J. O.: Tidewater glaciers of Svalbard: Recent changes and estimates of calving fluxes, Pol. Polar Res., 30, 85–142, 2009. a

Brock, B. W., Willis, I. C., Sharp, M. J., and Arnold, N. S.: Modelling seasonal and spatial variations in the surface energy balance of Haut Glacier d'Arolla, Switzerland, Ann. Glaciol., 31, 53–62,, 2000. a

Brock, B. W., Willis, I. C., and Sharp, M. J.: Measurement and parameterization of aerodynamic roughness length variations at Haut Glacier d'Arolla, Switzerland, J. Glaciol., 52, 281–297,, 2006. a, b, c

Casella, V. and Franzini, M.: Modelling steep surfaces by various configurations of nadir and oblique photogrammetry, ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci., III-1, Prague, Czech Republic, 175–182,, 2016. a

Chambers, J. R., Smith, M. W., Quincey, D. J., Carrivick, J. L., Ross, A. N., and James, M. R.: Glacial Aerodynamic Roughness Estimates: Uncertainty, Sensitivity, and Precision in Field Measurements, J. Geophys. Res.-Earth, 125, 1–19,, 2020. a, b, c, d, e, f, g

Chappell, A. and Heritage, G. L.: Using illumination and shadow to model aerodynamic resistance and flow separation: An isotropic study, Atmos. Environ., 41, 5817–5830,, 2007. a

Chudley, T. R., Christoffersen, P., Doyle, S. H., Abellan, A., and Snooke, N.: High-accuracy UAV photogrammetry of ice sheet dynamics with no ground control, The Cryosphere, 13, 955–968,, 2019. a, b

Dachauer, A.: Aerodynamic Roughness Length of Crevassed Tidewater Glaciers from UAV Mapping, MS thesis, p. 100, available at: (last access: 7 December 2021), 2020. a, b, c, d, e

Dachauer, A.: ArminDach/z0_UAVs: (Version v1), Zenodo [code],, 2021. a

Dachauer, A., Hann, R., and Hodson, A.: z0 of crevassed tidewater glaciers, DataverseNO [data set],, 2021. a

DJI: Phantom 4 Pro/Pro+, User Manual v1.4, Tech. rep., DJI, Shenzhen, China, 2017. a

DJI: Mavic 2 Enterprise Series, User Manual v1.8, Tech. rep., DJI, Shenzhen, China, 2019. a

Ericson, Y., Falck, E., Chierici, M., Fransson, A., and Kristiansen, S.: Marine CO2 system variability in a high arctic tidewater-glacier fjord system, Tempelfjorden, Svalbard, Cont. Shelf Res., 181, 1–13,, 2019. a

ESA: Copernicus Open Access Hub, available at:, last access: 23 August 2020. a, b

Esau, I. and Repina, I.: Wind climate in kongsfjorden, svalbard, and attribution of leading wind driving mechanisms through turbulence-resolving simulations, Adv. Meteorol., 2012, 568454,, 2012. a

Fausto, R. S., Van As, D., Box, J. E., Colgan, W., Langen, P. L., and Mottram, R. H.: The implication of nonradiative energy fluxes dominating Greenland ice sheet exceptional ablation area surface melt in 2012, Geophys. Res. Lett., 43, 2649–2658,, 2016. a

Federman, A., Santana Quintero, M., Kretz, S., Gregg, J., Lengies, M., Ouimet, C., and Laliberte, J.: UAV photgrammetric workflows: A best practice guideline, Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., Ottawa, Canada, XLII-2/W5, 237–244,, 2017. a

Fitzpatrick, N., Radić, V., and Menounos, B.: A multi-season investigation of glacier surface roughness lengths through in situ and remote observation, The Cryosphere, 13, 1051–1071,, 2019. a, b, c, d, e, f, g, h, i, j, k

Grimmond, C. S. B. and Oke, T. R.: Aerodynamic Properties of Urban Areas Derived from Analysis of Surface Form, J. Appl. Meteorol., 38, 1262–1292,<1262:APOUAD>2.0.CO;2, 1999. a

Hagen, J. O., Liestøl, O., Roland, E., and Jørgensen, T.: Glacier Atlas of Svalbard and Jan Mayen, Norsk polarinstitutt, Oslo, 1993. a

Hann, R., Altstädter, B., Betlem, P., Deja, K., Dragańska-Deja, K., Ewertowski, M., Hartvich, F., Jonassen, M., Lampert, A., Laska, M., Sobota, I., Storvold, R., Tomczyk, A., Wojtysiak, K., and Zagórski, P.: Scientific Applications of Unmanned Vehicles in Svalbard, SESS report 2020, Svalbard Integrated Arctic Earth Observing System, Zenodo,, 2021. a, b

Hanssen-Bauer, I., Førland, E. J., Hisdal, H., Mayer, S., Sandø, A. B., Sorteberg, A., Adakudlu, M., Andresen, J., Bakke, J., Beldring, S., Benestad, R., Bilt, W., Bogen, J., Borstad, C., Breili, K., Breivik, Ø., Børsheim, K. Y., Christiansen, H. H., Dobler, A., Engeset, R., Frauenfelder, R., Gerland, S., Gjelten, H. M., Gundersen, J., Isaksen, K., Jaedicke, C., Kierulf, H., Kohler, J., Li, H., Lutz, J., Melvold, K., Mezghani, A., Nilsen, F., Nilsen, I. B., Nilsen, J. E. Ø., Pavlova, O., Ravndal, O., Risebrobakken, B., Saloranta, T., Sandven, S., Schuler, T. V., Simpson, M. J. R., Skogen, M., Smedsrud, L. H., Sund, M., Vikhamar-Schuler, D., Westermann, S., and Wong, W. K.: Climate in Svalbard 2100, The Norwegian Centre for Climate Services (NCCS), Trondheim, p. 105, 2019. a, b

Hock, R.: Glacier melt: A review of processes and their modelling, Prog. Phys. Geog., 29, 362–391,, 2005. a

Irvine-Fynn, T. D., Sanz-Ablanedo, E., Rutter, N., Smith, M. W., and Chandler, J. H.: Measuring glacier surface roughness using plot-scale, close-range digital photogrammetry, J. Glaciol., 60, 957–969,, 2014. a, b, c

James, M. R. and Robson, S.: Mitigating systematic error in topographic models derived from UAV and ground-based image networks, Earth Surf. Proc. Land., 39, 1413–1420,, 2014. a

James, M. R., Robson, S., and Smith, M. W.: 3-D uncertainty-based topographic change detection with structure-from-motion photogrammetry: precision maps for ground control and directly georeferenced surveys, Earth Surf. Proc. Land., 42, 1769–1788,, 2017. a

Karner, F., Obleitner, F., Krismer, T., Kohler, J., and Greuell, W.: A decade of energy and mass balance investigations on the glacier Kongsvegen, Svalbard, J. Geophys. Res.-Atmos., 118, 3986–4000,, 2013. a

Lettau, H.: Note on aerodynamic roughness-parameter estimation on the basis of roughness-element description, J. Appl. Meteorol., 8, 828–832, 1969. a, b, c, d, e, f, g, h, i, j, k, l

Macdonald, R. W., Griffiths, R. F., and Hall, D. J.: An improved method for the estimation of surface roughness of obstacle arrays, Atmos. Environ., 32, 1857–1864,, 1998. a

Mansell, D., Luckman, A., and Murray, T.: Dynamics of tidewater surge-type glaciers in northwest Svalbard, J. Glaciol., 58, 110–118,, 2012. a

Miles, E. S., Steiner, J. F., and Brun, F.: Highly variable aerodynamic roughness length (z0) for a hummocky debris-covered glacier, J. Geophys. Res.-Atmos., 122, 8447–8466,, 2017. a, b

Munro, D. S.: Surface roughness and bulk heat transfer on a glacier: comparison with eddy correlation, J. Glaciol., 35, 343–348,, 1989. a, b, c, d, e, f

Murray, T., James, T., MacHeret, Y., Lavrentiev, I., Glazovsky, A., and Sykes, H.: Geometric changes in a tidewater glacier in svalbard during its surge cycle, Arct. Antarct. Alp. Res., 44, 359–367,, 2012. a

Nicholson, L. I., Pȩtlicki, M., Partan, B., and MacDonell, S.: 3-D surface properties of glacier penitentes over an ablation season, measured using a Microsoft Xbox Kinect, The Cryosphere, 10, 1897–1913,, 2016. a, b, c, d, e, f

Nield, J. M., Chiverrell, R. C., Darby, S. E., Leyland, J., Vircavs, L. H., and Jacobs, B.: Complex spatial feedbacks of tephra redistribution, ice melt and surface roughness modulate ablation on tephra covered glaciers, Earth Surf. Proc. Land., 38, 95–102,, 2013a. a, b

Nield, J. M., King, J., Wiggs, G. F. S., Leyland, J., Bryant, R. G., Chiverrell, R. C., Darby, S. E., Eckardt, F. D., Thomas, D. S. G., Vircavs, L. H., and Washington, R.: Estimating aerodynamic roughness over complex surface terrain, J. Geophys. Res.-Atmos., 118, 12948–12961,, 2013b. a

Obleitner, F.: The energy budget of snow and ice at Breidamerkurjokull, Vatnajokull, Iceland, Bound.-Lay. Meteorol., 97, 385–410,, 2000. a

Quincey, D., Smith, M., Rounce, D., Ross, A., King, O., and Watson, C.: Evaluating morphological estimates of the aerodynamic roughness of debris covered glacier ice, Earth Surf. Proc. Land., 42, 2541–2553,, 2017. a, b, c, d, e, f, g

Rees, W. G. and Arnold, N. S.: Scale-dependent roughness of a glacier surface: Implications for radar backscatter and aerodynamic roughness modelling, J. Glaciol., 52, 214–222,, 2006. a, b

Rounce, D. R., Quincey, D. J., and McKinney, D. C.: Debris-covered glacier energy balance model for Imja–Lhotse Shar Glacier in the Everest region of Nepal, The Cryosphere, 9, 2295–2310,, 2015. a, b, c

Ryan, J. C., Hubbard, A. L., Box, J. E., Todd, J., Christoffersen, P., Carr, J. R., Holt, T. O., and Snooke, N.: UAV photogrammetry and structure from motion to assess calving dynamics at Store Glacier, a large outlet draining the Greenland ice sheet, The Cryosphere, 9, 1–11,, 2015. a

Shepard, M. K., Campbell, B. A., Bulmer, M. H., Farr, T. G., Gaddis, L. R., and Plaut, J. J.: The roughness of natural terrain: A planetary and remote sensing perspective, J. Geophys. Res.-Planet., 106, 32777–32795,, 2001. a

Smeets, C. J., Duynkerke, P. G., and Vugts, H. F.: Observed wind profiles and turbulence fluxes over an ice surface with changing surface roughness, Bound.-Lay. Meteorol., 92, 101–123,, 1999. a, b, c

Smith, B. E., Raymond, C. F., and Scambos, T.: Anisotropic texture of ice sheet surfaces, J. Geophys. Res.-Earth, 111, 1–8,, 2006. a

Smith, M. W.: Roughness in the Earth Sciences, Earth-Sci. Rev., 136, 202–225,, 2014. a, b, c, d

Smith, M. W., Quincey, D. J., Dixon, T., Bingham, R. G., Carrivick, J. L., Irvine-Fynn, T. D. L., and Rippin, D. M.: Aerodynamic roughness of glacial ice surfaces derived from high-resolution topographic data, J. Geophys. Res.-Earth, 121, 748–766,, 2016. a, b, c, d, e, f, g, h, i, j, k

Smith, T., Smith, M. W., Chambers, J. R., Sailer, R., Nicholson, L., Mertes, J., Quincey, D. J., Carrivick, J. L., and Stiperski, I.: A scale-dependent model to represent changing aerodynamic roughness of ablating glacier ice based on repeat topographic surveys, J. Glaciol., pp. 1–15,, 2020. a, b

Uysal, M., Toprak, A. S., and Polat, N.: DEM generation with UAV Photogrammetry and accuracy analysis in Sahitler hill, Measurement, 73, 539–543,, 2015.  a, b

van Tiggelen, M., Smeets, P. C. J. P., Reijmer, C. H., Wouters, B., Steiner, J. F., Nieuwstraten, E. J., Immerzeel, W. W., and van den Broeke, M. R.: Mapping the aerodynamic roughness of the Greenland Ice Sheet surface using ICESat-2: evaluation over the K-transect, The Cryosphere, 15, 2601–2621,, 2021. a, b

Verhoeven, G.: Taking computer vision aloft – archaeological three-dimensional reconstructions from aerial photographs with photoscan, Archaeol. Prospect., 18, 67–73,, 2011. a

Short summary
This study investigated the aerodynamic roughness length (z0) – an important parameter to determine the surface roughness – of crevassed tidewater glaciers on Svalbard using drone data. The results point out that the range of z0 values across a crevassed glacier is large but in general significantly higher compared to non-crevassed glacier surfaces. The UAV approach proved to be an ideal tool to provide distributed z0 estimates of crevassed glaciers which can be used to model turbulent fluxes.