A Geometric-based Aerodynamic Roughness Length Formulation and Code Applied to Lidar-derived Snow Surfaces

1. Introduction

Shallow, seasonal snowpack is a critical reservoir for yearly water yield in the western US [1,2]. The snowpack acts as a temporary water storage system and is essential for all regional hydrology, influencing water availability for all natural ecosystems, agricultural, urban systems. Consequently, precise modeling of snowpack dynamics is crucial for accurate estimations of snowmelt [3]. Within 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 [4,5,6]. Because of this influence, $z_0$ has a substantial control on snowpack evolution, redistribution, melt rates, and melt water [2,6,7,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 [9,10]. Rough snow absorbs more solar energy than pure snow, which affects both local melt rates as well as the global climate energy budget [11]. 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 [2,12]. Accurately measuring snow surface roughness is essential to understanding this interaction with solar radiation and predicting its effects on snowmelt and climate.

Originally defined for studies of wind profiles, $z_0$ is effectively the height above a surface at which the wind velocity is zero [5,6,13]. Several wind velocities at varying heights are required to quantify $z_0$ using the wind profile method [5,14]. 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. 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 meters, 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 meters [2,15]. 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 [16,17]. Site-specific values capture variations in land cover, topography, and environmental factors, which generic or averaged estimates cannot adequately represent [17].

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 [2,14]. 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 [18]. Remote sensing techniques, particularly LiDAR [19] and more recently photogrammetry [20], can integrate the detailed topographic surface and the snowpack surface properties [17,18]. 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 [6,21]. For this study, the well-known formula developed by Lettau in 1969 [13] is utilized:

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

(1)

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 cm², and S is the “specific area", or “lot area” measured in cm² [13]. Figure 1 shows these physical parameters for an idealized case of simple pyramids.

This paper presents code to compute $z_0$ for a surface, provided that the data are on a regular grid. 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 use three different snow surfaces at different resolutions, as the snow surface, and thus its roughness, can evolve rapidly [8]. We use a fresh snow surface at 1 mm resolution, a snow surface at the time of peak accumulation at 1 m resolution, and a melting snow with sun cups at 5 mm resolution. 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$ on LiDAR data. 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 nonuniform 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\cdots 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:

$$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}}$$

(2)

We detail our approach below.

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 [22]. 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 3 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 [23], 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 histogram of $z_0$ values. Further, while the wind direction tends to be consistent, especially over a snow surface [24], 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 ([25,26,27,28]) and is included as an optional filter in some GIS software ([29,30]). The symmetric 2-D Gaussian smoothing kernel with standard deviation $\sigma$ is given as

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

(3)

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 ([31]).

3. Datasets and Preparation

3.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 1 and Figure 2). The fresh snow (FS-FC) dataset was collected during a new snowfall event with snow 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/). 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 ([32]; https://www.danielgm.net/cc/).

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/). 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. [33]. All data are available online [34].

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, b, c) fresh snow in Fort Collins (lines in 0.005 m intervals), (d, e, f) peak accumulation at Niwot Saddle (lines in 2 m intervals), and (g, h, 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 [35]. As such, black lines/areas are above the datum of zero while red lines/areas are below the datum of zero.

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, b, c) fresh snow in Fort Collins (lines in 0.005 m intervals), (d, e, f) peak accumulation at Niwot Saddle (lines in 2 m intervals), and (g, h, 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 [35]. As such, black lines/areas are above the datum of zero while red lines/areas are below the datum of zero.

The elevation range is larger for the more coarse surfaces (Figure 2d versus Figure 2g and Figure 2a). However, comparing the elevation range (Figure 2c, f and 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.

3.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 ([36]; https://www.goldensoftware.com/products/surfer). 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 [35]. Since the Niwot Saddle covered a hillslope with substantial terrain variation (see [37]), the surface (PA-NS) was detrended three times. First, there was a linear detrend to removed the overall slope. Then, a 3^rd-order polynomial was fitted to the surface to remove the saddle feature. Finally, a linear detrend removed the positioning bias.

4. Results

The value of $z_0$ varied as a function of resolution (Table 2). 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, and e), more than 80% of the $z_0$ value were in the lowest bin. The sun cup surface had the most undulations (Figure 2) c and d), and the most elements at the original and coarsened resolution (Table 2). 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 [38], it does vary [39].

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 2 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 data set. 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.

5. Discussion

Lettau’s formulation [13] of the aerodynamic roughness length $z_0$, implemented in this paper, is a geometric representation of the surface. There are other formulations [40,41], but they do not appear to be as useful for snow [18]. Lettau’s formulation [13] is widely used for snow and glacier surfaces [2,14,42]. With the propensity of spatial surface measurement approaches and data, especially for the snowpack [20,33,43,44,45], the geometric estimation of $z_0$ is a useful approach. As such, this paper provides the computer code to estimate Lettau’s [13]$z_0$ formulation.

Andreas [46] derived an empirical formulation of $z_0$ for sea ice based on the physical roughness of the surface. It used transect data, based on the dominant wind direction, to compute the roughness. The Andreas [46] 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 [39].

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 data set, $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.