Articles | Volume 15, issue 6
Research article
 | Highlight paper
11 Jun 2021
Research article | Highlight paper |  | 11 Jun 2021

Mapping the aerodynamic roughness of the Greenland Ice Sheet surface using ICESat-2: evaluation over the K-transect

Maurice van Tiggelen, Paul C. J. P. Smeets, Carleen H. Reijmer, Bert Wouters, Jakob F. Steiner, Emile J. Nieuwstraten, Walter W. Immerzeel, and Michiel R. van den Broeke

The aerodynamic roughness of heat, moisture, and momentum of a natural surface are important parameters in atmospheric models, as they co-determine the intensity of turbulent transfer between the atmosphere and the surface. Unfortunately this parameter is often poorly known, especially in remote areas where neither high-resolution elevation models nor eddy-covariance measurements are available. In this study we adapt a bulk drag partitioning model to estimate the aerodynamic roughness length (z0m) such that it can be applied to 1D (i.e. unidirectional) elevation profiles, typically measured by laser altimeters. We apply the model to a rough ice surface on the K-transect (west Greenland Ice Sheet) using UAV photogrammetry, and we evaluate the modelled roughness against in situ eddy-covariance observations. We then present a method to estimate the topography at 1 m horizontal resolution using the ICESat-2 satellite laser altimeter, and we demonstrate the high precision of the satellite elevation profiles against UAV photogrammetry. The currently available satellite profiles are used to map the aerodynamic roughness during different time periods along the K-transect, that is compared to an extensive dataset of in situ observations. We find a considerable spatio-temporal variability in z0m, ranging between 10−4 m for a smooth snow surface and 10−1 m for rough crevassed areas, which confirms the need to incorporate a variable aerodynamic roughness in atmospheric models over ice sheets.

1 Introduction

Between 1992 and 2018, the mass loss of the Greenland Ice Sheet (GrIS) contributed 10.8 ± 0.9 mm to global mean sea-level rise (Shepherd et al.2020). This mass loss is caused in approximately equal parts by an increase in ice discharge (Mouginot et al.2019; King et al.2020) and an increase in surface meltwater runoff (Noël et al.2019). Runoff occurs mostly in the low-lying ablation area of the GrIS, where bare ice is exposed to on-average positive air temperatures throughout summer (Smeets et al.2018; Fausto et al.2021). As a consequence, the downward turbulent mixing of warmer air towards the bare ice, the sensible heat flux, is an important driver of GrIS mass loss next to radiative fluxes (Fausto et al.2016; Kuipers Munneke et al.2018; Van Tiggelen et al.2020).

Although the strong vertical temperature gradient provides the required source of energy, it is the persistent katabatic winds that generate the turbulent mixing through wind shear (Forrer and Rotach1997; Heinemann1999). Additionally, the surface of the GrIS close to the ice edge is very rough (Yi et al.2005; Smeets and Van den Broeke2008). It is composed of closely spaced obstacles, such as ice hummocks, crevasses, melt streams, and moulins. Due to the effect of form drag (or pressure drag) τr, the magnitude of the turbulent fluxes increases with surface roughness (e.g. Garratt1992), thereby enhancing surface melt (Van den Broeke1996; Herzfeld et al.2006). As of today, the effect of form drag on the sensible heat flux over the GrIS, and therefore its impact on surface runoff, remains poorly known.

The first challenge in modelling this turbulent mixing resides in accurately modelling the surface shear stress, without the need to calculate the detailed air pressure distribution around each individual surface obstacle. Such bulk drag models have been developed by, e.g., Arya (1975) to estimate the drag caused by pressure ridges on Arctic pack ice. This model was extended by Hanssen-Bauer and Gjessing (1988) for varying sea-ice concentrations. A more general drag model was proposed by Raupach (1992), which was extended by Andreas (1995) for sastrugi, by Smeets et al. (1999) for rough ice, and by Shao and Yang (2008) for surfaces with higher obstacle density, such as urban areas. Lüpkes et al. (2012) and Lüpkes and Gryanik (2015) developed a bulk drag model for sea ice that is used in multiple atmospheric models. Over glaciers, semi-empirical approaches based on Lettau (1969) are often used, such as by Munro (1989), Fitzpatrick et al. (2019), and Chambers et al. (2019).

The second challenge is the application of such models in weather and climate models, which requires mapping small-scale obstacles over large areas, e.g. an entire glacier or ice sheet. Historically, the surveying of rough ice was spatially limited to areas accessible for instrument deployment, possibly introducing a bias when it comes to quantifying the overall roughness of a glacier. The recent development of airborne techniques, such as uncrewed aerial vehicle (UAV) photogrammetry and airborne lidar, opened up new possibilities for mapping surface roughness properties. While these techniques enable the high-resolution mapping of roughness obstacles, they often only cover portions of a glacier or ice sheet. On the other hand, satellite altimetry provides the means to cover entire ice sheets, though the horizontal resolution remains a limiting factor when mapping all the obstacles that contribute to form drag. Depending on the type of surface, parameterizations using available satellite products are possible, as presented for Arctic sea ice by Lüpkes et al. (2013), Petty et al. (2017), and Nolin and Mar (2019).

The third and final challenge is the experimental validation of bulk drag models over remote rough ice areas, which requires either in situ eddy-covariance or multi-level wind and temperature measurements. Long-term and continuous datasets remain scarce on the GrIS, although simplifying in situ methods can be applied for long-term monitoring of turbulent fluxes (Van Tiggelen et al.2020, T20).

In this paper, we address the first two challenges by applying the model of Raupach (1992) to 1 m resolution elevation profiles measured over the west GrIS by the ICESat-2 laser altimeter. We apply the bulk drag model to roughness information from UAV photogrammetry, and we address the third challenge by evaluating the modelled aerodynamic roughness against in situ eddy-covariance measurements. We then evaluate the ICESat-2 elevation profiles against UAV photogrammetry, and we finally apply the bulk drag model to the ICESat-2 profiles obtained over an extended area and during different time periods.

This paper is organized as follows. In Sect. 2 we describe the modifications in the bulk drag model, and in Sect. 3 we describe the elevation datasets used to force the model. We then evaluate the bulk drag model for one site in Sect. 4.1 and the roughness statistics derived from ICESat-2 at multiple sites in Sect. 4.3. In Sect. 4.4 we apply the model to map the aerodynamic roughness length (z0m) along the K-transect, on the west GrIS.

Figure 1(a) Map of the K-transect, with the location of the automatic weather stations and mass balance sites indicated by the pink diamonds. The black boxes A and B delineate the areas mapped by UAV photogrammetry. The large black box indicates the area covered in Figs. 6 and 9. The background image was taken by the MSI instrument (ESA, Sentinel-2) on 12 August 2019. Pixel intensity is manually adjusted over the ice sheet for increased contrast. The green solid lines denote the ICESat-2 laser tracks that are compared to the UAV surveys (Table 2). (b) Sites S5 (6 September 2019), S6 (6 September 2019), and S9 (3 September 2019) taken during the yearly maintenance. Note that no data from the AWS shown at S9 are used in this study. (c) Location of the K-transect on the Greenland Ice Sheet.

2 Model

2.1 Definition of the aerodynamic roughness length z0m

Atmospheric models assume that the lowest grid point above the surface is located in the inertial sublayer (or surface layer). In this layer, the eddy diffusivity for momentum increases linearly with height and decreases with atmospheric stability, which yields the semi-logarithmic vertical profile of horizontal wind speed. Over a rough surface, the pressure drag force on the obstacles acts as an additional sink of momentum, next to skin friction. Furthermore, the turbulent wakes generated by the flow separation enhance turbulent mixing. This may be approximated by an increase in the eddy diffusivity in the roughness sublayer (Garratt1992; Raupach1992; Harman and Finnigan2007). As such, the vertical profile of horizontal wind speed (u(z)) over a rough surface can be written as

(1) u ( z ) = u * κ ln z - d z 0m - Ψ m z - d L o + Ψ m z 0m L o + Ψ m ^ ( z ) ,

where z is the height above the surface, u*=τρ0.5 is the friction velocity, ρ the air density, κ=0.4 the von Kármán constant, τ the total surface shear stress, and z0m the roughness length for momentum. The average wind profile in Eq. (1) is shifted upwards by a displacement height d, which is defined as the centroid of the drag force profile on the roughness elements (Jackson1981). z0m is thus defined as the height above d, where u(z)=0. The dependency of the eddy diffusivity for momentum on the diabatic stability and on the turbulent wake diffusion are described as Ψmz-dLo and Ψm(z)^, respectively, where Lo is the Obukhov length. The hat notation is used for the roughness layer quantities, as in Harman and Finnigan (2007). Above the roughness sublayer, Ψm(z)^=0.

The problem we address is the estimation of z0m. Rewriting Eq. (1) and assuming neutral conditions (i.e. Ψm=0), yields

(2) z 0m = ( z - d ) exp κ u ( z ) u * - Ψ m ( z ) ^ ) - 1 .

Hence, the process of finding z0m is equivalent to finding d, u(z)u*, and Ψm(z)^ simultaneously.

2.2 Bulk drag model of z0m

The main task is to model the total surface shear stress τ=ρu*2, which for a rough surface is the sum of both form drag τr and skin friction τs:

(3) τ = τ r + τ s .

Both τr and τs are parameters of the flow but can be related to the geometry of the roughness obstacles using a bulk drag model. Two important parameters of the roughness obstacles are their height (H) and their frontal area index (λ), defined as

(4) λ = A f A l ,

with Af the frontal area of the roughness obstacles perpendicular to the flow and Al the total horizontal area. At this point, we will differ from the model by Shao and Yang (2008), who add an extra term in Eq. (3) in order to separate the skin friction at the roughness elements and the underlying surface. We also differ from the models by Lüpkes et al. (2012) and Lüpkes and Gryanik (2015), where skin friction over sea ice is separated between a component over open water and a component over ice floes. In the case of a rough ice surface, there is no clear distinction between the obstacles and the underlying surface. Therefore, we follow the model of Raupach (1992, R92), which is designed for surfaces with a moderate frontal area index (λ<0.2). As we will see in the next sections, the ice surfaces considered here do not exceed λ=0.2. As a comparison we will also consider the models from Lettau (1969, L69) and Macdonald et al. (1998, M98). The detailed equations of the bulk drag models can be found in Appendix A.

2.3 Definition of the height (H) and frontal area index (λ) over a rough ice surface

Here we introduce the type of surfaces that we are considering. Our aim is to model the aerodynamic roughness of a rough ice surface, including its dependence on wind direction (Van Tiggelen et al.2020). We will consider rectangular elevation profiles of length L=200 m, measured upwind from a point of interest (e.g. an automatic weather station, or AWS). This geometry is a strong simplification of the true fetch footprint, which is calculated for a specific wind direction at S5 in Fig. 2, after Kljun et al. (2015). Yet this simplification allows us to use 1D elevation datasets, such as profiles from the ICESat-2 satellite laser altimeter. Besides, the true fetch footprint depends on flow parameters such as the friction velocity (u*) and the boundary-layer height (Kljun et al.2015), which are not known a priori.

Figure 2(a) Measured elevation profiles for four different wind directions upwind of AWS S5. (b) Filtered elevation profiles and (c) orthomosaic true-colour image of AWS S5 and surroundings taken by UAV photogrammetry on 6 September 2019. The different coloured rectangles in (c) indicate the profiles shown in panel (a). The profiles have been vertically offset by 5 m in (a) and by 2 m in (b) for clarity. The black line in (a) denotes the low-frequency contribution of the profiles for a cut-off wavelength Λ=35 m. The pink arrow in (c) denotes the displacement vector of the AWS between the ICESat-2 overpass on 14 March 2019 and the UAV imagery on 6 September 2019. The estimated extent of the 50 % and 80 % fetch footprints for the data in September 2019 in a specific wind direction is shown by the black ovals.

Four measured elevation profiles and a high-resolution orthomosaic image are shown in Fig. 2. These were measured on 6 September 2019 at site S5 (67.094 N, 50.069 W, 560 m) in the locally prevailing wind directions, using UAV photogrammetry, of which the details will be given in Sect. 3. At this site, pyramidal ice hummocks with heights between 0.5 and 1.5 m are superimposed on larger domes of more than 50 m in diameter (see also Fig. 1b). The elevation profiles for different fetch directions illustrate three important issues: (1) the zero-referencing of the surface, (2) the identification of distinct roughness obstacles, and (3) the important variability of the surrounding topography, depending on the fetch direction. The obstacles being anisotropic, the surface appears rougher in the southerly directions than in the easterly directions. In addition, the ice ridges and troughs have variable heights and depths, which means that describing this rough ice surface with a few length scales (e.g. H, λ) in order to estimate the aerodynamic roughness will introduce some uncertainty. This is mainly because each individual obstacle has a different contribution to the total drag. Unfortunately, these individual drag contributions cannot be modelled, due to the unknown shape of the wind profile between the roughness elements. Without a universal theory of drag over complex surfaces, several simplifications need to be made.

We choose here to approximate the true surface as an array of f identical obstacles of height H in the profile of length L (Munro1989; Smith et al.2016) (Fig. 3). This avoids the use of empirical formulas for the estimation of z0m and allows us to apply the bulk drag models. The approach of approximating a natural surface by uniquely shaped obstacles is formally justified by Kean and Smith (2006), as most of the form drag is caused by the largest and steepest obstacles. On the other hand, large natural obstacles also tend to be wide, so their relatively small frontal area index considerably reduces their contribution to the total form drag (Fitzpatrick et al.2019). To remove the influence of the widest obstacles, the elevation profile of length L is linearly detrended, and the power spectral density of the detrended profile is computed in order to filter out all the wavelengths larger than the cutoff wavelength Λ=35 m. This value is found to give optimal results, which is shown in Appendix B. In order to avoid spectral leakage when applying Fourier statistics on short and aperiodic signals, we extend each input profile with the identical but mirrored profile before computing the power spectral density. This yields a symmetrical and thus periodic profile of length 2L, which is then high-pass filtered. The final statistics are then computed using the first half of the filtered profile of length L. Typical filtered profiles are shown in Figs. 2b and 3. These profiles only contain the ice hummocks, as the high-pass filter removed the influence of the large-scale domes.

Figure 3Steps in converting a measured digital elevation model to the modelled topography, where L is the length of the profile, f the number of obstacles, H the height of the obstacles, and w the width of the elevation profile. The location and height of AWS S5 are shown on top of the UAV elevation profile. The grey dots denote all the ATL03 photons, while the black dots denote the selected photons for the kriging procedure. The solid blue line denotes the 1 m resolution interpolated profile for ATL03 data, and the pink dots denote the 20 m resolution ATL06 signal.


The height of the roughness obstacles (H) is taken as

(5) H = 2 σ z ̃ ,

in which σz̃ is the standard deviation of the filtered elevation profile. This is an arbitrary but convenient choice, as the standard deviation of the topography captures all the scales in the filtered profile but remains insensitive to the height of the small-scale obstacles, which we assume to have a negligible influence on the overall drag. Unfortunately, the variance is sensitive to the height of the largest obstacles and thus to the chosen value for Λ.

Next, we define an obstacle as a group of consecutive positive values of filtered heights, after Munro (1989) (see also Smith et al.2016), which yields f, the number of obstacles. The obstacle frontal area index (λ) in the direction of the elevation profile is then computed as (Fig. 3)

(6) λ = f H w L w = f H L ,

where w is the width of the profile, set to 15 m. This value was chosen to match the approximate ICESat-2 footprint diameter, yet it is much smaller than the width of the real fetch footprint (Fig. 2). We assume that the obstacles and the elevation profile have the same width, which removes all information about the shape of the obstacles in the direction perpendicular to the wind direction. This simplification avoids the additional uncertainty regarding the aggregation of 2D datasets in the process of modelling z0m and allows us to apply the model to ICESat-2 profiles.

To summarize, a measured elevation profile is now completely defined by the height of the obstacles (H), and the frontal area index (λ), after high-pass filtering (see Fig. 3). This now allows us to apply a bulk drag model to estimate one value for z0m per 200 m profile. The exact placement of the obstacles is resolved in the process (Fig. 3) but does not serve as input for the drag models. Detailed equations of the bulk drag model can be found in Appendix A.

Table 1Description of z0m in situ datasets on the K-transect. The data from Meesters et al. (1997), Smeets and Van den Broeke (2008), Lenaerts et al. (2014), and Van Tiggelen et al. (2020) are denoted M97, SB08, L14, and T20 , respectively. The measurement methods are (p) profile, (sec) sonic eddy covariance, or (vpec) vertical propeller eddy covariance.

Download Print Version | Download XLSX

3 Datasets

3.1 Eddy-covariance measurements

Vertical propeller eddy-covariance (VPEC; see also T20) measurements are available at sites S5 (67.094 N, 50.069 W, 560 m) and S6 (67.079 N, 49.407 W, 1010 m) since 2016, while AWS observations are available since 1993 and 1995 for each site (Smeets et al.2018). For this study we use eddy-covariance measurements acquired during September 2019 at site S5 and also site SHR (67.097 N, 49.957 W, 710 m), and from September 2018 to August 2019 at site S6. All these sites are situated in the lower ablation area of the K-transect, which is a 140 km transect of AWS and mass balance observations on the western part of the GrIS (van de Wal et al.2012; Smeets et al.2018). It extends from the ice edge up to 1850 m elevation and therefore covers many contrasting types of surfaces, ranging from the rough crevassed bare ice close to the ice edge, to the year-round firn-covered surface at the highest locations (see Figs. 1 and 6). At the end of the melting season, the bare ice surface at S5 and SHR is characterized by densely packed hummocks up to 1.5 m height, while at S6 it is characterized by more sparsely packed hummocks of 0.6 m average height. The datasets include 30 min observations of the friction velocity u*(z) and wind speed u(z) at the same height above the surface (z=3.7 m). Two independent techniques were used at S5 and SHR. The first technique is the sonic eddy-covariance (SEC) method, which uses measurements from a sonic anemometer (CSAT3B, Campbell scientific, Logan, USA) sampled at 10 Hz. The second technique is the VPEC method, which relies on measurements of a vertical propeller, horizontal propeller, and fine-wire thermocouple sampled at 5 Hz. At S6, only the VPEC method was used with a sampling interval of 4 s. For both methods, the roughness length (z0m) is calculated using Eq. (2).

We only select data taken during near-neutral conditions (z/Lo<0.1), and we assume that the measurements are taken above the roughness layer, i.e. Ψm(z)^=0. The latter is a reasonable assumption, given that the height of the obstacles (H) at these sites is less than 1.5 m, which means that the roughness layer unlikely exceeds 3 m (Smeets et al.1999; Harman and Finnigan2007). On the other hand, when applying the drag model to estimate z0m (Appendix A), the correction factor Ψm(z)^ is taken into account. The reason is that the obstacles are located in the roughness layer, where the vertical wind profiles deviate from the inertial sublayer wind profiles, according to Eq. (1). Details about the processing steps and further data selection strategies can be found in T20. The data selection strategy removes all data points with wind directions outside the [80;200] interval. In the following sections we average ln(z0m), which is of interest for the determination of the vertical profile of horizontal wind speed (Eq. 1). Additional in situ averaged z0m measurements obtained during different time periods and at several locations along the K-transect are taken from Meesters et al. (1997, M97), Smeets and Van den Broeke (2008, SB08), Lenaerts et al. (2014, L14), and T20 and are summarized in Table 1.

3.2 UAV structure from motion

The high-resolution elevation maps are derived using a structure-from-motion workflow using UAV imagery. Two crevassed areas close to the ice edge were mapped using an eBee fixed-wing UAV from Sensefly, while the area surrounding S5 was mapped using a Mavic Pro quadcopter UAV from DJI. Multiple overlapping true-colour images of the surface are processed in Agisoft Photoscan to produce 3D elevation maps. Detailed information about this workflow can be found in Immerzeel et al. (2014), Kraaijenbrink et al. (2016), and references therein. Briefly, the same surface features are identified on different images and are used to reconstruct the 3D geometry between the surface and the camera position. The resulting point cloud of the surface is then gridded and finally geo-referenced using the information of the UAV GPS, which yields a digital elevation model (DEM) of the surface. No additional ground-control points were used for the elevation maps, which is of little relevance in this study, as we are not interested in the exact absolute elevation but in relative obstacle heights. Details about the UAV DEMs are provided in Table 2.

Table 2Description of DEMs obtained by UAV photogrammetry and of the corresponding overlapping ICESat-2 laser beams.

Download Print Version | Download XLSX

The elevation profiles are then extracted by projecting all the DEM points in a 200 m × 15 m rectangle on the centre line, followed by averaging the projected points in 1 m bins (see Fig. 3). The aim of this averaging method is to mimic ICESat-2 profiles.

3.3 ICESat-2 laser altimeter

Launched in September 2018 by the National Aeronautics and Space Administration (NASA), ICESat-2 (Ice, Cloud, and Elevation Satellite-2) carries a laser altimeter system in near-polar orbit (Markus et al.2017). The altimeter relies on a photon-counting system, which in combination with both the spacecraft's position and its pointing orientation, enables the retrieval of the 3-D position of individual backscattered photons (Neumann et al.2019). Our hypothesis is that the small footprint diameter ( 15 m) and short along-track spacing between these footprints (0.7 m) allows for an accurate estimation of land ice aerodynamic roughness properties.

A typical geolocated photon measurement ATL03 (Neumann et al.2019) can be seen in Fig. 3 for site S5 and in Fig. 4a for area A. Details about which ICESat-2 measurements are compared against the UAV surveys are provided in Table 2. Not more than one ICESat-2 measurement exactly overlaps each UAV survey. This is mainly due to the presence of clouds and due to changes in laser pointing orientations in other ICESat-2 measurements, but also due to changes in the studied locations due to ice flow. The global geolocated photon product ATL03 requires some processing steps before the roughness statistics can be computed. These steps mainly involve the selection of valid photons, aggregating the 3-D photon positions on a regular along-track grid, and finally correcting for remaining biases. The standard algorithm used to derive an accurate estimate of land ice height product from ATL03 to ATL06 is described in detail by Smith et al. (2019). Unfortunately the 20 m along-track resolution of the ATL06 land ice height product is too coarse for aerodynamic roughness calculations for two reasons. First, in the ATL06 product there are only 10 points in 200 m sections, which is not enough to apply the high-pass filter. And second, on this scale the number of roughness obstacles (f) would be greatly underestimated as can be seen in Figs. 3 and 4. Fortunately, information smaller than the footprint diameter can be extracted from the ATL03 product, as shown by Herzfeld et al. (2020), in which a density–dimension algorithm is used that facilitates surface height determination at the 0.7 m nominal along-tack resolution. In the following part we describe a method to produce a 1 m resolution along-track surface height estimation from the ATL03 raw photons signal.

Figure 4(a) Elevation profile at site A measured by the UAV and by ICESat-2 (solid lines), selected ICESat-2 photons (black dots), and ICESat-2 ATL06 height (pink dashed line). The UAV and ATL06 profiles have been vertically offset by 2 m for clarity. (b) Filtered profiles (solid lines) and residual photon elevations after filtering per 200 m windows (grey dots), where the UAV and ATL03 filtered profiles have also been vertically offset. (c) Probability density function of the filtered ICESat-2 profile (blue dashed line), UAV profile (orange solid line), and residual photon elevations (black line).


The first step involves selecting all the ATL03 photons that have been flagged as either low, medium, or high confidence by the ATL03 algorithm. All the selected photons are projected on the along-track segment, and a median absolute difference filter is used to remove all the photon heights which deviate too much from the local ensemble median,

(7) z - q low 0.6745 | z - z | z z + q high 0.6745 | z - z | ,

where z denotes the median of z within a moving window. We choose qlow=1 and qhigh=2 in order to filter more photons below than above the median. We assume that the highest detected photons are more likely to be first surface reflections, while the lower photons are more likely to be delayed by scattering. We set the window length to 50 m. The previous selection strategy could also be applied for retrieving the surface in the case of multiple reflections (e.g. shallow supraglacial lakes), but this was not tested.

The second step involves interpolating the irregular photon locations on a regular, 1 m resolution, along-track grid. The overlap between the individual footprints means that the geolocated photon heights in close vicinity must be correlated to each other, with a correlation diameter similar to twice the footprint diameter ( 30 m). We take advantage of this feature to interpolate the ATL03 photons using a k-nearest neighbour, one-dimensional, ordinary kriging algorithm, of which details can be found in Hengl (2009). In essence, the interpolation weights depend on the covariance with the nearby measurements, which is assumed to decrease over distance. A Gaussian covariance function with a radius of 15 m is found to fit the experimental semi-variograms best. For computational efficiency, only the 100 closest geolocated photons within a quarter footprint diameter (3.75 m) of each grid point are used for the interpolation. We only choose the high confidence photons, but if there is less than 1 photon per 0.7 m, we also select the medium-confidence photons. If there are not enough medium-confidence photons, we increase the search radius to half the footprint diameter (7.5 m) or even up to a footprint diameter (15 m). The low-confidence photons are only used as a last resort. If in a 15 m footprint diameter there are still not enough photons present, the height on that grid point is not estimated, which results in a gap. A sensitivity experiment using different photon selection strategies and different kriging parameters is found in Appendix B.

The last step involves grouping the interpolated elevation measurements in 200 m along-track windows and the high-pass filtering using a cut-off wavelength of Λ=35 m (Sect. 2). The height of the obstacles (H) is defined as twice the standard deviation of the filtered signal (Eq. 5).

Although 1 m resolution is still too coarse to capture all the small-scale obstacles that contribute to form drag, we expect that most of the form drag over rough ice is caused by the larger obstacles that are resolved by the ICESat-2 altimeter. Furthermore, the small-scale information is still indirectly present in the scatter of the surrounding photons to the closest grid point, which is a measure of both the instrumental error and the surface slope, but also of the surface roughness (Gardner1982).

An alternative approach that does not require gridding the ATL03 product to 1 m resolution would be to use the standard deviation of the raw photon signal detrended for the resolved 20 m resolution ATL06 data, as in Yi et al. (2005) and Kurtz et al. (2008). However, as we will see in the following sections, this would overestimate the height of the roughness obstacles. In addition, the frontal area index (λ) would remain unknown.

When working with the 1 m interpolation profile, we model the standard deviation of the unresolved topography (σsub) according to

(8) σ sub = σ ph,res 2 - σ i 2 0.5 / 2 ,

where σph,res is the standard deviation of the photon residual elevations, defined as the signal of the selected photons minus the interpolated 1 m resolution profile (Fig. 4), σi=0.13 m is the standard deviation due to the instrumental precision (Brunt et al.2019). We calculate σsub for each 200 m profile. The total variance of the surface elevation measured by the laser altimeter in 200 m intervals is the sum of both the resolved and unresolved variance:

(9) σ tot = σ res ̃ 2 + σ sub 2 ,

in which σres̃ is the resolved standard deviation of the filtered 1 m resolution profile. The height of the roughness obstacles, corrected for the unresolved topography, is then estimated according to

(10) H corr = 2 σ tot .

The obstacle frontal area index (λ) is finally computed using Eq. (6), where the number of obstacles (f) is estimated from the filtered profiles. Both H and λ are then used as input for the bulk drag model (Appendix A), which results in one value for z0m per 200 m profile.

The filtered ICESat-2 signal and residual photon elevations at site A are shown in Fig. 4b, and their probability density functions are shown in Fig. 4c. At this site, the filtered ICESat-2 signal at 1 m resolution captures most of the information present in the UAV signal. On the other hand, the residual photon elevations, defined as the selected photons detrended for the interpolated profile under Eq. (9), still contain much larger scatter than the UAV elevation profile. This demonstrates that roughness is not the only factor explaining the scatter in the raw altimeter signal. Therefore using the residual scatter (Eq. 9) will overestimate the height of the roughness obstacles. In the next sections, we will analyse the uncorrected height of the obstacles (H), unless stated otherwise.

4 Results

4.1 Evaluation of the bulk drag model forced with a UAV DEM

Bulk drag models are often used as a convenient way to estimate the aerodynamic properties of natural surfaces. Nevertheless, the number of quantitative evaluations of these models for rough snow and ice surfaces is very limited. Brock et al. (2006) found that z0m modelled using the method by L69 (Eq. A13) agrees well with observations over melting ice on a mountain glacier, although they used shorter profiles, up to 15 m in length, and sampled in the orientation perpendicular to the wind direction. On the other hand, Van den Broeke (1996) found that L69 overestimates z0m at site S4 at the K-transect (lowest site in Fig. 6). The same overestimation was found by Smeets et al. (1999), Fitzpatrick et al. (2019), and Chambers et al. (2019) for rough glacier ice but also by Miles et al. (2017) for a debris-covered glacier. These studies all use different methods at different sites to estimate H and λ, which illustrates the limited suitability of the model by L69 for realistic snow and ice surfaces.

Figure 5Modelled z0m at site S5 using three different bulk drag models, Lettau (1969, L69, blue lines), Macdonald et al. (1998, M98, green lines), and Raupach (1992, R92, red lines), and using two different values for the drag coefficient for form drag: Cd=0.25 (solid lines) and Cd=0.1 (dashed lines). Solid grey symbols are measurements from sonic eddy covariance (SEC) or vertical propeller eddy covariance (VPEC). Additional data are from Van Tiggelen et al. (2020, T20). Pink circles are the model results forced with H and λ from UAV photogrammetry, using the R92 model and Cd parameterized using Eq. (A2).


To verify the suitability of several drag models (see Appendix A), we use the eddy-covariance observations at site S5 as independent validation (Sect. 3). Different values of z0m are calculated for different fetch directions as depicted in Fig. 2. Figure 5 compares both the estimated z0m from in situ observations and the modelled z0m at the end of the ablation season, as a function of the measured obstacle frontal area index λ. The L69 model (Eq. A13) overestimates z0m for λ<0.04 at this location (Fig. 5, blue line). In accordance with L69, the drag coefficient of an individual obstacle Cd=0.25 is likely too high for naturally streamlined obstacles. Furthermore, L69 does not consider the displacement height, which means that the height of the obstacles (H) relevant for form drag is overestimated. Nevertheless, L69 still yields a reasonable estimate of z0m for λ>0.04, which can be explained by the neglect of the displacement that is compensated for by too small Cd for these fetch directions. The method by M98 (Eq. A14) does account for the displacement height, and, while using the same drag coefficient Cd=0.25, it gives improved results for λ<0.04 compared to L69 (Fig. 5, green line). The same holds for the model by R92 (Fig. 5, red line). M98 is expected to fail for very small λ, due to the absence of skin friction. Using Cd=0.1, all three models perform better for λ<0.05 but perform poorly for λ>0.04 (Fig. 5, dashed lines). This is a strong indication that Cd is not constant, but varies with the wind direction, depending on the exact placement and shape of the obstacles. In Sect. 4.3 we estimate the values for Cd required to fit the model to the observations; these values vary between 0.1 and 0.3 and show a weak relationship with H. The parametrization for Cd from Garbrecht et al. (1999) (Eq. A2), for which Cd increases with H, yields most acceptable results when used in combination with the R92 model (Fig. 5). Note that Lüpkes et al. (2012) use a constant value for Cd.

Figure 6Estimated height of the roughness obstacles (H) from ICESat-2 between 16 October 2018 and 6 September 2020 in the lower part of the K-transect, West Greenland. The locations of the automatic weather stations are given by the pink diamonds. The black boxes A and B delineate the areas mapped by UAV photogrammetry. A hillshade of ArcticDEM (Porter et al.2018) is shown as background over the ice sheet.

The R92 model with the parametrization for Cd allows for some variability in modelled z0m for the same λ but is still not able to reproduce the eddy-covariance observations (Fig. 5). We attribute this to the parametrization of Cd, which was derived for sea-ice pressure ridges and therefore likely less suitable for rough ice hummocks. Nevertheless, the overall error between model and observation is acceptable, given the simplicity of the bulk drag models that were designed for idealized roughness geometries. As pointed out by L69, realistic modelling of total drag over a complex natural surface should intuitively require a complete variance spectrum of the topography. Linking variance spectra to the total drag has been investigated recently through numerical simulations (Yang et al.2016; Zhu and Anderson2019; Li et al.2020), but a universal and physically based relationship for complex surfaces is still lacking. In the next sections, we therefore use the model of R92 with a parametrized Cd for mapping z0m using either UAV or ICESat-2 profiles.

4.2 Height of the roughness obstacles (H) estimated from ICESat-2

The estimated height of the obstacles (H) using 2 years of ICESat-2 measurements (16 October 2018–6 September 2020) crossing the lower part of the K-transect is shown in Fig. 6. H ranges between less than 0.1 m at the higher locations and more than 3 m in rough crevassed areas near the ice edge. At first glance a clear pattern of roughness emerges, in which ice dynamics and elevation seem to be the controlling factors. Low-lying bare-ice areas are rougher, while the higher, firn-covered areas are smooth. Nevertheless, the roughness is very variable locally due to isolated crevasses and melt channels. In addition, we expect a seasonal variability that is not yet captured in this analysis.

4.3 Evaluation of ICESat-2 roughness statistics against UAV DEMs

Climate models and satellite altimeter corrections require information about the larger-scale spatial variability of surface (aerodynamic) roughness. This motivated us to compare the roughness statistics acquired with high-resolution UAV photogrammetry to the statistics estimated from the ICESat-2 laser altimeter.

Figure 7(a1, a2) Estimated height of the roughness obstacles (H). (b1, b2) Estimated frontal area index. (c1, c2) Estimated aerodynamic roughness length (z0m). (a1, b1, c1) Area A. (a2, b2 c2) Area B. The black lines denote the roughness statistics estimated from the ATL03 filtered profile, with or without accounting for the residual photon elevations (dashed and solid lines, respectively). The orange line denotes estimates using UAV elevation profiles, and the pink line denotes the height of obstacles estimated using the scatter of ATL03 photons detrended for the ATL06 signal.


The elevation profile from the UAV survey in box A (Fig. 6) was already compared to the overlapping ICESat-2 profiles in Fig. 4a, while H, λ, and z0m are compared in Fig. 7. In box A, the UAV and ICESat-2 profiles were taken 11 d apart at the end of the ablation season. The height (H) and frontal area index (λ) of the roughness obstacles are estimated for 200 m intervals, with each interval centre separated by 50 m. Overall, the uncorrected 1 m profile from ICESat-2 (Fig. 7, solid black line) clearly captures all the largest obstacles and the large-scale variability but still slightly underestimates both the height (H) and the frontal area index (λ) of the obstacles, compared to the UAV surveys (Fig. 7, orange line). This is expected, given the size of the laser footprint and the low-pass-filtering properties of the kriging procedure. On the other hand, the correction using the standard deviation of the photon distribution (Eq. 9) overestimates H (Fig. 7, dashed black line). This can be explained by additional processes that affect the local photon distribution but that we did not consider, such as the forward scattering in the atmosphere (Kurtz et al.2008), the penetration of photons in the ice layer (Cooper et al.2021), or simply the presence of outliers that passed the median absolute difference filter (Eq. 7). Furthermore, the obstacle frontal area index (λ) is underestimated by the ICESat-2 altimeter, since we do not account for unresolved obstacles when counting the number of obstacles (f). In addition, using the standard deviations of the ATL03 product detrended for the 20 m resolution ATL06 signal results in an even greater overestimation of H (Fig. 7, purple line). This is due to the fact that besides the additional processes broadening the altimeter signal, the scatter of this signal also contains the large-scale variability at wavelengths larger than Λ=35 m. We assumed that such large wavelengths can be neglected in the drag calculations; therefore they are removed in the filtered UAV and ICESat-2 profiles.

Two more UAV surveys were performed in September 2019 in area B and around S5, but the overlapping ICESat-2 profiles were measured during winter (see Table 2). The comparison of H, λ, and modelled z0m is given in Fig. 7. In area B, crevassed and slightly rougher than A, the elevation was measured in December, 3 months after the UAV survey. The uncorrected ICESat-2 profiles show a slightly more pronounced underestimation of H compared to area A, which we relate to snowfall reducing the height of the roughness obstacles. On the other hand, the corrected ICESat-2 profiles overestimate H by 0.06 m, which translates into an overestimation of z0m by approximately 2.5×10-3 m (Fig. 7). On average, the uncorrected ICESat-2 values underestimate z0m by 2.9×10-3 m for area A and 9×10-3 m for area B, which corresponds to  40 % and  36 % of the average z0m estimated by the UAV at these two sites.

Figure 8(a) Drag model evaluation at site S5. (b) Drag coefficient for form drag (Cd) used in the bulk drag model (black line) or required to perfectly fit the observations. The solid orange line is the modelled z0m using the R92 model and UAV photogrammetry on 6 September 2019, while the dashed orange line is the orange line shifted down by a factor of 10. Solid symbols are measurements from sonic eddy-covariance (SEC) or vertical propeller eddy-covariance (VPEC). Additional data are from Van Tiggelen et al. (2020, T20). The vertical dashed line denotes the direction sampled by the ICESat-2 laser beam on 14 March 2019. The error bar denotes the range between the uncorrected and corrected ICESat-2 measurements.


At site S5, UAV elevation profiles and eddy-covariance measurements are available in September 2019, while the ICESat-2 elevation profile was measured in March (Table 2). Both the satellite and the UAV profiles are shown in Fig. 3. Although the UAV profile is too short to statistically compare H and λ to the ICESat-2 altimeter, the qualitative comparison between the two confirms that the satellite altimeter is very capable of detecting most of the obstacles that are smaller than 20 m in width. Interestingly, some depressions in the UAV DEM are not captured by ICESat-2, most likely as a result of snow filling them in March. Furthermore, the bending (or “doming”) of the UAV profile is visible near the edges, which is a consequence of the lack of reference ground control points in the UAV data processing, which is a common issue with UAV data processing (James and Robson2014). Both H and λ are smaller in the satellite profile than in the UAV profile, but the modelled z0m agrees qualitatively with that estimated using observations from the AWS S5 during March–April. During this time period, z0m is approximately a factor of 10 smaller than during the end of the ablation season (Fig. 8, dashed orange line). Unfortunately, the track direction of the satellite altimeter rarely coincides with the wind direction measured by the anemometers at this location, due to the katabatic forcing. This prevents a direct comparison of ICESat-2 roughness to in situ observations, as the aerodynamic roughness strongly depends on the wind direction (Van Tiggelen et al.2020). The z0m value estimated from ICESat-2 profiles must thus be interpreted as the aerodynamic roughness in the wind direction along the direction of the ground laser track.

Only a high-resolution, two-dimensional DEM, e.g. obtained using a UAV, allows for an accurate description of the aerodynamic roughness around a point of interest in multiple directions. An example of such an analysis is shown for site S5 in Fig. 8. The R92 model applied to the UAV elevation profiles reproduces the considerable variability of the estimated z0m using in situ observations. Across all the wind directions available in the measurements, z0m using the UAV profiles is underestimated by 7.6×10-3 m, or 28 % of the average value estimated by the SEC method in September 2019 (2.65×10-2 m). The comparison improves when comparing the modelled z0m to VPEC measurements from September 2016 and August 2017 (T20), the model now overestimates the estimated z0m by in situ observations by 1.1×10-3 m or 9 % of the observed value (1.25×10-2 m). As these data contain more wind directions, the overestimation of z0m in the southerly fetch directions is compensated for by an underestimation in the easterly directions (Fig. 8). The difference between different in situ data highlights the variability in z0m in time but also the uncertainty in the field measurements. The difference in averaged estimated z0m using in situ observations during the overlapping period across all wind directions is 12 % between the VPEC and the SEC methods.

The ice hummocks seen in the easterly directions have smaller H and λ, which results in a smaller z0m than in the southerly directions. This is due to the anisotropic nature of the ice hummocks and is confirmed by the eddy-covariance observations, regardless of the season. The extent of the UAV survey allows the application of the drag model for wind directions that rarely occur during the measurement period. This is particularly useful for the development of z0m parametrizations in atmospheric models. Interestingly, the topography at site S5 translates into a wavy pattern of z0m as a function of wind direction, with two local minima at fetch directions of 90 and 180 (Fig. 8).

Figure 9Estimated aerodynamic roughness length (z0m), without accounting for the residual photon backscatter, from ICESat-2 between 16 October 2018 and 6 September 2020 in the lower part of the K-transect, West Greenland. The locations of the automatic weather stations are given by the pink diamonds. The black boxes A and B delineate the areas mapped by UAV photogrammetry. A hillshade of ArcticDEM (Porter et al.2018) is shown as background over the ice sheet.

To summarize, three independent but co-located methods, namely UAV photogrammetry, ICESat-2 laser altimetry, and in situ eddy-covariance measurements, allow us to estimate the aerodynamic roughness of a rough ice surface at a specific site. The comparison confirms our two initial hypotheses: (1) the variability of estimated z0m using in situ observations as a function of wind direction found by T20 is indeed a consequence of the anisotropic topography, and (2) the ICESat-2 data are very well suited to estimate z0m of a rough ice surface in both space and time. Without correcting for the residual scatter in photon elevations, the 1 m resolution ICESat-2 profiles most likely provide a lower bound of roughness, as they underestimate z0m by almost a factor of 2 at the two rough ice locations in areas A and B. On the other hand, an attempt to account for this residual scatter may lead to an overestimated z0m, by a factor that depends on the noise in the raw altimeter data. Nevertheless, given the fact that z0m varies over several orders of magnitude, we deem this method useful to understand the spatio-temporal variability of the aerodynamic roughness length over the GrIS.

4.4 Results: mapping the roughness length z0m using ICESat-2

In this section we apply the elevation profile filtering described in Sect. 3 and the R92 model with parameterized Cd (see Appendix A) to ICESat-2 ATL03 data to model and map the aerodynamic roughness (z0m) over the K-transect. We process nearly 2 years of ICESat-2 ATL03 measurements, taken between 16 October 2018 and 6 September 2020. The results, without accounting for the unresolved photon scatter, are presented in Fig. 9. Within a distance of 10 km from the ice margin, z0m ranges between 10−3 and 10−1 m. There is a clear transition of z0m values that separates the rough (S4, S5, KAN_L, SHR) and smooth (S6 and higher) surface in the ablation zone. Within a distance of several kilometres, z0m can vary by more than 1 order of magnitude, e.g. north of locations S5 and KAN_L (Fig. 9). In order to quantify the variability in time, we group the z0m values from ICESat-2 in two groups, July–September and October–June, which correspond to ICESat-2 cycles 1–3 & 5–8 and 4 & 8 respectively. The average z0m value for the two groups in each 200 m elevation bin is presented in Fig. 10. During summer, the average z0m value is 1×10-2 m below 600 m a.s.l., while it is around 6×10-3 m during the other months. The average roughness approaches its minimum value of 10−4 m above 1000 m a.s.l., regardless of the time period. When the ICESat-2 altimeter does not detect any obstacle, the bulk drag model only accounts for skin friction, which is prescribed as a constant in the model. Interestingly, z0m decreases very near the ice margin, which might be explained by the decreasing ice velocity at the margin, as most of the glaciers in this area are land-terminating.

The measurements described in Table 1 are also included in Fig. 10. The comparison indicates that the satellite product captures the overall variability along the K-transect (Fig. 10). In particular at lower elevations, the modelled z0m is within the range of in situ observations. The in situ roughness z0m can vary due to changing wind direction, but also due to instrumental uncertainty. Especially the smooth sites where the profile method has been used can exhibit large variability, such as at site S9 (see Table 1).

Figure 10Estimated aerodynamic roughness length (z0m) from ICESat-2 between 16 October 2018 until 6 September 2020 along the K-transect. The data were averaged over 200 m elevation bins and two time periods: summer (July–September) and winter (October–June). The variability range denotes the two-sided standard deviation within one elevation class. The in situ observations are described in Table 1.


Unfortunately, the ICESat-2 altimeter is not able to detect obstacles that contribute to form drag at sites S6 (1010 m a.s.l.) and higher. At S6, the surface is flat during winter but becomes rough during summer with ice hummocks with 0.6 m average height (SB08). Unfortunately the horizontal extent of these obstacles is smaller than the ICESat-2 footprint diameter ( 15 m). Higher up, the ice hummocks become even smaller and the surface eventually becomes snow-covered year-round. Nevertheless, snow sastrugi, known to reach up to 0.5 m height at site S10 from photographic evidence, still contribute to form drag. This results in a maximum observed value of z0m=7×10-4 m at sites S9 and S10 (Fig. 10). Using a rough estimate for both H and λ at S6 and S10, based on photographs taken during the end of the ablation season, yields more realistic values for z0m (Fig. A1) than using H and λ from the ICESat-2 elevation profiles. Therefore we conclude that the roughness obstacles are not properly resolved at these locations in the ATL03 data using the algorithm presented in this study, even when the correction using the residual photon scatter is applied. This is mainly due to the limited footprint of the ICESat-2 measurements, but also due to the orientation of the surface features, which limits the detectability of highly anisotropic features from 1D profile measurements. These limitations in the ICESat-2 measurements result in a uniform prescribed value of z0m=1×10-4 m for elevations above 1000 m a.s.l. The algorithm described in Sect. 3 could be adapted to extract these features from the ATL03 data. For instance, smaller-scale obstacles could be retrieved in multiple directions at cross-over points, using the information from multiple ICESat-2 tracks. However, this is beyond the scope of this study, which is to map the aerodynamic roughness of rough ice at large scales.

For now, the implications of these findings for the sensible heat flux, and thus surface ablation, remain to be investigated. The areas with high z0m are in the low-lying ablation zone close to the ice edge, where the highest melt rates are observed. Accounting for the variable roughness might shed light on the drivers of these high melt rates.

5 Conclusions

The aerodynamic roughness of a surface (z0m) in part defines the magnitude of the surface turbulent energy fluxes, yet is often poorly known for glaciers and ice sheets. We adapt the bulk drag partitioning model from Raupach (1992) such that it can be applied to 1D elevation profiles. Forcing this model with 1 m resolution elevation profiles taken from the ICESat-2 satellite laser altimeter z0m becomes a quantifiable and mappable quantity. The model assumes that the surface is composed of regularly spaced, identical obstacles, which all have the same drag coefficient. Despite the fact that the drag coefficient for each individual obstacle remains poorly known, the evaluation in this study against different in situ observations, using different techniques and for different locations and time periods, demonstrates the validity of this model. On the other hand, the use of the model of Lettau (1969) is not recommended over a rough ice surface, as it does not separate the form drag and the skin friction and neglects both the effects of the displacement height and of inter-obstacle sheltering.

Mapping surface obstacles at 1 m resolution using the ICESat-2 altimeter data proves possible, as long as the roughness obstacles are large enough (e.g. crevasses, ice hummocks). Obstacles that are small compared to the ICESat-2 footprint diameter of  15 m, such as ice hummocks found above 1000 m elevation in summer, or snow sastrugi expected year-round at even higher locations on the ice sheet from photographic evidence, are not resolved by the ICESat-2 measurements when used in combination with the methods presented in this study. This translates into a lower bound of z0m10-3 m that can realistically be mapped using this method. Furthermore, accounting for the scatter in the unresolved altimeter signal leads to overestimates of the aerodynamic roughness, as this scatter is a consequence of many different processes that must individually be modelled.

The methods presented in this paper can effectively be used to map z0m at ice sheet elevations below 1000 m. This lower ablation area is also where the contribution of turbulent heat fluxes to surface ablation, and thus runoff, is the largest. As a consequence of the orientation of the ICESat-2 orbit, the modelled z0m must be interpreted as the roughness that would be felt by air flowing in the direction parallel to each laser track. Surfaces of glaciers are often anisotropic, and z0m can vary by over 1 order of magnitude depending on the local wind direction.

The implications of the highly variable aerodynamic roughness for turbulent heat fluxes, and thus surface ablation, remain to be investigated. As current regional climate models typically use constant values for z0m, these implications can be significant, especially in the lower ablation area where most of the surface runoff is generated. This analysis revealed for instance that highly crevassed areas have aerodynamic roughness values over 10−1 m, 2 orders of magnitude larger than typically used in regional climate models over bare ice.

Appendix A: Bulk drag model

For a single roughness element of height H and frontal area Af, placed on a horizontal area Al, Raupach (1992, R92) models the form drag as

(A1) lim λ 0 τ r = F d A l = ρ C d A f A l u ( H ) 2 = ρ C d λ u ( H ) 2 ,

where Fd is the pressure drag force exerted on the obstacle, H the obstacle height, λ the frontal area index, and Cd the drag coefficient of the obstacle. An important uncertainty resides in choosing an accurate value for Cd, due to its dependence on the shape of the obstacles, on the Reynolds number, and on the surface texture. Based on the analysis by Garbrecht et al. (2002) for sea-ice pressure ridges, we choose the following parameterization,

(A2) C d = 1 2 ( 0.185 + 0.147 H ) if H 2.5 m 1 2 0.22 log H 0.2 if H > 2.5 m .

Note that the factor 1/2 is a consequence of a different definition for Cd in Garbrecht et al. (2002) than Eq. (A1).

Similarly, R92 models the skin friction for an unobstructed flat surface as

(A3) lim λ 0 τ s = ρ C s ( z ) u ( z ) 2 ,

where Cs(z) is the drag coefficient of the flat surface, referenced at a height z. Following Andreas (1995), Cs(z) is estimated from the 10 m drag coefficient Cs(10) measured over a flat surface, according to

(A4) C s ( z ) = C s ( 10 ) - 0.5 - 1 κ ln 10 - d z - d - Ψ m ^ ( z ) - 2 ,

with Cs(10)=1.2071×10-3, which yields z0m=10-4 m for a perfectly flat surface in this model. This value was chosen as it is the minimum value estimated using in situ observations by Smeets and Van den Broeke (2008) during winter over different snow surfaces.

In reality, the total surface shear stress is the sum of both the form drag on each individual obstacle (τr) and the skin friction on the underlying surface (τs) (Eq. 3). Furthermore, an additional complexity arises at increasing obstacle frontal area index (λ), as the obstacles may effectively shelter a part of the surface and each other, thereby reducing both the skin friction and the form drag. Based on the previous work of Arya (1975), and on scaling arguments of the effective shelter volume, R92 includes sheltering and models the total surface shear stress over multiple obstacles as


where c=0.25 is an empirical constant that determines the sheltering efficiency. The latter equation may be written in the form

(A7) X e - X = a ,



which is solved iteratively using Xi+1=eXi and X0=a, after R92. The solution yields u(H)u*.

The conversion of u(H)u* to z0m,R92 is finally possible using the semi-logarithmic wind profile Eq. (2) and referencing the wind speed at z=H. However, an expression for the displacement height d and the roughness sublayer profile function Ψm(z)^ is still required. For the displacement height, the simplified expression by Raupach (1994) is used:

(A10) d = H 1 - 1 - exp ( - c d 1 λ ) c d 1 λ ,

with cd1=7.5, which is then used to derive the value for Ψr at height z=H using the following expression:

(A11) Ψ m ( H ) ^ = log c w - 1 + c w ,


(A12) c w = z * - d H - d ,

in which z* is the upper height of the roughness layer. Raupach (1994) empirically determined that cw=2, which yields Ψm(H)^=0.193.

To summarize, the aerodynamic roughness length z0m of an elevation profile of length L is modelled according to the following steps.

  • The elevation profile is high-pass filtered using a cut-off wavelength Λ=35 m.

  • The obstacle height (H) is set to twice the standard deviation of the filtered profile.

  • Each group of consecutive positive heights is defined as a single obstacle, which yields the number of obstacles (f) per profile length (L).

  • The frontal area index (λ) is calculated using Eq. (6).

  • The displacement height d is estimated with Eq. (A10).

  • Cs(z) is evaluated at z=H using Eq. (A4), while Cd is parameterized using Eq. (A2).

  • u(H)u* is estimated from Eqs. (A7)–(A9).

  • z0m/H is estimated by evaluating the logarithmic wind profile at a height z=H, using Eq. (2)

Following the steps above, z0m can be estimated for any H and λ, which is done in Fig. A1. At areas A and B and site S5, H and λ are estimated from the UAV surveys and from ICESat-2 data. At site S6, we assume that H=0.6±0.1 m and λ=0.045±0.015, based on photographs taken during the end of the ablation season. At the highest site S10, we assume that H=0.3±0.2 m and λ=0.02±0.01, which are typical values for sastrugi (Andreas1995).

Figure A1Estimated z0m using the R92 model with parameterized Cd (Appendix A), as a function of obstacle height H and frontal area index λ. The solid squares denote the estimated H and λ at three sites using UAV surveys. The dashed squares are estimates based on photographs taken at the end of the ablation season. See Fig. 1 for the location of each site.


Other attempts have been made to relate z0m to the geometry of multiple surface roughness elements. For instance Lettau (1969, L69) empirically relate z0m to the average frontal area index of the roughness obstacles, which has been adapted by Munro (1989) for the surface of a glacier:

(A13) z 0m,L69 = 2 C d H A f A l = 2 C d H λ .

Note that Eq. (A2) was adapted in order to be consistent with the definition of Cd in Eq. (A1). In fact, Macdonald et al. (1998, M98) have shown that Eq. (A13) can be obtained by assuming that there is only form drag and by setting d=0, Ψm(z)^=0 and Cd=0.25. By including the displacement height d, M98 is able to reproduce the non-linear feature of the z0mH=f(λ) curve:

(A14) z 0m,M98 = ( H - d ) exp - C d κ 2 λ 1 - d H - 0.5 .
Appendix B: Sensitivity experiments

B1 Cutoff wavelength Λ

We find that the optimal value of the cutoff wavelength for the high-pass filter is Λ=35 m. This may be explained by the fact that the resulting filtered topography using Λ=35 m still contains most ( 80 %) of the total variance of the slope spectrum. The latter is defined as the power spectral density of the first derivative of the elevation profile. A sensitivity experiment using different values for Λ at S5 can be found in Fig. B1. Changing the value for Λ strongly impacts the estimated H (Fig. B1c), as the elevation profiles considered here contain information at all wavelengths (Fig. B1a). On the other hand, increasing the value for Λ above 35 m does not significantly affect the estimated frontal area index λ (Fig. B1b). Overall, increasing Λ from 10 to 50 m increases the modelled z0m from 7.6×10-4 to 2.8×10-2 m at S5, in the direction 184 that matches the ICESat-2 track (Fig. B1d).

B2 ATL03 kriging parameters

In order to interpolate the geolocated photon product ATL03 in a regular 1 m resolution elevation profile, a fixed set of interpolation parameters was used, referred to as the default set. These are the median filter coefficients in Eq. (7) qlow=1 and qhigh=2, the median filter window length of 50 m, the choice of a Gaussian covariance function with a radius of 15 m in the kriging equations, and the maximum distance of photon distance to each regular grid point of 15 m.

This default parameter set was found to give robust results, even when only medium- or low-confidence photons are present in the ATL03 data. A sensitivity experiment by varying each parameter separately in a 200 m portion of areas A and B is given in Fig. B2. While the interpolated ATL03 elevation still misses small-scale features present in the UAV data, varying each parameter does not give improved results (Fig. B2).

Figure B1(a) Filtered elevation profile in direction 186, (b) estimated obstacle frontal area index, (c) estimated obstacle height, and (d) modelled aerodynamic roughness length at site S5 for different high-pass cutoff wavelengths Λ. See Fig. 8 in main text for the labels in panel (d).


Figure B2Elevation profiles in a 200 m portion of area A (a1, b1) and area B (a2, b2). The top panels contain the ATL03 data sorted in confidence levels (dots), the ATL06 data (pink triangles), the profiles measured by UAV photogrammetry (orange line), and the 1 m interpolated ATL03 data using the default settings used in the main text (blue line). The bottom panels contain the 1 m interpolated AT03 data using different origins and photon filtering settings.


Data availability

The data used in the figures can be downloaded here: (Van Tiggelen2020). Additional data can be obtained from the authors without conditions.

Author contributions

PS, MRvB, CHR, and WI secured funding for this study and for the field experiments. JFS piloted the UAV in the field. PS and MvT acquired and processed the eddy-covariance data. EJN, JFS, and WI made the UAV elevation models. MvT and BW processed the ICESat-2 data. MvT, with the help of all the co-authors, designed the study and wrote the first draft. All authors participated in the review and in the editing of the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


The authors thank all the people and institutes that help maintain the instruments in the field. We are grateful to Nanna Karlsson, Dirk van As, and Giorgio Cover for their support in the field. This work is funded by the Utrecht University and by the Netherlands Polar Programme (NPP), of the Netherlands Organisation of Scientific Research, section Earth and Life Sciences (NWO/ALWOP.431). This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. Jakob F. Steiner and Walter W. Immerzeel acknowledge support by the Netherlands Organization for Scientific Research NWO (016.181.308) and European Research Council (676819). The views and interpretations in this publication are those of the authors and are not necessarily attributable to ICIMOD.

Financial support

This research has been supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant no. ALWOP.431).

Review statement

This paper was edited by Benjamin Smith and reviewed by Christof Lüpkes, Shane Grigsby, and one anonymous referee.


Andreas, E. L.: Air-ice drag coefficients in the western Weddell Sea 2. A model based on form drag and drifting snow, J. Geophys. Res., 100, 4833–4843,, 1995. a, b, c

Arya, S. P. S.: A Drag Partition Theory for Determining the Large-Scale Roughness Parameter and Wind Stress on the Arctic Pack Ice, J. Geophys. Res., 80, 3447–3454, 1975. 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

Brunt, K. M., Neumann, T. A., and Smith, B. E.: Assessment of ICESat-2 Ice Sheet Surface Heights, Based on Comparisons Over the Interior of the Antarctic Ice Sheet, Geophys. Res. Lett., 46, 13072–13078,, 2019. 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 Surf., 125, e2019JF005167,, 2019. a, b

Cooper, M. G., Smith, L. C., Rennermalm, A. K., Tedesco, M., Muthyala, R., Leidman, S. Z., Moustafa, S. E., and Fayne, J. V.: Spectral attenuation coefficients from measurements of light transmission in bare ice on the Greenland Ice Sheet, The Cryosphere, 15, 1931–1953,, 2021. a

Fausto, R. S., Van As, D., Box, J. E., Colgan, W., and Langen, P. L.: Quantifying the surface energy fluxes in South Greenland during the 2012 high melt episodes using in-situ observations, Front. Earth Sci., 4, 1–9,, 2016. a

Fausto, R. S., van As, D., Mankoff, K. D., Vandecrux, B., Citterio, M., Ahlstrøm, A. P., Andersen, S. B., Colgan, W., Karlsson, N. B., Kjeldsen, K. K., Korsgaard, N. J., Larsen, S. H., Nielsen, S., Pedersen, A. Ø., Shields, C. L., Solgaard, A. M., and Box, J. E.: PROMICE automatic weather station data, Earth Syst. Sci. Data Discuss. [preprint],, in review, 2021. 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

Forrer, J. and Rotach, M. W.: On the turbulence structure in the stable boundary layer over the Greenland ice sheet, Bound.-Lay. Meteorol., 85, 111–136,, 1997. a

Garbrecht, T., Lüpkes, C., Augstein, E., and Wamser, C.: Influence of a sea ice ridge on low-level airflow, J. Geophys. Res.-Atmos., 104, 24499–24507,, 1999. a

Garbrecht, T., Lüpkes, C., Hartmann, J., and Wolff, M.: Atmospheric drag coefficients over sea ice – Validation of a parameterisation concept, Tellus A, 54, 205–219,, 2002. a, b

Gardner, C. S.: Target signatures for laser altimeters: an analysis, Appl. Opt., 21, 448,, 1982. a

Garratt, J. R.: The atmospheric boundary layer, Cambridge University Press, Cambridge, 1992. a, b

Hanssen-Bauer, I. and Gjessing, Y. T.: Observations and model calculations of aerodynamic drag on sea ice in the Fram Strait, Tellus A, 40, 151–161,, 1988. a

Harman, I. N. and Finnigan, J. J.: A simple unified theory for flow in the canopy and roughness sublayer, Bound.-Lay. Meteorol., 123, 339–363,, 2007. a, b, c

Heinemann, G.: The KABEG'97 field experiment: An aircraft-based study of Katabatic wind dynamics over the Greenland ice sheet, Bound.-Lay. Meteorol., 93, 75–116,, 1999. a

Hengl, T.: A practical guide to geostatistical mapping, Office for Official Publications of the European Communities, Luxembourg, 2009. a

Herzfeld, U. C., Box, J. E., Steffen, K., Mayer, H., Caine, N., and Losleben, M. V.: A Case Study or the Influence of Snow and Ice Surface Roughness on Melt Energy, Z. Gletscherkd. Glazialgeol, 39, 1–42, 2006. a

Herzfeld, U. C., Trantow, T., Lawson, M., Hans, J., and Medley, G.: Surface heights and crevasse morphologies of surging and fast-moving glaciers from ICESat-2 laser altimeter data - Application of the density-dimension algorithm (DDA- ice) and evaluation using airborne altimeter and Planet SkySat data, Sci. Remote Sens., 3, 104743,, 2020. a

Immerzeel, W. W., Kraaijenbrink, P. D. A., Shea, J. M., Shrestha, A. B., Pellicciotti, F., Bierkens, M. F. P., and Jong, S. M. D.: Remote Sensing of Environment High-resolution monitoring of Himalayan glacier dynamics using unmanned aerial vehicles, Remote Sens. Environ., 150, 93–103,, 2014. a

Jackson, P.: On the displacement height in the logarithmic velocity profile, J. Fluid Mech., 111, 15–25, 1981. a

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

Kean, J. W. and Smith, J. D.: Form drag in rivers due to small-scale natural topographic features: 2. Irregular sequences, J. Geophys. Res.-Earth Surf., 111, 1–15,, 2006. a

King, M. D., Howat, I. M., Candela, S. G., Noh, M. J., Jeong, S., Noël, B. P. Y., Broeke, M. R. V. D., Wouters, B., and Negrete, A.: Dynamic ice loss from the Greenland Ice Sheet driven by sustained glacier retreat, Commun. Earth Environ., 1, 1–7,, 2020. a

Kljun, N., Calanca, P., Rotach, M. W., and Schmid, H. P.: A simple two-dimensional parameterisation for Flux Footprint Prediction (FFP), Geosci. Model Dev., 8, 3695–3713,, 2015. a, b

Kraaijenbrink, P. D., Shea, J. M., Pellicciotti, F., Jong, S. M., and Immerzeel, W. W.: Object-based analysis of unmanned aerial vehicle imagery to map and characterise surface features on a debris-covered glacier, Remote Sens. Environ., 186, 581–595,, 2016. a

Kuipers Munneke, P., Smeets, C. J. P. P., Reijmer, C. H., Oerlemans, J., van de Wal, R. S. W., and van den Broeke, M. R.: The K-transect on the western Greenland Ice Sheet: Surface energy balance (2003–2016), Arctic, Antarct. Alp. Res., 50, S100003,, 2018. a

Kurtz, N. T., Markus, T., Cavalieri, D. J., Krabill, W., Sonntag, J. G., and Miller, J.: Comparison of ICESat data with airborne laser altimeter measurements over arctic sea ice, IEEE T. Geosci. Remote, 46, 1913–1924,, 2008. a, b

Lenaerts, J. T. M., Smeets, C. J. P. P., Nishimura, K., Eijkelboom, M., Boot, W., van den Broeke, M. R., and van de Berg, W. J.: Drifting snow measurements on the Greenland Ice Sheet and their application for model evaluation, The Cryosphere, 8, 801–814,, 2014. a, b

Lettau, H.: Note on Aerodynamic Roughness-Parameter Estimation on the Basis of Roughness-Element Description, J. Appl. Meteorol. Climatol., 8, 828–832, 1969. a, b, c, d, e

Li, Q., Bou‐Zeid, E., Grimmond, S., Zilitinkevich, S., and Katul, G.: Revisiting the Relation Between Momentum and Scalar Roughness Lengths of Urban Surfaces, Q. J. Roy. Meteor. Soc., 146, 3144–3164,, 2020. a

Lüpkes, C. and Gryanik, V. M.: A stability-dependent parametrization of transfer coefficients formomentum and heat over polar sea ice to be used in climate models, J. Geophys. Res., 120, 552–581,, 2015. a, b

Lüpkes, C., Gryanik, V. M., Hartmann, J., and Andreas, E. L.: A parametrization, based on sea ice morphology, of the neutral atmospheric drag coefficients for weather prediction and climate models, J. Geophys. Res.-Atmos., 117, D13112,, 2012. a, b, c

Lüpkes, C., Gryanik, V. M., Rösel, A., Birnbaum, G., and Kaleschke, L.: Effect of sea ice morphology during Arctic summer on atmospheric drag coefficients used in climate models, Geophys. Res. Lett., 40, 446–451,, 2013. a

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, b, c

Markus, T., Neumann, T., Martino, A., Abdalati, W., Brunt, K., Csatho, B., Farrell, S., Fricker, H., Gardner, A., Harding, D., Jasinski, M., Kwok, R., Magruder, L., Lubin, D., Luthcke, S., Morison, J., Nelson, R., Neuenschwander, A., Palm, S., Popescu, S., Shum, C. K., Schutz, B. E., Smith, B., Yang, Y., and Zwally, J.: The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2): Science requirements, concept, and implementation, Remote Sens. Environ., 190, 260–273,, 2017. a

Meesters, A. G., Bink, N. J., Vugts, H. F., Cannemeijer, F., and Henneken, E. A.: Turbulence observations above a smooth melting surface on the Greenland ice sheet, Bound.-Lay. Meteorol., 85, 81–110,, 1997. a, b

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

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244,, 2019. a

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

Neumann, T. A., Martino, A. J., Markus, T., Bae, S., Bock, M. R., Brenner, A. C., Brunt, K. M., Cavanaugh, J., Fernandes, S. T., Hancock, D. W., Harbeck, K., Lee, J., Kurtz, N. T., Luers, P. J., Luthcke, S. B., Magruder, L., Pennington, T. A., Ramos-Izquierdo, L., Rebold, T., Skoog, J., and Thomas, T. C.: The Ice, Cloud, and Land Elevation Satellite – 2 mission: A global geolocated photon product derived from the Advanced Topographic Laser Altimeter System, Remote Sens. Environ., 233, 111325,, 2019. a, b

Noël, B., van de Berg, W. J., Lhermitte, S., and van den Broeke, M. R.: Rapid ablation zone expansion amplifies north Greenland mass loss, Sci. Adv., 5, eaaw0123,, 2019. a

Nolin, A. W. and Mar, E.: Arctic sea ice surface roughness estimated from multi-angular reflectance satellite imagery, Remote Sens., 11, 1–12,, 2019. a

Petty, A. A., Tsamados, M. C., and Kurtz, N. T.: Atmospheric form drag coefficients over Arctic sea ice using remotely sensed ice topography data, spring 2009–2015, J. Geophys. Res.-Earth Surf., 122, 1472–1490,, 2017. a

Porter, C., Morin, P., Howat, I., Noh, M.-J., Bates, B., Peterman, K., Keesey, S., Schlenk, M., Gardiner, J., Tomko, K., Willis, M., Kelleher, C., Cloutier, M., Husby, E., Foga, S., Nakamura, H., Platson, M., Wethington, Michael, J., Williamson, C., Bauer, G., Enos, J., Arnold, G., Kramer, W., Becker, P., Doshi, A., D'Souza, C., Cummens, P., Laurier, F., and Bojesen, M.: ArcticDEM, Harvard Dataverse,, 2018. a, b

Raupach, M.: Drag and drag partition on rough surfaces, Bound.-Lay. Meteorol., 60, 375–395, 1992. a, b, c, d, e, f, g

Raupach, M. R.: Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index, Bound.-Lay. Meteorol., 71, 211–216,, 1994. a, b

Shao, Y. and Yang, Y.: A theory for drag partition over rough surfaces, J. Geophys. Res.-Earth Surf., 113, 1–9,, 2008. a, b

Shepherd, A., Ivins, E., Rignot, E., Smith, B., van den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N., A, G., Agosta, C., Ahlstrøm, A., Babonis, G., Barletta, V. R., Bjørk, A. A., Blazquez, A., Bonin, J., Colgan, W., Csatho, B., Cullather, R., Engdahl, M. E., Felikson, D., Fettweis, X., Forsberg, R., Hogg, A. E., Gallee, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P. L., Lecavalier, B., Loomis, B., Luthcke, S., McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mottram, R., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noël, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., Sandberg Sørensen, L., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K. W., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., van de Berg, W. J., van der Wal, W., van Wessem, M., Vishwakarma, B. D., Wiese, D., Wilton, D., Wagner, T., Wouters, B., and Wuite, J.: Mass balance of the Greenland Ice Sheet from 1992 to 2018, Nature, 579, 233–239,, 2020. a

Smeets, C. and Van den Broeke, M. R.: Temporal and spatial variations of the aerodynamic roughness length in the ablation zone of the greenland ice sheet, Bound.-Lay. Meteorol., 128, 315–338,, 2008.  a, b, c, d

Smeets, C. J. P. P., Duynkerke, P. G., and Vugts, H. F.: Observed Wind Profiles and Turbulence Fluxes over an ice Surface with Changing Surface Roughness, Boundary-Lay. Meteorol., 92, 99–121,, 1999. a, b, c

Smeets, P. C. J. P., Kuipers Munneke, P., van As, D., van den Broeke, M. R., Boot, W., Oerlemans, H., Snellen, H., Reijmer, C. H., and van de Wal, R. S. W.: The K-transect in west Greenland: automatic weather station data (1993–2016), Arctic, Antarct. Alp. Res., 50, S100002,, 2018. a, b, c

Smith, B., Fricker, H. A., Holschuh, N., Gardner, A. S., Adusumilli, S., Brunt, K. M., Csatho, B., Harbeck, K., Huth, A., Neumann, T., Nilsson, J., and Siegfried, M. R.: Land ice height-retrieval algorithm for NASA's ICESat-2 photon-counting laser altimeter, Remote Sens. Environ., 233, 111352,, 2019. a

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 Surf., 121, 748–766,, 2016. a, b

van de Wal, R. S. W., Boot, W., Smeets, C. J. P. P., Snellen, H., van den Broeke, M. R., and Oerlemans, J.: Twenty-one years of mass balance observations along the K-transect, West Greenland, Earth Syst. Sci. Data, 4, 31–35,, 2012. a

Van den Broeke, M. R.: Characteristics of the lower ablation zone of the West Greenland ice sheet for energy-balance modelling, Ann. Glaciol., 23, 7–13,, 1996. a, b

Van Tiggelen, M.: mvantiggelen/TC_2020: Release 001 (Version 001), Zenodo [Data set],, 2020. a

Van Tiggelen, M., Smeets, P. C., Reijmer, C. H., and van den Broeke, M. R.: A Vertical Propeller Eddy-Covariance Method and Its Application to Long-term Monitoring of Surface Turbulent Fluxes on the Greenland Ice Sheet, Bound.-Lay. Meteorol., 176, 441–463,, 2020. a, b, c, d, e, f

Yang, X. I., Sadique, J., Mittal, R., and Meneveau, C.: Exponential roughness layer and analytical model for turbulent boundary layer flow over rectangular-prism roughness elements, J. Fluid Mech., 789, 127–165,, 2016. a

Yi, D., Zwally, H. J., and Sun, X.: ICESat measurement of Greenland ice sheet surface slope and roughness, Ann. Glaciol., 42, 83–89,, 2005. a, b

Zhu, X. and Anderson, W.: Turbulent flow over urban-like fractals: Prognostic roughness model for unresolved generations, J. Turbul., 19, 995–1016,, 2019. a

Short summary
We developed a method to estimate the aerodynamic properties of the Greenland Ice Sheet surface using either UAV or ICESat-2 elevation data. We show that this new method is able to reproduce the important spatiotemporal variability in surface aerodynamic roughness, measured by the field observations. The new maps of surface roughness can be used in atmospheric models to improve simulations of surface turbulent heat fluxes and therefore surface energy and mass balance over rough ice worldwide.