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

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 (z_0m) 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 z_0m, ranging between 10⁻⁴ m for a smooth snow surface and 10⁻¹ 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 Rotach, 1997; Heinemann, 1999). Additionally, the surface of the GrIS close to the ice edge is very rough (Yi et al., 2005; Smeets and Van den Broeke, 2008). 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. Garratt, 1992), thereby enhancing surface melt (Van den Broeke, 1996; 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 (z_0m) 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 z_0m

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 (Garratt, 1992; Raupach, 1992; Harman and Finnigan, 2007). As such, the vertical profile of horizontal wind speed (u(z)) over a rough surface can be written as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}u\left(z\right)={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}& \left[\ln \left({\displaystyle \frac{z-d}{z_{\text{0m}}}}\right)-{\mathrm{\Psi }}_{\text{m}}\left({\displaystyle \frac{z-d}{L_{\text{o}}}}\right)\right.\\ & \left.+{\mathrm{\Psi }}_{\text{m}}\left({\displaystyle \frac{z_{\text{0m}}}{L_{\text{o}}}}\right)+\widehat{{\mathrm{\Psi }}_{\text{m}}}\left(z\right)\right],\end{array}\end{array}$$

where z is the height above the surface, $u_{*}={\left(\frac{\mathit{\tau }}{\mathit{\rho }}\right)}^{\mathrm{0.5}}$ is the friction velocity, ρ the air density, κ=0.4 the von Kármán constant, τ the total surface shear stress, and z_0m 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 (Jackson, 1981). z_0m 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 ${\mathrm{\Psi }}_{\text{m}}\left(\frac{z-d}{L_{\text{o}}}\right)$ and $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}$, respectively, where L_o is the Obukhov length. The hat notation is used for the roughness layer quantities, as in Harman and Finnigan (2007). Above the roughness sublayer, $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}=\mathrm{0}$.

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

$$\begin{array}{}\text{(2)}& z_{\text{0m}}=(z-d){\left[\exp \left(\mathit{\kappa }{\displaystyle \frac{u\left(z\right)}{u_{*}}}-\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)})\right)\right]}^{-\mathrm{1}}.\end{array}$$

Hence, the process of finding z_0m is equivalent to finding d, $\frac{u\left(z\right)}{u_{*}}$, and $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}$ simultaneously.

2.2 Bulk drag model of z_0m

The main task is to model the total surface shear stress $\mathit{\tau }=\mathit{\rho }{u}_{*}^{\mathrm{2}}$, which for a rough surface is the sum of both form drag τ_r and skin friction τ_s:

$$\begin{array}{}\text{(3)}& \mathit{\tau }={\mathit{\tau }}_{\text{r}}+{\mathit{\tau }}_{\text{s}}.\end{array}$$

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

$$\begin{array}{}\text{(4)}& \mathit{\lambda }={\displaystyle \frac{A_{\text{f}}}{A_{\text{l}}}},\end{array}$$

with A_f the frontal area of the roughness obstacles perpendicular to the flow and A_l 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 (Munro, 1989; Smith et al., 2016) (Fig. 3). This avoids the use of empirical formulas for the estimation of z_0m 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 3 Steps 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

$$\begin{array}{}\text{(5)}& H=\mathrm{2}\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_z},\end{array}$$

in which $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_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)

$$\begin{array}{}\text{(6)}& \mathit{\lambda }=f{\displaystyle \frac{Hw}{Lw}}={\displaystyle \frac{fH}{L}},\end{array}$$

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 z_0m 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 z_0m 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 1 Description of z_0m 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.

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 (z_0m) is calculated using Eq. (2).

We only select data taken during near-neutral conditions ($z/L_{\text{o}}<\mathrm{0.1}$), and we assume that the measurements are taken above the roughness layer, i.e. $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}=\mathrm{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 Finnigan, 2007). On the other hand, when applying the drag model to estimate z_0m (Appendix A), the correction factor $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}$ 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 $[\mathrm{80}{}^{\circ };\mathrm{200}{}^{\circ }]$ interval. In the following sections we average ln(z_0m), which is of interest for the determination of the vertical profile of horizontal wind speed (Eq. 1). Additional in situ averaged z_0m 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 2 Description of DEMs obtained by UAV photogrammetry and of the corresponding overlapping ICESat-2 laser beams.

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,

$$\begin{array}{}\text{(7)}& \langle z\rangle -{\displaystyle \frac{q_{\text{low}}}{\mathrm{0.6745}}}\langle |z-\langle z\rangle |\rangle \le z\le \langle z\rangle +{\displaystyle \frac{q_{\text{high}}}{\mathrm{0.6745}}}\langle |z-\langle z\rangle |\rangle ,\end{array}$$

where 〈z〉 denotes the median of z within a moving window. We choose q_low=1 and q_high=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 (Gardner, 1982).

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

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{\text{sub}}={\left({\mathit{\sigma }}_{\text{ph,res}}^{\mathrm{2}}-{\mathit{\sigma }}_{\text{i}}^{\mathrm{2}}\right)}^{\mathrm{0.5}}/\mathrm{2},\end{array}$$

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:

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_{\text{tot}}=\sqrt{{\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\text{res}}}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{sub}}^{\mathrm{2}}},\end{array}$$

in which $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\text{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

$$\begin{array}{}\text{(10)}& H_{\text{corr}}=\mathrm{2}{\mathit{\sigma }}_{\text{tot}}.\end{array}$$

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 z_0m 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 z_0m 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 z_0m 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 5 Modelled z_0m 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: C_d=0.25 (solid lines) and C_d=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 C_d 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 z_0m are calculated for different fetch directions as depicted in Fig. 2. Figure 5 compares both the estimated z_0m from in situ observations and the modelled z_0m at the end of the ablation season, as a function of the measured obstacle frontal area index λ. The L69 model (Eq. A13) overestimates z_0m for λ<0.04 at this location (Fig. 5, blue line). In accordance with L69, the drag coefficient of an individual obstacle C_d=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 z_0m for λ>0.04, which can be explained by the neglect of the displacement that is compensated for by too small C_d for these fetch directions. The method by M98 (Eq. A14) does account for the displacement height, and, while using the same drag coefficient C_d=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 C_d=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 C_d 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 C_d 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 C_d from Garbrecht et al. (1999) (Eq. A2), for which C_d 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 C_d.

Figure 6 Estimated 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 C_d allows for some variability in modelled z_0m for the same λ but is still not able to reproduce the eddy-covariance observations (Fig. 5). We attribute this to the parametrization of C_d, 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 Anderson, 2019; 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 C_d for mapping z_0m 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 (z_0m). (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 z_0m 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 z_0m 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 z_0m by approximately $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ m (Fig. 7). On average, the uncorrected ICESat-2 values underestimate z_0m by $\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{3}}$ m for area A and $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{3}}$ m for area B, which corresponds to ≈ 40 % and ≈ 36 % of the average z_0m estimated by the UAV at these two sites.

Figure 8 (a) Drag model evaluation at site S5. (b) Drag coefficient for form drag (C_d) used in the bulk drag model (black line) or required to perfectly fit the observations. The solid orange line is the modelled z_0m 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 Robson, 2014). Both H and λ are smaller in the satellite profile than in the UAV profile, but the modelled z_0m agrees qualitatively with that estimated using observations from the AWS S5 during March–April. During this time period, z_0m 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 z_0m 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 z_0m using in situ observations. Across all the wind directions available in the measurements, z_0m using the UAV profiles is underestimated by $\mathrm{7.6}\times {\mathrm{10}}^{-\mathrm{3}}$ m, or 28 % of the average value estimated by the SEC method in September 2019 ($\mathrm{2.65}\times {\mathrm{10}}^{-\mathrm{2}}$ m). The comparison improves when comparing the modelled z_0m to VPEC measurements from September 2016 and August 2017 (T20), the model now overestimates the estimated z_0m by in situ observations by $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{3}}$ m or 9 % of the observed value ($\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{2}}$ m). As these data contain more wind directions, the overestimation of z_0m 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 z_0m in time but also the uncertainty in the field measurements. The difference in averaged estimated z_0m 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 z_0m 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 z_0m parametrizations in atmospheric models. Interestingly, the topography at site S5 translates into a wavy pattern of z_0m as a function of wind direction, with two local minima at fetch directions of 90 and 180^∘ (Fig. 8).

Figure 9 Estimated aerodynamic roughness length (z_0m), 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 z_0m 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 z_0m 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 z_0m 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 z_0m, by a factor that depends on the noise in the raw altimeter data. Nevertheless, given the fact that z_0m 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 z_0m using ICESat-2

In this section we apply the elevation profile filtering described in Sect. 3 and the R92 model with parameterized C_d (see Appendix A) to ICESat-2 ATL03 data to model and map the aerodynamic roughness (z_0m) 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, z_0m ranges between 10⁻³ and 10⁻¹ m. There is a clear transition of z_0m 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, z_0m 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 z_0m 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 z_0m value for the two groups in each 200 m elevation bin is presented in Fig. 10. During summer, the average z_0m value is $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ m below 600 m a.s.l., while it is around $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}$ m during the other months. The average roughness approaches its minimum value of 10⁻⁴ 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, z_0m 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 z_0m is within the range of in situ observations. The in situ roughness z_0m 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 10 Estimated aerodynamic roughness length (z_0m) 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 $z_{\text{0m}}=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{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 z_0m (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 $z_{\text{0m}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{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 z_0m 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 (z_0m) 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 z_0m 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 $z_{\text{0m}}\approx {\mathrm{10}}^{-\mathrm{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 z_0m 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 z_0m 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 z_0m 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 z_0m, 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⁻¹ 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 A_f, placed on a horizontal area A_l, Raupach (1992, R92) models the form drag as

$$\begin{array}{}\text{(A1)}& \underset{\mathit{\lambda }\to \mathrm{0}}{lim}{\mathit{\tau }}_{\text{r}}={\displaystyle \frac{F_{\text{d}}}{A_{\text{l}}}}=\mathit{\rho }{C}_{\text{d}}{\displaystyle \frac{A_{\text{f}}}{A_{\text{l}}}}u(H{)}^{\mathrm{2}}=\mathit{\rho }{C}_{\text{d}}\mathit{\lambda }u(H{)}^{\mathrm{2}},\end{array}$$

where F_d is the pressure drag force exerted on the obstacle, H the obstacle height, λ the frontal area index, and C_d the drag coefficient of the obstacle. An important uncertainty resides in choosing an accurate value for C_d, 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,

$$\begin{array}{}\text{(A2)}& C_{\text{d}}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{0.185}+\mathrm{0.147}H)& \text{if}H\le \mathrm{2.5}\text{m}\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{0.22}\log \left(\frac{H}{\mathrm{0.2}}\right)\right)& \text{if}H>\mathrm{2.5}\text{m}.\end{array}\right.\end{array}$$

Note that the factor $\mathrm{1}/\mathrm{2}$ is a consequence of a different definition for C_d in Garbrecht et al. (2002) than Eq. (A1).

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

$$\begin{array}{}\text{(A3)}& \underset{\mathit{\lambda }\to \mathrm{0}}{lim}{\mathit{\tau }}_{\text{s}}=\mathit{\rho }{C}_{\text{s}}\left(z\right)u(z{)}^{\mathrm{2}},\end{array}$$

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

$$\begin{array}{}\text{(A4)}& C_{\text{s}}\left(z\right)={\left[C_{\text{s}}(\mathrm{10}{)}^{-\mathrm{0.5}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }}}\left(\ln \left({\displaystyle \frac{\mathrm{10}-d}{z-d}}\right)-\widehat{{\mathrm{\Psi }}_{\text{m}}}\left(z\right)\right)\right]}^{-\mathrm{2}},\end{array}$$

with $C_{\text{s}}\left(\mathrm{10}\right)=\mathrm{1.2071}\times {\mathrm{10}}^{-\mathrm{3}}$, which yields $z_{\text{0m}}={\mathrm{10}}^{-\mathrm{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

$$\begin{array}{}\text{(A5)}& {\displaystyle }\mathit{\tau }\left(\mathit{\lambda }\right)={\mathit{\tau }}_{\text{s}}\left(\mathit{\lambda }\right)+{\mathit{\tau }}_{\text{r}}\left(\mathit{\lambda }\right)\text{(A6)}& {\displaystyle }\begin{array}{rl}& =\mathit{\rho }u(H{)}^{\mathrm{2}}\left[C_{\text{s}}\left(H\right)\exp \left(-c\mathit{\lambda }{\displaystyle \frac{u\left(H\right)}{u_{*}}}\right)\right.,\\ & \left.+\mathit{\lambda }{C}_{\text{d}}\exp \left(-c\mathit{\lambda }{\displaystyle \frac{u\left(H\right)}{u_{*}}}\right)\right],\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A7)}& X{e}^{-X}=a,\end{array}$$

with

$$\begin{array}{}\text{(A8)}& {\displaystyle }& {\displaystyle }X={\displaystyle \frac{c\mathit{\lambda }}{\mathrm{2}}}{\displaystyle \frac{u\left(H\right)}{u_{*}}},\text{(A9)}& {\displaystyle }& {\displaystyle }a={\displaystyle \frac{c\mathit{\lambda }}{\mathrm{2}}}{\left(C_{\text{s}}+\mathit{\lambda }{C}_{\text{d}}\right)}^{-\mathrm{0.5}},\end{array}$$

which is solved iteratively using $X_{i+\mathrm{1}}=e^{X_i}$ and X₀=a, after R92. The solution yields $\frac{u\left(H\right)}{u_{*}}$.

The conversion of $\frac{u\left(H\right)}{u_{*}}$ to z_0m,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 $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}$ is still required. For the displacement height, the simplified expression by Raupach (1994) is used:

$$\begin{array}{}\text{(A10)}& d=H\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\exp (-\sqrt{c_{d\mathrm{1}}\mathit{\lambda }})}{\sqrt{c_{d\mathrm{1}}\mathit{\lambda }}}}\right],\end{array}$$

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

$$\begin{array}{}\text{(A11)}& \widehat{{\mathrm{\Psi }}_{\text{m}}\left(H\right)}=\log \left(c_{\text{w}}\right)-\mathrm{1}+c_{\text{w}},\end{array}$$

where

$$\begin{array}{}\text{(A12)}& c_{\text{w}}={\displaystyle \frac{z_{*}-d}{H-d}},\end{array}$$

in which z_* is the upper height of the roughness layer. Raupach (1994) empirically determined that c_w=2, which yields $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(H\right)}=\mathrm{0.193}$.

To summarize, the aerodynamic roughness length z_0m 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).

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

• $\frac{u\left(H\right)}{u_{*}}$ is estimated from Eqs. (A7)–(A9).

• $z_{\text{0m}}/H$ is estimated by evaluating the logarithmic wind profile at a height z=H, using Eq. (2)

Following the steps above, z_0m 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=\mathrm{0.6}\pm \mathrm{0.1}$ m and $\mathit{\lambda }=\mathrm{0.045}\pm \mathrm{0.015}$, based on photographs taken during the end of the ablation season. At the highest site S10, we assume that $H=\mathrm{0.3}\pm \mathrm{0.2}$ m and $\mathit{\lambda }=\mathrm{0.02}\pm \mathrm{0.01}$, which are typical values for sastrugi (Andreas, 1995).

Figure A1 Estimated z_0m using the R92 model with parameterized C_d (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 z_0m to the geometry of multiple surface roughness elements. For instance Lettau (1969, L69) empirically relate z_0m to the average frontal area index of the roughness obstacles, which has been adapted by Munro (1989) for the surface of a glacier:

$$\begin{array}{}\text{(A13)}& z_{\text{0m,L69}}=\mathrm{2}{C}_{\text{d}}H{\displaystyle \frac{A_{\text{f}}}{A_{\text{l}}}}=\mathrm{2}{C}_{\text{d}}H\mathit{\lambda }.\end{array}$$

Note that Eq. (A2) was adapted in order to be consistent with the definition of C_d 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, $\widehat{{\mathrm{\Psi }}_{\text{m}}\left(z\right)}=\mathrm{0}$ and C_d=0.25. By including the displacement height d, M98 is able to reproduce the non-linear feature of the $\frac{z_{\text{0m}}}{H}=f\left(\mathit{\lambda }\right)$ curve:

$$\begin{array}{}\text{(A14)}& z_{\text{0m,M98}}=(H-d)\exp \left(-{\left[{\displaystyle \frac{C_{\text{d}}}{{\mathit{\kappa }}^{\mathrm{2}}}}\mathit{\lambda }\left(\mathrm{1}-{\displaystyle \frac{d}{H}}\right)\right]}^{-\mathrm{0.5}}\right).\end{array}$$

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 z_0m from $\mathrm{7.6}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{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) q_low=1 and q_high=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 B2 Elevation 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.