A New Formulation and Code to Compute Aerodynamic Roughness Length for Gridded Geometry—Tested on Lidar-Derived Snow Surfaces

Abstract

The roughness of the Earth’s surface dictates the nature of air flow across it. Detailed meteorological data that are necessary to access the aerodynamic roughness ($z_0$) are not widely collected and, as such, the geometry of a surface can be used to estimate $z_0$. Here, we present a novel formulation, and the corresponding computer code, to compute $z_0$ based on the Lettau (1969) geometric approach. The new code produces a mean $z_0$, as well as a histogram of all $z_0$ values for each individual roughness element (e.g., 10 s of thousand for the 1000 × 1000 grids) discretized using watersheds, as well as directional $z_0$ diagrams, which can be matches with the wind rose for the location. The formulation includes two parameters that may optionally be applied to smooth the surface before calculating $z_0$. By calculating $z_0$ as a function of these two parameters, we demonstrate the sensitivity of the $z_0$ value to these parameter choices. Since a large portion of the Earth’s surface is snow covered during some parts of the year, and the roughness of the snow surface varies over the snow season and over space, we apply the code to three snow surface datasets. Each surface is during a different phases of the snowpack. Each surface is evaluated at two resolutions). These surfaces are: fresh snow accumulation (1 m² at 1 and 10 mm), peak accumulation (1 km² at 1 and 10 m) and ablation sun cups (25 m² at 5 and 50 mm).

1. Introduction

Within climate, meteorologic, hydrologic, and cryospheric modeling, the aerodynamic roughness length, $z_0$, is considered a control on the latent and sensible heat transfer processes due to its reliance on wind speed, direction, and underlying topographic features [1,2,3,4,5]. Originally defined for studies of wind profiles, $z_0$ is effectively the height above a surface at which the wind velocity is zero [3,5,6]. Several wind velocities at varying heights are required to quantify $z_0$ using the wind profile method [3,7]. The measured wind velocities are then plotted against the natural log of the height, resulting in the natural log of $z_0$ as the y-intercept of the best linear fit [1].

Calculating a site-specific $z_0$ value has historically been challenging due to the reliance on a complex, multi-level wind tower to measure the wind profile [7,8]. These traditional methods, while effective, are labor-intensive, spatially limited, and provide sparse data points. However, advancements in remote sensing technologies have significantly improved the feasibility and accuracy of estimating $z_0$ over large surface areas [9]. Remote sensing techniques, particularly LiDAR [10,11] and more recently photogrammetry [12,13], can integrate the detailed topographic surface [9,14]. This geometrically based approach enables the calculation of $z_0$ values across extensive and complex landscapes, providing higher spatial resolution and broader coverage than previously possible. This new methodology of $z_0$ calculation may supplement or replace traditional techniques of determining $z_0$.

Aerodynamic roughness length is a calculated parameter which characterizes the loss of wind momentum due to surface roughness [5,15]. For this study, the well-known formula developed by Lettau in 1969 [6] is utilized:

(1)$z_0=0.5h{\displaystyle \frac{s}{S}},$

where h is the “average vertical extent”, or “effective obstacle height”, measured in cm, s is the “silhouette area of the average obstacle”, measured in ${\mathrm{cm}}^2$, and S is the “specific area”, or “lot area” measured in ${\mathrm{cm}}^2$ [6]. Figure 1 shows these physical parameters for an idealized case of simple pyramids.

This paper presents a new method, and the associated code [16], to compute $z_0$ for a surface, provided that the data are on a regular grid. Various papers (e.g., [7,9]) use the Lettau [6] method to compute $z_0$; however, the exact method of discretizing a surface into roughness elements is not specified. Other methods use or assume regularly spaced roughness elements [17,18]. Here, we use watersheds to segment specific elements of whatever size and shape (Figure 1) and then compute $z_0$ for each element. We also include filtering options to smooth a surface that allows the user to reduce the total number of elements, such as elements on top of elements.

The code can be applied to any surface and at any resolution, and produces the distribution of $z_0$ values with directionality to consider variable wind directions. Wind direction comes into play in the Formula (1) in the determination of the silhouette area s. Here, we evaluate the code on three different surfaces, each at two resolutions. The three surfaces are initially assessed for an area of 1000 by 1000 pixels and then the resolution is coarsened by an order of magnitude through averaging to a 100 by 100 pixel grid.

2. Methodology

In this section, we describe the method of computing $z_0$ from gridded data, such as those that can be obtained from LiDAR. The $z_0$ formulation is not new, but the approach is novel. This code takes as input elevation data measured on a regular grid along with a vector indicating the direction the wind is blowing. Since the surface is comprised of obstacles that are non-uniform in shape and size, we approach the $z_0$ computation element-wise. That is, we consider each obstacle as an individual roughness element.

To compute $z_0$ we first to segment the surface into roughness elements (the obstacles) ${\mathcal{O}}_i$, $i=1\dots N$. We then find the lot area $S_i$, silhouette area $s_i$, and height $h_i$ of each obstacle. The value of $z_0$ for the surface as given by Equation (1) is the sum:

(2)$z_0=0.5{\Sigma }_{i=1}^N{\displaystyle \frac{h_i}{N}}{\displaystyle \frac{s_i/N}{S_i/N}}=0.5{\displaystyle \frac{1}{N}}{\Sigma }_{i=1}^N{h}_i{\displaystyle \frac{s_i}{S_i}}$

We detail our approach below, and summarized it in Table 1.

2.1. Segmentation into Roughness Elements

One of the challenges when computing $z_0$ on real-world data is separating individual roughness elements when such elements are of non-uniform shape and size. We approach this problem using the watershed algorithm [19]. The watershed algorithm segments a topographic relief into catchment basins centered around each local minimum. As roughness elements correspond intuitively to local maxima of the surface, we first invert the surface. For a surface with points $(x,y,f(x,y\left)\right)$, the inverted surface consists of the set of points $(x,y,-f(x,y\left)\right)$. Maxima (minima) of the inverted surface correspond to minima (maxima) of the original surface. Hence, the watershed algorithm applied to the inverted surface yields a segmentation into elements associated to each local maximum of the original surface. These segments are our individual roughness elements, or obstacles.

Intuitively, the watershed algorithm places a marker at each local minimum (maximum in the original data), and floods the relief. The level of the flood is the same increases with time at the same rate in all catchment basins. As the level of the flood increases, eventually the floods originating at two distinct catchment basins will meet, at which point a watershed line is added. The process is continued until the entire relief is flooded and, therefore, segmented. In our case, each watershed segmentation is associated with a single local maximum, and we consider it to be a single roughness element.

2.2. Computing Areas and Heights

The vertical silhouette area of an obstacle represents the area that is upwind, i.e., the silhouette of the region hit directly by the blowing wind. First, we compute the surface normal vector for every data point using a standard bicubic interpolation in all three directions. This takes into account not just immediate neighbors, but nearby data points as well. Projecting each normal vector onto the wind direction vector and retaining only the points that result in a nonzero negative projection indicates data points on the surface that will be hit directly by the wind.

In a single watershed element, we identify all data points which comprise the region hit by the wind and will use them to determine the vertical silhouette area of the obstacle. For each data point determined to be in the upwind region, a general unit square the extent of a pixel is transformed using Rodrigues’ formula [20], so that the normal vector of the square is aligned with the normal vector associated with the data point. The square is then projected onto the subspace perpendicular to the direction of the wind. The area of the resulting parallelogram is computed and contributes to the vertical silhouette. To compute the vertical extent, we again restrict to the data points associated to upwind regions and subtract the minimum height from the maximum height. The lot area is the area of the roughness element, which is easily scaled from the number of points segmented into each watershed.

Using the above quantities, a value for $z_0$ is computed for each roughness element and averaged over the entire surface. Since there are thousands to tens of thousands of roughness elements per million pixels, the output is also a histogram of $z_0$ values. Further, while the wind direction tends to be consistent, such as seen over snow surfaces [21], the computation can performed directionally.

2.3. Smoothing of the Surface

Two optional modifications may be applied to the data before separating individual roughness features of the surface via the watershed algorithm. An optional Gaussian filter may be applied to the surface (in our code) that smoothens the surface so that there are fewer roughness features associated with small watershed regions. This may help, for instance, in removing small roughness features coming from measurement errors of the surface. This filter smooths by taking a weighted average of neighboring points, where the weight is determined by a symmetric two-dimensional Gaussian (with standard deviation $\sigma$) centered at each point. This dampens or removes high frequency features, i.e., small perturbations, and acts as a low-pass filter. This is a common filter choice for smoothing topographical data ([22,23,24,25]) and is included as an optional filter in some GIS software ([26,27]). The symmetric 2-D Gaussian smoothing kernel with standard deviation $\sigma$ is given as

(3)$G(x,y)-{\displaystyle \frac{1}{2\pi {\sigma }^2}}{e}^{-{\displaystyle \frac{x^2+y^2}{2{\sigma }^2}}},$

where the size of the filter is given by $2\lceil 2\sigma \rceil +1$. Additionally, maxima less than a specified percentage higher than the surrounding area may be suppressed, resulting in fewer maxima and a simplified watershed profile. This is a common technique in image analysis to reduce over-segmentation ([28]).

3. Testing of the Formulation and Code on Seasonal Snow

The code [16] is tested on the snow surface since it can evolve rapidly [29], and thus, its roughness changes quickly [9,14]. Further, $z_0$ has a substantial control on snowpack evolution, redistribution, melt rates, and melt water [5,29,30,31]. The magnitude of $z_0$ varies with surface roughness; as surface roughness increases, so does $z_0$. Fresh snow on a flat, undisturbed surface can have a $z_0$ as low as 0.0002 m, reflecting its highly smooth and uniform characteristics. In contrast, the surface of a debris-covered glacier presents a heterogeneous texture, leading to $z_0$ values that typically range between 0.005 and 0.5 m [8,32]. These variations indicate the impact of surface properties on aerodynamic roughness, which directly influences snow processes. Given this wide variability, accurately determining $z_0$ for specific sites is critical to improving the precision of hydrologic and cryospheric models [14,33]. Site-specific values capture variations in land cover, topography, and environmental factors, which generic or averaged estimates cannot adequately represent [14]. Although $z_0$ for snow has been shown to vary a lot [8,14,32], changes in $z_0$ by less than a factor of 2 can cause substantially different sublimation estimates [33,34].

We care about the snow since the seasonal snowpack is a critical temporary water storage reservoir for annual water yield for many parts of the world, from the shallow snow in the plains that cover a large portion of the globe [35,36] to the deep snows in the mountains [37]. Consequently, precise modeling of snowpack dynamics is crucial for accurate estimations of snowmelt [8]. Additionally, snow is the most reflective naturally occurring surface type on Earth and it covers a significant portion of the globe, including as much as 50% of the land mass of the Northern Hemisphere each season. In its pure state, snow reflects as much as 98% of incident solar radiation [38,39]. Rough snow absorbs more solar energy than pure snow, which affects both local melt rates as well as the global climate energy budget [40]. As the snow surface roughness increases, it absorbs more solar energy, increasing the snow melt rate and the amount of energy retained in the thermodynamic system. As such, the snowpack surface plays an important role in the Earth’s total energy budget through reflection and absorption of solar energy [31,41]. Accurately measuring snow surface roughness is essential to understanding this interaction with solar radiation and predicting its effects on snowmelt and climate.

4. Datasets and Preparation

4.1. Snow Surface Datasets

The $z_0$ value was computed for three resolutions of data representing different time periods of the snowpack evolution (accumulation, peak, and ablation) (Table 2 and Figure 2). The fresh snow (FS-FC) dataset was collected during a new snowfall event with snowcover existing on the ground. The FS-FC raw point cloud was at a resolution finer than 1 mm. The ablation-sun cup (SC-PH) dataset were collected when substantial sun cups were present during snowmelt. The SC-PH raw point cloud was at a resolution of about 3 mm. Both datasets were collected with a Faro Focus3D X 130 model Terrestrial lidar Scanner (TLS) (https://www.faro.com/ (accessed on 3 April 2025)). The ablation-sun cup data were collected as three separate scans due to the nature of the surface, i.e., sun cups. The three point cloud were sewn together with Cloud Compare, an open-source software for 3D point cloud and mesh analysis ([42]; https://www.danielgm.net/cc/ (accessed on 3 April 2025)).

The peak accumulation (PA-NS) data were from the Niwot Saddle behind Boulder, Colorado U.S.A., that is part of the Niwot Long Term Ecological Research Program (https://nwt.lternet.edu/ (accessed on 3 April 2025)). The point clouds were not available publicly, and the data were already interpolated to a 1-m grid. The data were obtained as per Harpold et al. [44]. All data are available online [45].

The elevation range is larger for the more coarse surfaces (Figure 2d vs. Figure 2g and Figure 2a). However, comparing the elevation range (Figure 2c,f,i) to the resolution illustrates less prominence in roughness features as the resolutions becomes larger. Coarsening the resolution of the individual datasets did not alter the elevation distribution, except for the sun cup surface (Figure 2i). These surfaces present three different snowpack conditions with features highlighted by the varying resolutions.

4.2. Data Preparation

Following the collection of the lidar scans (fresh snow and ablation-sun cup datasets), each raw point cloud (three sewn together for SC-PH) was processed and evaluated using Cloud Compare. The scans were then cropped to the area of interest, ensuring that the dataset was limited to the preferred geographic region or feature under study. Additionally, any rogue points, such as noise or erroneous data points often caused by environmental elements, were identified and removed to enhance the accuracy and clarity of the dataset. The cropped point cloud was imported into Golden Software’s Surfer ([46]; https://www.goldensoftware.com/products/surfer (accessed on 3 April 2025)). The software allows the user to interpolate the data using pre-defined methodology, such as kriging, used for this study, and user defined resolutions. The kriging method was chosen because of its capacity to account for the spatial arrangement of measured points, which is relevant to the surface roughness analysis. The result is a smoothed, gridded, surface, instead of the point cloud, which is required to be used within the code.

Each gridded dataset was detrended in the X and Y directions to reduce bias from sloping surfaces [43]. Since the Niwot Saddle covered a hillslope with substantial terrain variation (see [47]), the surface (PA-NS) was detrended three times. First, there was a linear detrend to removed the overall slope. Then, a 3rd-order polynomial was fitted to the surface to remove the saddle feature. Finally, a linear detrend removed the positioning bias.

5. Results

The value of $z_0$ varied as a function of resolution (Table 3). This is both for the resolution of the original data, i.e., fresh snow versus peak accumulation versus ablation sun cups, and coarsening of the data. There are fewer elements as we coarsen the data, and thus, fewer elements become more granular. For the finest resolution surface (fresh snow), the elements are the smallest.

Most of the roughness elements were small and thus had a low $z_0$ value (a left skew in Figure 3). For the finer resolution grids (Figure 3a,c,e), more than 80% of the $z_0$ value were in the lowest bin. The sun cup surface had the most undulations (Figure 2c,d), and the most elements at the original and coarsened resolution (Table 3). At smaller scales, there are many more features with a low $z_0$ value.

The code computes $z_0$ with directionality, which does vary as a function of coarsening (Figure 4). The smoothest (lowest $z_0$ values) are not the same direction between the original and coarsened datasets. This is relevant because, while wind tends to blow consistently from the same direction due to terrain influences [48], it does vary [49].

As discussed in Section 2.3, two optional parameters, $\sigma$ and % of maximum height, may be chosen to smoothen the surface before calculating $z_0$. For the $z_0$ values reported in Table 3 and Figure 3 and Figure 4, no smoothing was implemented (both smoothing parameters were set to 0). To illustrate how $z_0$ depends on these parameters, we calculated $z_0$ for ranges parameter values for the Niwot 1 m dataset. Figure 5 shows the results. For small values of $\sigma$, $z_0$ is insensitive to the choice of $\sigma$, as seen in the left panel of Figure 5. For larger values of $\sigma$, the surface is smoothed to the extent that essential roughness features are lost. Figure 6 illustrates how the watersheds associated with peaks of the surface decrease as $\sigma$ increases. However, $z_0$ does depend more sensitively to the % parameter for small values of $\sigma$, as seen in the right panel of Figure 5.

6. Discussion

Lettau’s formulation [6] of the aerodynamic roughness length $z_0$, implemented in this paper, is a geometric representation of the surface. There are other formulations [17,50], but they do not appear to be as useful for snow [9]. Lettau’s formulation [6] is widely used for snow and glacier surfaces [4,7,8]. With the propensity of spatial surface measurement approaches and data, especially for the snowpack [11,12,13,44], the geometric estimation of $z_0$ is a useful approach. As such, this paper provides the computer code [16] to estimate Lettau’s [6] $z_0$ formulation.This paper is not a comparison between the anemometric and geometric approaches to compute $z_0$ [9].

Andreas [51] derived an empirical formulation of $z_0$ for sea ice based on the physical roughness of the surface. This work used transect data, based on the dominant wind direction, to compute the roughness. The Andreas [51] approach could be applied to spatial datasets to consider changes in wind direction (Figure 4). Wind directions can vary seasonally and even daily based on weather [49].

The theoretical foundation of Lettau’s formula [6] allows it to be applied to any surface type; however, this study focused on the application on snow-covered terrain. This is not due to the limitation of the method, but from the availability of the LiDAR datasets used. The high-resolution LiDAR scan made it an ideal test case for analyzing the geometry of the snow surface. Additionally, the snow surface is spatially and temporally variable in nature and, therefore, this methodology is applicable throughout an entire snow-covered season. However, this approach and code are applicable to other surface types where adequate spatial surface data are available.

Following common approaches in image analysis to reduce effects of measurement noise or oversegmentation, the code includes two optional parameters that smoothen the surface before calculating $z_0$. The first parameter is the standard deviation $\sigma$ of the smoothing filter. The second parameter is the % of the maximum height of the surface; peaks with height lower than this percentage height are not included as obstacles contributing to the $z_0$ calculation. We showed the choice of these parameter values does significantly affect the calculated value of $z_0$. However, in the case of the Niwot 1 m dataset, $z_0$ is insensitive to the choice $\sigma$ for small values of $\sigma$. The values of these smoothing parameters should be reported in any calculation of $z_0$ using the Lettau method as implemented in this paper. Further work should elucidate how to choose these parameter values to best match an application, such as using the aerodynamic roughness in a simulation of fluid flow across a surface.

Scale and resolution play a critical role in selecting an appropriate value of $z_0$, as highlighted by Sanow [14]. Since the aerodynamic roughness length is determined by the geometric characteristics of the surface, the spatial resolution of the input data must be aligned with the physical scale of the surface features. Therefore, users of this method and the accompanying code should consider matching the resolution of their data with the spatial extent and variability of the study area.

For small-scale studies, (1 × 1 m, 10 × 10 m) resolutions in the range of 1 cm to 1 mm are recommended to accurately capture the microtopography. Conversely, for larger domains (1 × 1 km and greater) a coarser resolution of >1 m may be appropriate, as finer-scale variations become less influential at broader scales. However, finer is not always better, extremely high-resolution data over large areas may introduce unnecessary computational load without significantly improving $z_0$ estimates. Ultimately, the preferred resolution depends on both the spatial scale of the study and the objectives of the analysis. For instance, if the goal is to model surface-atmosphere interactions at the microscale, high-resolution data are necessary. For regional-scale modeling, a coarser resolution that captures broader topographic patterns may suffice.

The method and code [16] presented here has broad applicability, beyond snow. We focus on snow roughness due to it dynamic nature [29]. Computing $z_0$ from surface geometry can rely on the terrestrial [10,11] and airborne lidar [44,52] data collection methods; spaceborne lidar [53,54] is also available to measure the snow and other surfaces. Further, photogrammetry and structure-from-motion methods are now available from unmanned aerial vehicles (UAVs) [13], from airplanes [12], and satellites [55]. Hybrid lidar-photogrammetry is available over small areas at a 1 to 2-cm resolution with the iPad Pro and the 3-D scanner app [56]. All these new data are mapping the snow surface at a variety of resolutions and extents and thus when gridded, they can be used to estimated the geometric-based $z_0$. We recommend that the user consider the balance between the resolution and extent of the data they collect as well as working across the different scales.

Author Contributions

Conceptualization, R.P., S.R.F., I.O., R.A.N., P.D.S. and J.E.S.; methodology, R.P., R.A.N., P.D.S., S.R.F. and J.E.S.; software, R.P., R.A.N., P.D.S. and J.E.S.; validation, R.A.N., P.D.S., J.E.S. and S.R.F.; formal analysis, R.A.N., P.D.S., J.E.S. and S.R.F.; data curation, S.R.F. and J.E.S.; writing—original draft preparation, R.P., S.R.F., J.E.S., R.A.N. and P.D.S.; writing—review and editing, S.R.F., R.A.N., P.D.S. and J.E.S.; visualization, R.A.N., P.D.S., S.R.F. and R.P.; supervision, S.R.F., P.D.S. and I.O.; project administration, S.R.F., I.O. and P.D.S.; funding acquisition, I.O., S.R.F. and P.D.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The Peak Accumulation snow surface from Niwot Saddle is available from the Niwot Long Term Ecological Research Program (https://nwt.lternet.edu/ (accessed on 3 April 2025)), and were identified in Harpold et al. [44]. The Fresh Snow, Peak Accumulation, and Sun Cups snow surfaces are available through Dyrad [45] (https://doi.org/10.5061/dryad.qv9s4mwr2 (accessed on 3 April 2025)). The code presented in this paper is available at [16].

Acknowledgments

Graham Sexstone of the U.S. Geological Survey acquired the Niwot Sadddle dataset for this work.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

AOI	is the area of interest h is the average vertical extent or effective obstacle height, measured in cm
s	is the silhouette area of the average obstacle, measured in ${\mathrm{cm}}^2$
S	is the specific area or lot area, measured in ${\mathrm{cm}}^2$
$z_0$	is the (geometric) aerodynamic roughness length, measured in m

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Physical parameters used in Lettau’s $z_0$ formula given by Equation (1). AOI is the Area Of Interest.

Figure 2 The three gridded datasets used in the $z_0$ computation as 1000 by 1000 pixels (a,d,g) and coarsened by a factor of 10 (b,e,h), plus the histograms of the distribution of elevation for each surface (c,f,i) for (a–c) fresh snow in Fort Collins (lines in 0.005 m intervals), (d–f) peak accumulation at Niwot Saddle (lines in 2 m intervals), and (g–i) ablation sun cups in the Poudre Headwaters (lines in 0.05 m intervals). The surfaces have been detrended to remove any slope-based bias [43]. As such, black lines/areas are above the datum of zero while red lines/areas are below the datum of zero.

Figure 3 Histogram of $z_0$ values for the peak accumulation on fresh snow in Fort collins at (a) 0.001 m and (b) 0.01 m resolutions, the Niwot Saddle (c) at 1 m and (d) 10 m resolutions, and the ablation sun cups at the Poudre Headwaters at (e) 0.005 m and (f) 0.05 m resolutions for an eastward wind.

Figure 4 The mean $z_0$ value varies with the direction from which the wind hits the surface. Shown above is the mean $z_0$ as it varies with wind direction from 24 different directions for the peak accumulation on fresh snow in Fort Collins at (a) 0.001 m and (b) 0.01 m resolutions, the Niwot Saddle (c) at 1m and (d) 10 m resolutions, the ablation sun cups at the Poudre Headwaters at (e) 0.005 m and (f) 0.05 m resolutions.

Figure 5 Shown is the sensitivity of the mean $z_0$ of the Niwot Saddle at 1m resolutions for an eastward wind to changes in the parameter $\sigma$ (left) in the Gaussian smoothing and to changes in both the % peak suppression parameter and to $\sigma$ (right) in the preliminary filtering step.

Figure 6 The watershed of a small subset of the Niwot Saddle at 1m resolution for the (a) raw data and data filtered via Gaussian smoothing with (b) $\sigma =0.5$, (c) $\sigma =1$, and (d) $\sigma =2$.

Table 1

Summary of the encoded algorithm for computing $z_0$ for a single wind direction. The user must specify the dataset, scale, smoothening parameters (if applicable), and directions of the wind. All wind directions may be optionally explored. When looping over all wind directions, the smoothening and watershed computations are consistent for all directions of the wind.

Step		Details
Smoothing (Optional)		· Gaussian
		· Small maxima suppression
Segment roughness elements		via watershed algorithm on the inverted surface
For each obstacle	Determine lot area	from the extent of the watershed
	Determine silhouette area	restricting to the upwind component of obstacle
	Determine height	restricting to the upwind component of obstacle
	Compute $z_0$
Compute summary statistics		over all obstacles (for a single single wind direction)

Table 2

Summary of the three datasets used in the $z_0$ computation. All are for a 1000 by 1000 pixel area at the stated resolution. Data were obtained in the state of Colorado, USA.

Dataset	Fresh Snow	Peak Accumulation	Ablation—Sun Cups
short name/code	FS-FC	PA-NS	SC-PH
location	Fort Collins	Niwot Saddle	Poudre headwaters
latitude, longitude	40.6, −105.1	40.0547, −105.5890	40.4396, −105.7739
date acquired	19 April 2021	20 May 2010	13 June 2017
resolution (m)	0.001	1	0.05

Table 3

Number of elements identified for each datasets and resolution with the mean $z_0$ value (in cm) for wind blowing from west to east.

Surface	Fresh Snow	-	Peak Accumulation	-	Ablation—Sun Cups	-
resolution (m)	0.001	0.010	1	10	0.005	0.050
elements	7366	195	21090	74	36259	380
mean $z_0$ ($\times {10}^{-3}$)	0.00054	0.0075	9.4	421	0.012	0.14

References

1. Oke, T.R. Boundary Layer Climates; 2nd ed. Cambridge University Press: Cambridge, UK, 1987; ISBN 0-415-04319-0

2. Andreas, E. Parameterizing scalar transfer over snow and ice: A review. J. Hydrometeorol.; 2002; 3, pp. 417-432. [DOI: https://dx.doi.org/10.1175/1525-7541(2002)003<0417:PSTOSA>2.0.CO;2]

3. Gromke, C.; Manes, C.; Walter, B.; Lehning, M.; Guala, M. Aerodynamic roughness length of fresh snow. Bound.-Layer Meteorol.; 2011; 141, pp. 21-34. [DOI: https://dx.doi.org/10.1007/s10546-011-9623-3]

4. Nield, J.M.; King, J.; Wiggs, G.F.; Leyland, J.; Bryant, R.G.; Chiverrell, R.C.; Darby, S.E.; Eckardt, F.D.; Thomas, D.S.; Vircavs, L.H. . Estimating aerodynamic roughness over complex surface terrain. J. Geophys. Res. Atmos.; 2013; 118, pp. 12948-12961. [DOI: https://dx.doi.org/10.1002/2013JD020632]

5. Smith, M.W. Roughness in the earth sciences. Earth Sci. Rev.; 2014; 136, pp. 202-225. [DOI: https://dx.doi.org/10.1016/j.earscirev.2014.05.016]

6. Lettau, H. Note on aerodynamic roughness-parameter estimation on the basis of roughness-element description. J. Appl. Meteorol.; 1969; 8, pp. 828-832. [DOI: https://dx.doi.org/10.1175/1520-0450(1969)008<0828:NOARPE>2.0.CO;2]

7. Munro, D.S. Surface roughness and bulk heat transfer on a glacier: Comparison with eddy correlation. J. Glaciol.; 1989; 35, pp. 343-348. [DOI: https://dx.doi.org/10.3189/S0022143000009266]

8. Miles, E.S.; Steiner, J.F.; Brun, F. Highly variable aerodynamic roughness length for a hummocky debris-covered glacier. J. Geophys. Res. Atmos.; 2017; 122, pp. 8447-8466. [DOI: https://dx.doi.org/10.1002/2017JD026510]

9. Sanow, J.E.; Fassnacht, S.R.; Kamin, D.J.; Sexstone, G.A.; Bauerle, W.L.; Oprea, I. Geometric versus anemometric surface roughness for a shallow accumulating snowpack. Geosciences; 2018; 8, 463. [DOI: https://dx.doi.org/10.3390/geosciences8120463]

10. Hood, J.L.; Hayashi, M. Assessing the application of a laser rangefinder for determining snow depth in inaccessible alpine terrain. Hydrol. Earth Syst. Sci.; 2010; 14, pp. 901-910. [DOI: https://dx.doi.org/10.5194/hess-14-901-2010]

11. Revuelto, J.; López-Moreno, J.I.; Azorín-Molina, C.; Zabalza, J.; Arguedas, G.; Vicente-Serrano, S.M. Mapping the annual evolution of snow depth in a small catchment in the Pyrenees using the long-range terrestrial laser scanning. J. Maps; 2014; 10, pp. 379-393. [DOI: https://dx.doi.org/10.1080/17445647.2013.869268]

12. Nolan, M.; Larsen, C.; Sturm, M. Mapping snow depth from manned aircraft on landscape scales at centimeter resolution using structure-from-motion photogrammetry. Cryosphere; 2015; 9, pp. 1445-1463. [DOI: https://dx.doi.org/10.5194/tc-9-1445-2015]

13. Revuelto, J.; Alonso-Gonzalez, E.; Vidaller-Gayan, I.; Lacroix, E.; Izagirre, E.; Rodríguez-López, G.; López-Moreno, J.I. Intercomparison of UAV platforms for mapping snow depth distribution in complex alpine terrain. Cold Reg. Sci. Technol.; 2021; 190, 103344. [DOI: https://dx.doi.org/10.1016/j.coldregions.2021.103344]

14. Sanow, J.E.; Fassnacht, S.R.; Suzuki, K. How does a dynamic surface roughness affect snowpack modelling?. Polar Sci.; 2024; 41, 101110. [DOI: https://dx.doi.org/10.1016/j.polar.2024.101110]

15. Jacobson, M.Z. Fundamentals of Atmospheric Modeling; 2nd ed. Cambridge University Press: Cambridge, UK, 2005; ISBN 978-0-521-54865-6

16. z₀ Lettau LiDAR Watershed Code. [GitHub Repository]. 2025; Available online: https://github.com/rneville/Z0-Lettau-LiDAR-Watershed (accessed on 4 June 2025).

17. Counihan, J. Wind tunnel determination of the roughness length as a function of the fetch and the roughness density of three-dimensional roughness elements. Atmos. Environ.; 1971; 5, pp. 637-642. [DOI: https://dx.doi.org/10.1016/0004-6981(71)90120-X]

18. Petersen, R.L. A wind tunnel evaluation of methods for estimating surface roughness length at industrial facilities. Atmos. Environ.; 1997; 31, pp. 45-57. [DOI: https://dx.doi.org/10.1016/S1352-2310(96)00154-9]

19. Meyer, F. Topographic distance and watershed lines. Signal Process.; 1994; 38, pp. 113-125. [DOI: https://dx.doi.org/10.1016/0165-1684(94)90060-4]

20. Encyclopedia of Mathematics: Rodrigues Formula. Available online: https://encyclopediaofmath.org/wiki/Rodriguesformula (accessed on 30 December 2024).

21. Winstral, A.; Elder, K.; Davis, R.E. Spatial snow modeling of wind-redistributed snow using terrain-based parameters. J. Hydrometeorol.; 2002; 3, pp. 524-538. [DOI: https://dx.doi.org/10.1175/1525-7541(2002)003<0524:SSMOWR>2.0.CO;2]

22. Lindsay, J.B. Whitebox GAT: A case study in geomorphometric analysis. Comput. Geosci.; 2016; 95, pp. 75-84. [DOI: https://dx.doi.org/10.1016/j.cageo.2016.07.003]

23. Lindsay, J.B.; Francioni, A.; Cockburn, J.M.H. LiDAR DEM Smoothing and the Preservation of Drainage Features. Remote Sens.; 2016; 11, 1926. [DOI: https://dx.doi.org/10.3390/rs11161926]

24. O’Neil, G.L.; Saby, L.; Band, L.E.; Goodall, J.L. Effects of LiDAR DEM smoothing and conditioning techniques on a topography-based wetland identification model. Water Resour. Res.; 2019; 55, pp. 4343-4363. [DOI: https://dx.doi.org/10.1029/2019WR024784]

25. Erdbrügger, J.; van Meerveld, I.; Bishop, K.; Seibert, J. Effect of DEM-smoothing and-aggregation on topographically-based flow directions and catchment boundaries. J. Hydrol.; 2021; 602, 126717. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126717]

26. Whitebox Tools. Available online: https://www.whiteboxgeo.com/ (accessed on 4 June 2025).

27. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Böhner, J. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev.; 2015; 8, pp. 1991-2007. [DOI: https://dx.doi.org/10.5194/gmd-8-1991-2015]

28. The MathWorks, Inc. Imhin Documentation. mathworks.com. Available online: https://www.mathworks.com/help/images/ref/imhmin.html (accessed on 5 October 2024).

29. Fassnacht, S.R.; Williams, M.W.; Corrao, M.V. Changes in the surface roughness of snow from millimetre to metre scales. Ecol. Complex.; 2009; 6, pp. 221-229. [DOI: https://dx.doi.org/10.1016/j.ecocom.2009.05.003]

30. Andreas, E.L. A theory for the scalar roughness and the scalar transfer coefficients over snow and sea ice. Bound.-Layer Meteorol.; 1987; 38, pp. 159-184. [DOI: https://dx.doi.org/10.1007/BF00121562]

31. Musselman, K.N.; Clark, M.P.; Liu, C.; Ikeda, C.; Rasmussen, R. Slower snowmelt in a warmer world. Nat. Clim. Chang.; 2017; 7, pp. 214-220. [DOI: https://dx.doi.org/10.1038/nclimate3225]

32. Brock, B.; Willis, I.; Sharp, M. Measurement and parameterization of aerodynamic roughness length variations at Haut Glacier d’Arolla, Switzerland. J. Glaciol.; 2006; 52, pp. 281-297. [DOI: https://dx.doi.org/10.3189/172756506781828746]

33. Fassnacht, S.R. Temporal changes in small scale snowpack surface roughness length for sublimation estimates in hydrological modelling. Cuad. Investig. Geog.; 2010; 36, 43. [DOI: https://dx.doi.org/10.18172/cig.1226]

34. Hultstrand, D.M.; Fassnacht, S.R. The Sensitivity of Snowpack Sublimation Estimates to Instrument and Measurement Uncertainty Perturbed in a Monte Carlo Framework. Front. Earth Sci.; 2018; 12, pp. 728-738. [DOI: https://dx.doi.org/10.1007/s11707-018-0721-0]

35. Brown, R.D. Northern Hemisphere Snow Cover Variability and Change, 1915–97. J. Clim.; 2000; 13, pp. 2339-2355. [DOI: https://dx.doi.org/10.1175/1520-0442(2000)013<2339:NHSCVA>2.0.CO;2]

36. Mialon, A.; Fily, M.; Royer, A. Seasonal snow cover extent from microwave remote sensing data: Comparison with existing ground and satellite based measurements. EARSeL EProceedings; 1994; 4, pp. 215-225.

37. Immerzeel, W.W.; Lutz, A.F.; Andrade, M.; Bahl, A.; Biemans, H.; Bolch, T.; Hyde, S.; Brumby, S.; Davies, B.J.; Elmore, A.C. . Importance and vulnerability of the world’s water towers. Nature; 2020; 577, pp. 364-369. [DOI: https://dx.doi.org/10.1038/s41586-019-1822-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31816624]

38. Kukko, A.; Anttila, K.; Manninen, T.; Kaasalainen, S.; Kaartinen, H. Snow surface roughness from mobile laser scanning data. Cold Reg. Sci. Technol.; 2013; 96, pp. 23-35. [DOI: https://dx.doi.org/10.1016/j.coldregions.2013.09.001]

39. Zhuravleva, T.B.; Kokhanovsky, A. Influence of surface roughness on the reflective properties of snow. J. Quant. Spectrosc. Radiat. Transf.; 2011; 112, pp. 1352-1368. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2011.01.004]

40. Amory, C.; Naaim-Bouvet, F.; Gallee, H.; Vignon, E. Brief communication: Two well-marked cases of aerodynamic adjustment of sastrugi. Cryosphere; 2016; 10, pp. 743-750. [DOI: https://dx.doi.org/10.5194/tc-10-743-2016]

41. Thackeray, C.W.; Fletcher, C.G. Snow albedo feedback: Current knowledge, importance, outstanding issues and future directions. Prog. Phys. Geogr. Earth Environ.; 2016; 40, pp. 392-408. [DOI: https://dx.doi.org/10.1177/0309133315620999]

42. CloudCompare. Available online: https://www.danielgm.net/cc/ (accessed on 4 June 2025).

43. Fassnacht, S.R.; Stednick, J.D.; Deems, J.S.; Corrao, M.V. Metrics for assessing snow surface roughness from digital imagery. Water Resour. Res.; 2009; 45, W00D31. [DOI: https://dx.doi.org/10.1029/2008WR006986]

44. Harpold, A.A.; Guo, Q.; Molotch, N.; Brooks, P.D.; Bales, R.; Fernandez-Diaz, J.C.; Musselman, K.N.; Swetnam, T.L.; Kirchner, P.; Meadows, M.W. . LiDAR-derived snowpack data sets from mixed conifer forests across the Western United States. Water Resour. Res.; 2024; 50, pp. 2749-2755. [DOI: https://dx.doi.org/10.1002/2013WR013935]

45. Fassnacht, S.R.; Sanow, J.E. Fresh Snow and Ablation-Sun Cup Snow Surface Datasets for Evaluation of Geometry Aerodynamic Roughness Code [Dataset]. Dryad; 2025. [DOI: https://dx.doi.org/10.5061/dryad.qv9s4mwr2]

46. Golden Software—Surfer. Available online: https://www.goldensoftware.com/products/surfer/ (accessed on 30 December 2024).

47. Sexstone, G.A.; Fassnacht, S.R.; López-Moreno, J.I.; Hiemstra, C.A. Subgrid snow depth coefficient of variation spanning alpine to sub-alpine mountainous terrain. Cuad. Investig. Geográfica/Geogr. Res. Lett.; 2022; 48, pp. 79-96. [DOI: https://dx.doi.org/10.18172/cig.4951]

48. Winstral, A.; Marks, D. Simulating wind fields and snow redistribution using terrain-based parameters to model snow accumulation and melt over a semi-arid mountain catchment. Hydrol. Processes; 2002; 16, pp. 3585-3603. [DOI: https://dx.doi.org/10.1002/hyp.1238]

49. Moron, V.; Robertson, A.W.; Qian, J.-H.; Ghil, M. Weather types across the Maritime Continent: From the diurnal cycle to interannual variations. Front. Environ. Sci.; 2015; 2, 65. [DOI: https://dx.doi.org/10.3389/fenvs.2014.00065]

50. Macdonald, R.W.; Griffiths, R.F.; Hall, D.J. An improved method for the estimation of surface roughness of obstacle arrays. Atmos. Environ.; 1998; 32, pp. 1857-1864. [DOI: https://dx.doi.org/10.1016/S1352-2310(97)00403-2]

51. Andreas, E.L. A relationship between the aerodynamic and physical roughness of winter sea ice. Q. J. R. Meteorol. Soc.; 2011; 137, pp. 1581-1588. [DOI: https://dx.doi.org/10.1002/qj.842]

52. Painter, T.H.; Berisford, D.F.; Boardman, J.W.; Bormann, K.J.; Deems, J.S.; Gehrke, F.; Hedrick, A.; Joyce, M.; Laidlaw, R.; Marks, D. . The Airborne Snow Observatory: Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo. Remote Sens. Environ.; 2016; 184, pp. 139-152. [DOI: https://dx.doi.org/10.1016/j.rse.2016.06.018]

53. Singh, U.M.; Ismail, S.; Kavaya, M.J.; Winker, D.M.; Amzajerdian, F. Space-Based Lidar. Laser Remote Sensing; 1st ed. CRC Press: Boca Raton, FL, USA, 2005; pp. 799-900.

54. Kwon, Y.; Yoon, Y.; Forman, B.A.; Kumar, S.V.; Wang, L. Quantifying the observational requirements of a space-borne LiDAR snow mission. J. Hydrol.; 2021; 601, 126709. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126709]

55. Marti, R.; Gascoin, S.; Berthier, E.; de Pinel, M.; Houet, T.; Laffly, D. Mapping snow depth in open alpine terrain from stereo satellite imagery. Cryosphere; 2016; 10, pp. 1361-1380. [DOI: https://dx.doi.org/10.5194/tc-10-1361-2016]

56. Fassnacht, S.R.; Herrero, J.; Sanow, J.E. Consequences of Snow Surface Roughness Variability for Different Melting Conditions. Proceedings of the 92nd Annual Western Snow Conference; Bozeman, MT, USA, 19–22 May 2025.