the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Mapping the aerodynamic roughness of the Greenland Ice Sheet surface using ICESat2: evaluation over the Ktransect
Paul C. J. P. Smeets
Carleen H. Reijmer
Bert Wouters
Jakob F. Steiner
Emile J. Nieuwstraten
Walter W. Immerzeel
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 codetermine the intensity of turbulent transfer between the atmosphere and the surface. Unfortunately this parameter is often poorly known, especially in remote areas where neither highresolution elevation models nor eddycovariance 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 Ktransect (west Greenland Ice Sheet) using UAV photogrammetry, and we evaluate the modelled roughness against in situ eddycovariance observations. We then present a method to estimate the topography at 1 m horizontal resolution using the ICESat2 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 Ktransect, that is compared to an extensive dataset of in situ observations. We find a considerable spatiotemporal variability in z_{0m}, 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.
Between 1992 and 2018, the mass loss of the Greenland Ice Sheet (GrIS) contributed 10.8 ± 0.9 mm to global mean sealevel 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 lowlying ablation area of the GrIS, where bare ice is exposed to onaverage 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 HanssenBauer and Gjessing (1988) for varying seaice 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, semiempirical 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 smallscale 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 highresolution 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 eddycovariance or multilevel wind and temperature measurements. Longterm and continuous datasets remain scarce on the GrIS, although simplifying in situ methods can be applied for longterm 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 ICESat2 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 eddycovariance measurements. We then evaluate the ICESat2 elevation profiles against UAV photogrammetry, and we finally apply the bulk drag model to the ICESat2 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 ICESat2 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 Ktransect, on the west GrIS.
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 semilogarithmic 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
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{zd}{{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
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}:
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
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 ICESat2 satellite laser altimeter. Besides, the true fetch footprint depends on flow parameters such as the friction velocity (u_{*}) and the boundarylayer height (Kljun et al., 2015), which are not known a priori.
Four measured elevation profiles and a highresolution 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 zeroreferencing 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 highpass 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 highpass filter removed the influence of the largescale domes.
The height of the roughness obstacles (H) is taken as
in which $\stackrel{\mathrm{\u0303}}{{\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 smallscale 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)
where w is the width of the profile, set to 15 m. This value was chosen to match the approximate ICESat2 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 ICESat2 profiles.
To summarize, a measured elevation profile is now completely defined by the height of the obstacles (H), and the frontal area index (λ), after highpass 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.
3.1 Eddycovariance measurements
Vertical propeller eddycovariance (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 eddycovariance 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 Ktransect, 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 yearround firncovered 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 eddycovariance (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 finewire 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 nearneutral 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 Ktransect 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 highresolution elevation maps are derived using a structurefrommotion workflow using UAV imagery. Two crevassed areas close to the ice edge were mapped using an eBee fixedwing UAV from Sensefly, while the area surrounding S5 was mapped using a Mavic Pro quadcopter UAV from DJI. Multiple overlapping truecolour 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 georeferenced using the information of the UAV GPS, which yields a digital elevation model (DEM) of the surface. No additional groundcontrol 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.
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 ICESat2 profiles.
3.3 ICESat2 laser altimeter
Launched in September 2018 by the National Aeronautics and Space Administration (NASA), ICESat2 (Ice, Cloud, and Elevation Satellite2) carries a laser altimeter system in nearpolar orbit (Markus et al., 2017). The altimeter relies on a photoncounting system, which in combination with both the spacecraft's position and its pointing orientation, enables the retrieval of the 3D position of individual backscattered photons (Neumann et al., 2019). Our hypothesis is that the small footprint diameter (≈ 15 m) and short alongtrack 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 ICESat2 measurements are compared against the UAV surveys are provided in Table 2. Not more than one ICESat2 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 ICESat2 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 3D photon positions on a regular alongtrack 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 alongtrack 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 highpass 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 alongtack resolution. In the following part we describe a method to produce a 1 m resolution alongtrack surface height estimation from the ATL03 raw photons signal.
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 alongtrack segment, and a median absolute difference filter is used to remove all the photon heights which deviate too much from the local ensemble median,
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, alongtrack 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 knearest neighbour, onedimensional, 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 semivariograms 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 mediumconfidence photons. If there are not enough mediumconfidence photons, we increase the search radius to half the footprint diameter (7.5 m) or even up to a footprint diameter (15 m). The lowconfidence 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 alongtrack windows and the highpass filtering using a cutoff 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 smallscale 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 ICESat2 altimeter. Furthermore, the smallscale 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
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:
in which $\stackrel{\mathrm{\u0303}}{{\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
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 ICESat2 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 ICESat2 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.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 Ktransect (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 debriscovered 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.
To verify the suitability of several drag models (see Appendix A), we use the eddycovariance 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}.
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 eddycovariance observations (Fig. 5). We attribute this to the parametrization of C_{d}, which was derived for seaice 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 ICESat2 profiles.
4.2 Height of the roughness obstacles (H) estimated from ICESat2
The estimated height of the obstacles (H) using 2 years of ICESat2 measurements (16 October 2018–6 September 2020) crossing the lower part of the Ktransect 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. Lowlying bareice areas are rougher, while the higher, firncovered 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 ICESat2 roughness statistics against UAV DEMs
Climate models and satellite altimeter corrections require information about the largerscale spatial variability of surface (aerodynamic) roughness. This motivated us to compare the roughness statistics acquired with highresolution UAV photogrammetry to the statistics estimated from the ICESat2 laser altimeter.
The elevation profile from the UAV survey in box A (Fig. 6) was already compared to the overlapping ICESat2 profiles in Fig. 4a, while H, λ, and z_{0m} are compared in Fig. 7. In box A, the UAV and ICESat2 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 ICESat2 (Fig. 7, solid black line) clearly captures all the largest obstacles and the largescale 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 lowpassfiltering 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 ICESat2 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 largescale 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 ICESat2 profiles.
Two more UAV surveys were performed in September 2019 in area B and around S5, but the overlapping ICESat2 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 ICESat2 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 ICESat2 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 ICESat2 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.
At site S5, UAV elevation profiles and eddycovariance measurements are available in September 2019, while the ICESat2 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 ICESat2 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 ICESat2, 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 ICESat2 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 ICESat2 profiles must thus be interpreted as the aerodynamic roughness in the wind direction along the direction of the ground laser track.
Only a highresolution, twodimensional 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 eddycovariance 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).
To summarize, three independent but colocated methods, namely UAV photogrammetry, ICESat2 laser altimetry, and in situ eddycovariance 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 ICESat2 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 ICESat2 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 spatiotemporal variability of the aerodynamic roughness length over the GrIS.
4.4 Results: mapping the roughness length z_{0m} using ICESat2
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 ICESat2 ATL03 data to model and map the aerodynamic roughness (z_{0m}) over the Ktransect. We process nearly 2 years of ICESat2 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^{−3} and 10^{−1} 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 ICESat2 in two groups, July–September and October–June, which correspond to ICESat2 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^{−4} m above 1000 m a.s.l., regardless of the time period. When the ICESat2 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 landterminating.
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 Ktransect (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).
Unfortunately, the ICESat2 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 ICESat2 footprint diameter (≈ 15 m). Higher up, the ice hummocks become even smaller and the surface eventually becomes snowcovered yearround. 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 ICESat2 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 ICESat2 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 ICESat2 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, smallerscale obstacles could be retrieved in multiple directions at crossover points, using the information from multiple ICESat2 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 lowlying 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.
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 ICESat2 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 interobstacle sheltering.
Mapping surface obstacles at 1 m resolution using the ICESat2 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 ICESat2 footprint diameter of ≈ 15 m, such as ice hummocks found above 1000 m elevation in summer, or snow sastrugi expected yearround at even higher locations on the ice sheet from photographic evidence, are not resolved by the ICESat2 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 ICESat2 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^{−1} m, 2 orders of magnitude larger than typically used in regional climate models over bare ice.
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
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 seaice pressure ridges, we choose the following parameterization,
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
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
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
where c=0.25 is an empirical constant that determines the sheltering efficiency. The latter equation may be written in the form
with
which is solved iteratively using ${X}_{i+\mathrm{1}}={e}^{{X}_{i}}$ and X_{0}=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 semilogarithmic 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:
with c_{d1}=7.5, which is then used to derive the value for Ψ_{r} at height z=H using the following expression:
where
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 highpass filtered using a cutoff 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 ICESat2 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).
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:
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 nonlinear feature of the $\frac{{z}_{\text{0m}}}{H}=f\left(\mathit{\lambda}\right)$ curve:
B1 Cutoff wavelength Λ
We find that the optimal value of the cutoff wavelength for the highpass 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 ICESat2 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 lowconfidence 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 smallscale features present in the UAV data, varying each parameter does not give improved results (Fig. B2).
The data used in the figures can be downloaded here: https://doi.org/10.5281/zenodo.4386867 (Van Tiggelen, 2020). Additional data can be obtained from the authors without conditions.
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 eddycovariance data. EJN, JFS, and WI made the UAV elevation models. MvT and BW processed the ICESat2 data. MvT, with the help of all the coauthors, designed the study and wrote the first draft. All authors participated in the review and in the editing of the manuscript.
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 einfrastructure 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.
This research has been supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant no. ALWOP.431).
This paper was edited by Benjamin Smith and reviewed by Christof Lüpkes, Shane Grigsby, and one anonymous referee.
Andreas, E. L.: Airice drag coefficients in the western Weddell Sea 2. A model based on form drag and drifting snow, J. Geophys. Res., 100, 4833–4843, https://doi.org/10.1029/94JC02016, 1995. a, b, c
Arya, S. P. S.: A Drag Partition Theory for Determining the LargeScale 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 ICESat2 Ice Sheet Surface Heights, Based on Comparisons Over the Interior of the Antarctic Ice Sheet, Geophys. Res. Lett., 46, 13072–13078, https://doi.org/10.1029/2019GL084886, 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, https://doi.org/10.1029/2019jf005167, 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, https://doi.org/10.5194/tc1519312021, 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 insitu observations, Front. Earth Sci., 4, 1–9, https://doi.org/10.3389/feart.2016.00082, 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], https://doi.org/10.5194/essd202180, in review, 2021. a
Fitzpatrick, N., Radić, V., and Menounos, B.: A multiseason investigation of glacier surface roughness lengths through in situ and remote observation, The Cryosphere, 13, 1051–1071, https://doi.org/10.5194/tc1310512019, 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, https://doi.org/10.1023/A:1000466827210, 1997. a
Garbrecht, T., Lüpkes, C., Augstein, E., and Wamser, C.: Influence of a sea ice ridge on lowlevel airflow, J. Geophys. Res.Atmos., 104, 24499–24507, https://doi.org/10.1029/1999JD900488, 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, https://doi.org/10.1034/j.16000870.2002.01253.x, 2002. a, b
Gardner, C. S.: Target signatures for laser altimeters: an analysis, Appl. Opt., 21, 448, https://doi.org/10.1364/ao.21.000448, 1982. a
Garratt, J. R.: The atmospheric boundary layer, Cambridge University Press, Cambridge, 1992. a, b
HanssenBauer, I. and Gjessing, Y. T.: Observations and model calculations of aerodynamic drag on sea ice in the Fram Strait, Tellus A, 40, 151–161, https://doi.org/10.3402/tellusa.v40i2.11789, 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, https://doi.org/10.1007/s1054600691456, 2007. a, b, c
Heinemann, G.: The KABEG'97 field experiment: An aircraftbased study of Katabatic wind dynamics over the Greenland ice sheet, Bound.Lay. Meteorol., 93, 75–116, https://doi.org/10.1023/A:1002009530877, 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 fastmoving glaciers from ICESat2 laser altimeter data  Application of the densitydimension algorithm (DDA ice) and evaluation using airborne altimeter and Planet SkySat data, Sci. Remote Sens., 3, 104743, https://doi.org/10.1016/j.srs.2020.100013, 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 Highresolution monitoring of Himalayan glacier dynamics using unmanned aerial vehicles, Remote Sens. Environ., 150, 93–103, https://doi.org/10.1016/j.rse.2014.04.025, 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 groundbased image networks, Earth Surf. Process. Land., 39, 1413–1420, https://doi.org/10.1002/esp.3609, 2014. a
Kean, J. W. and Smith, J. D.: Form drag in rivers due to smallscale natural topographic features: 2. Irregular sequences, J. Geophys. Res.Earth Surf., 111, 1–15, https://doi.org/10.1029/2006JF000490, 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, https://doi.org/10.1038/s4324702000012, 2020. a
Kljun, N., Calanca, P., Rotach, M. W., and Schmid, H. P.: A simple twodimensional parameterisation for Flux Footprint Prediction (FFP), Geosci. Model Dev., 8, 3695–3713, https://doi.org/10.5194/gmd836952015, 2015. a, b
Kraaijenbrink, P. D., Shea, J. M., Pellicciotti, F., Jong, S. M., and Immerzeel, W. W.: Objectbased analysis of unmanned aerial vehicle imagery to map and characterise surface features on a debriscovered glacier, Remote Sens. Environ., 186, 581–595, https://doi.org/10.1016/j.rse.2016.09.013, 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 Ktransect on the western Greenland Ice Sheet: Surface energy balance (2003–2016), Arctic, Antarct. Alp. Res., 50, S100003, https://doi.org/10.1080/15230430.2017.1420952, 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, https://doi.org/10.1109/TGRS.2008.916639, 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, https://doi.org/10.5194/tc88012014, 2014. a, b
Lettau, H.: Note on Aerodynamic RoughnessParameter Estimation on the Basis of RoughnessElement 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, https://doi.org/10.1002/qj.3839, 2020. a
Lüpkes, C. and Gryanik, V. M.: A stabilitydependent parametrization of transfer coefficients formomentum and heat over polar sea ice to be used in climate models, J. Geophys. Res., 120, 552–581, https://doi.org/10.1002/2014JD022418, 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, https://doi.org/10.1029/2012JD017630, 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, https://doi.org/10.1002/grl.50081, 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, https://doi.org/10.1016/S13522310(97)004032, 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 Satellite2 (ICESat2): Science requirements, concept, and implementation, Remote Sens. Environ., 190, 260–273, https://doi.org/10.1016/j.rse.2016.12.029, 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, https://doi.org/10.1023/A:1000463626745, 1997. a, b
Miles, E. S., Steiner, J. F., and Brun, F.: Highly variable aerodynamic roughness length (z0) for a hummocky debriscovered glacier, J. Geophys. Res.Atmos., 122, 8447–8466, https://doi.org/10.1002/2017JD026510, 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.: Fortysix years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244, https://doi.org/10.1073/pnas.1904242116, 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., RamosIzquierdo, 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, https://doi.org/10.1016/j.rse.2019.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, https://doi.org/10.1126/sciadv.aaw0123, 2019. a
Nolin, A. W. and Mar, E.: Arctic sea ice surface roughness estimated from multiangular reflectance satellite imagery, Remote Sens., 11, 1–12, https://doi.org/10.3390/rs11010050, 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, https://doi.org/10.1002/2017JF004209, 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, https://doi.org/10.7910/DVN/OHHUKH, 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 zeroplane displacement as functions of canopy height and area index, Bound.Lay. Meteorol., 71, 211–216, https://doi.org/10.1007/BF00709229, 1994. a, b
Shao, Y. and Yang, Y.: A theory for drag partition over rough surfaces, J. Geophys. Res.Earth Surf., 113, 1–9, https://doi.org/10.1029/2007JF000791, 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, https://doi.org/10.1038/s4158601918552, 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, https://doi.org/10.1007/s1054600892910, 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, BoundaryLay. Meteorol., 92, 99–121, https://doi.org/10.1023/A:1001899015849, 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 Ktransect in west Greenland: automatic weather station data (1993–2016), Arctic, Antarct. Alp. Res., 50, S100002, https://doi.org/10.1080/15230430.2017.1420954, 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 heightretrieval algorithm for NASA's ICESat2 photoncounting laser altimeter, Remote Sens. Environ., 233, 111352, https://doi.org/10.1016/j.rse.2019.111352, 2019. a
Smith, M. W., Quincey, D. J., Dixon, T., Bingham, R. G., Carrivick, J. L., IrvineFynn, T. D. L., and Rippin, D. M.: Aerodynamic roughness of glacial ice surfaces derived from highresolution topographic data, J. Geophys. Res.Earth Surf., 121, 748–766, https://doi.org/10.1002/2015JF003759, 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.: Twentyone years of mass balance observations along the Ktransect, West Greenland, Earth Syst. Sci. Data, 4, 31–35, https://doi.org/10.5194/essd4312012, 2012. a
Van den Broeke, M. R.: Characteristics of the lower ablation zone of the West Greenland ice sheet for energybalance modelling, Ann. Glaciol., 23, 7–13, https://doi.org/10.3189/s0260305500013392, 1996. a, b
Van Tiggelen, M.: mvantiggelen/TC_2020: Release 001 (Version 001), Zenodo [Data set], https://doi.org/10.5281/zenodo.4386867, 2020. a
Van Tiggelen, M., Smeets, P. C., Reijmer, C. H., and van den Broeke, M. R.: A Vertical Propeller EddyCovariance Method and Its Application to Longterm Monitoring of Surface Turbulent Fluxes on the Greenland Ice Sheet, Bound.Lay. Meteorol., 176, 441–463, https://doi.org/10.1007/s10546020005367, 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 rectangularprism roughness elements, J. Fluid Mech., 789, 127–165, https://doi.org/10.1017/jfm.2015.687, 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, https://doi.org/10.3189/172756405781812691, 2005. a, b
Zhu, X. and Anderson, W.: Turbulent flow over urbanlike fractals: Prognostic roughness model for unresolved generations, J. Turbul., 19, 995–1016, https://doi.org/10.1080/14685248.2019.1568446, 2019. a