Geometric Versus Anemometric Surface Roughness for a Shallow Accumulating Snowpack

When applied to a snow-covered surface, aerodynamic roughness length, z0, is typically considered as a static parameter within energy balance equations. However, field observations show that z0 changes spatially and temporally, and thus z0 incorporated as a dynamic parameter may greatly improve models. To evaluate methods for characterizing snow surface roughness, we compared concurrent estimates of z0 based on (1) terrestrial light detection and ranging derived surface geometry of the snowpack surface (geometric, z0G) and (2) vertical wind profile measurements (anemometric, z0A). The value of z0G was computed from Lettau’s equation and underestimated z0A but compared well when scaled by a factor of 2.34. The Counihan method for computing z0G was found to be unsuitable for estimating z0 on a snow surface. During snowpack accumulation in early winter, z0 varied as a function of the snow-covered area (SCA). Our results show that as the SCA increases, z0 decreases, indicating there is a topographic influence on this relation.


Introduction
In the Northern Hemisphere, a seasonal snowpack can cover over 50% of the land area with the snow surface often the interface between the atmosphere and the earth [1]. The roughness of a snow surface is an important control on air-snow heat transfer [2], and changes in the snow surface can have substantial effects on the energy balance at this interface. Snow is a complicated surface with rapidly evolving physical roughness characteristics due to changing atmospheric conditions, the metamorphism of snow crystals, melting and freezing processes and redistribution by wind, especially in open areas [3]. Roughness characteristics also influence the air-surface momentum transfer on the snowpack due to wind [4]. The changes in wind momentum can reduce the energy budget, influence the formation of roughness features, and affect the aeolian movement of snow [4]. Heat flux modeling has typically used the aerodynamic roughness length (z 0 ) as a static parameter, in hydrologic, snowpack, and climate models [5,6], with z 0 only varying as a function of land cover type. For example, the Community Land Model version 4.0 (CLM4; http://www.cesm.ucar.edu/models/ccsm4.0/clm/)

Materials and Methods
The capability of a rough surface to absorb momentum from a turbulent boundary layer can be quantified by z 0 , which is a measure of the vertical turbulence that occurs when a horizontal wind flows over a rough surface [32]. In general, z 0 is a quantity that is computed from the Reynolds number and the roughness geometry of the surface [29]. For rough, turbulent regimes occurring in the atmospheric boundary layer, dependence on the Reynolds number vanishes and z 0 is only a function of roughness geometry [33]. Various relations have been found to relate the geometry of roughness elements with z 0 [2,29]. For example, the dependence of z 0 on the size, shape, density, and distribution of surface elements has been studied using wind tunnels, analytical investigations, numerical modeling, and field observations [34,35]. Smith [36] provides a comprehensive review of the different approaches and models developed to analyze surface roughness and highlights that almost all models were developed for simplistic natural surfaces (i.e., regular arrays of roughness elements).The lack of a clear method for calculating z 0 as a function of surface roughness is due to the complexity of surfaces that exists in nature and the direction, spatial, and temporal dependence.
The most robust approach for estimating z 0 is from the anemometric method used to generate a logarithmic wind profile and solve for z 0 [32]. The anemometric method can be used for any surface with any arrangement of roughness elements but requires a meteorological tower of at least two vertically spaced wind, temperature, and humidity measurements that can be used to approximate the respective gradients. The measurements integrate over a footprint area rather than the single-point location of the sensors based on the distance from measurement source, elevation of sensor, meteorological conditions, turbulent boundary layer, and atmospheric stability. All of these factors can potentially create turbulent fluctuations affecting the downwind measurements of the wind profile [37,38]. The anemometric method is also very sensitive to the wind measurement heights; Munro [2] found that adding 0.1 m to any of the heights can alter z 0 by an order of magnitude.
In contrast, the geometric method uses algorithms relating z 0 to characteristics of surface roughness elements and thus does not require tower instrumentation but only a measure of the geometry of the surface [29].
Anemometric data are used to estimate z 0 from the logarithmic wind profile through an empirical relation that describes the vertical distribution of horizontal wind speeds within the lowest portion of the planetary boundary layer [39]. The wind speed (U z in m/s) at height z (in m) above a surface is given by: where U* is the friction velocity (m/s), k is the Von Kármán constant (~0.40), and ψ is a stability term, and L is the Monin-Obukhov stability parameter. This equation is only valid through the hypothesis of stationarity and horizontal homogeneity. Under neutral stability conditions, z/L tends towards zero, and ψ can be neglected. The most common geometric method for estimating z 0 is simply a function of the height of the elements: where z h is the mean height of roughness elements in meters, and f 0 is an empirical coefficient derived from observation [28]. The frontal area index, which combines mean height and breadth (all in meters), and density of the roughness elements, is defined as roughness area density given by: where Ly is the mean breadth of the roughness elements perpendicular to the wind direction, and ρel is the density or number (n) of roughness elements per unit area [40]. Lettau [29] developed a formula for z 0 based on the geometry of the surface for irregular arrays of reasonably homogenous elements: In the Lettau formula, the coefficient 0.5 represents an average drag coefficient for the roughness elements, which was determined experimentally. Other geometric methods have been developed, especially to consider more regularly-shaped and distributed roughness elements, such as buildings in an urban setting [41,42]. The Counihan equation provides a geometric estimate of z 0 as: where A f is the total area in square meters silhouetted by the roughness elements, and A d is the total area covered by roughness elements. A meteorological tower was erected at Colorado State University Agricultural Research, Development and Education Center (ARDEC) South (http://aes-ardec.agsci.colostate.edu/), (40.629680, −104.99699) on a flat field that had no obstructions at least 100 m in the prevailing wind direction. The fetch was 40 m wide with the tower placed in the middle, leaving 100 unobstructed, homogenous meters upwind. Ten anemometers and five temperature and relative humidity sensors were placed vertically at different heights on the tower. The accuracy of the air temperature and relative humidity sensors (METER VP-3) was variable across a range of ±0.25-0.50 • C and ±4%, respectively (see http://manuals.decagon.com/Manuals/14053_VP-3_Web.pdf for more information). The METER Davis Cup Anemometers have a wind direction accuracy of ±7 • and a speed accuracy within ±5% (see http://manuals.decagon.com/Manuals/). Data were collected from February 2014 through March 2015. In mid-March 2014, the flat field was plowed to create additional underlying roughness, specifically furrows and troughs were formed perpendicular to the dominant wind direction at an approximate spacing of two meters. The approximate amplitude of the troughs and furrows was 25 cm deep and 50 cm wide.
Meteorological data were recorded every five minutes based on the average of one-minute observations. Anemometric data were evaluated for 153 instances when wind speeds were faster than 4 m/s to ensure neutral stability [8] and when the log-linear fit had an r 2 greater than 0.95. The height of the instruments was calculated based on the depth of snow, which did not exceed 10 cm.
This study estimated z 0 values from anemometric measurements and used them as a reference to evaluate concurrent geometric methods. The z 0A values were calculated using Equation (1) from logarithmic anemometer wind profile data. Surface elevation was measured using a FARO Focus3D X 130 model TLS (https://www.faro.com/products/). This lidar tool generates a point cloud scan of a given area with an error of ±2 mm and a resolution of approximately 7.5 mm. The TLS was set up in 2-3 locations around the area of interest with 6 reference spheres to match the images using the FARO Scene Software. The data were generated into a point cloud and interpolated to a solid surface with 10 mm resolution with the kriging method using the Golden Software's Surfer 8 (https://www.goldensoftware.com/products/surfer). The gridded data were de-trended in the x-y plane to remove the bias in slope of the field or the angle of the lidar. Gaps in the scans tended to be small (<100 mm), and the kriging interpolation eliminated them. Individual roughness elements were identified and for each element the silhouette lot area and obstacle height were determined using a MATLAB code (https://www.mathworks.com/products/matlab.html). This was required to compute the Lettau formula (Equation (4)). The 1000-m 2 area around the tower was scanned on 12 occasions when the concurrent anemometric and geometric measurements were acquired. One concurrent measurement set was made with no snow cover for each of the unplowed and plowed scenarios; seven concurrent measurement sets were made with partial snow-covered area (SCA) and three with complete snow cover. SCA was determined from digital photos taken from the TLS unit.
Both the Counihan and Lettau methods were used to calculate z 0G (Equations (4) and (5), respectively). The Counihan method was appropriate for this study because the roughness elements (furrows) in this study site were semi-regular. During each concurrent anemometric and geometric measurement set, the percentage of the area covered in snow, or SCA, was estimated from photographs.

Results
The unplowed versus plowed field yielded different z 0A values ( Figure 1). On average, the plowed field was almost 20 times as rough as the unplowed field, yet the coefficient of variation (COV) was essentially the same (0.67 and 0.62, respectively) ( Figure 1). The smallest z 0A values for the plowed field were of the same magnitude as some of the largest z 0A values for the unplowed field, in the range of 1 to 3 × 10 −3 m.
The Counihan method estimated z 0G values that were 1.39 times larger and had greater variation than the estimated z 0A values ( Figure 2). We used the Nash-Sutcliffe coefficient of efficiency (NSCE), which is a performance statistic based on a comparison of the data fit to the 1:1 line, to evaluate how estimates of z 0G compared with z 0A [43]. The NSCE of the Counihan z 0G was −1.18, and the Lettau z 0G was 0.14, indicating the Lettau method compared more favorably with the z 0A . A linear regression between both z 0G estimates (Counihan z 0G and Lettau z 0G ) and z 0A was fit through the data origin to evaluate if the bias between the two methods could be removed through simple linear scaling (Figure 2). When the Counihan z 0G values were scaled by 0.721 (1/1.39), the NSCE value only increased to 0.07. However, the NSCE increased to 0.88 when the Lettau z 0G values were scaled by 2.34 (1/0.428).

193
The estimated z0 values were found to vary as a function of the amount of SCA present ( Figure   194 3). As SCA increases, z0 decreases, with variability based on the calculation method (Figure 3a)  The estimated z 0 values were found to vary as a function of the amount of SCA present ( Figure 3). As SCA increases, z 0 decreases, with variability based on the calculation method (Figure 3a). A linear regression between SCA and each of the z 0 estimates showed r 2 values that were 0.01, 0.7, and 0.88 for the Counihan, anemometric, and Lettau methods of z 0 calculation, respectively. There were noticeable differences in z 0 depending whether SCA was increasing because snow was accumulating versus when SCA was decreasing because the snow was melting. For periods of snow accumulation, removing snow measurements that were not immediately following a snow event (the yellow boxes in Figure 3b that represent non-accumulation values) improved the linear relation between accumulating SCA and z 0 (R 2 = 0.94).  221 Figure 3. Linear relation between z 0 and snow-covered area (SCA as a %) for (a) all datasets (scaled Lettau and Counihan geometry-based and anemometric-based) with anemometric-based z 0 for the pre-plowed (orange circles) and plowed fields (green triangles) highlighted, and (b) the scaled Lettau geometry-based z 0 using a factor of 2.34 (see Figure 2). Lettau geometry-based z 0 measurements with non-accumulation snow measurements were removed. Lines are based on the best-fit linear regression of the data. Snow had been on the ground for numerous days prior to the two concurrent measurements (yellow squares) taken on 22 March 2014 (SCA = 100%) and 13 April 2014 (SCA = 70%). The snowfall was fresh for all other measurements.

Discussion
Geometrically estimated z 0 , although easier to measure, produced different values when compared to the anemometric derived values. The Counihan method overestimated by a factor of 0.721, whereas the Lettau method underestimated anemometric z 0 ( Figure 2) by a factor of 2.34. The Lettau method (Equation (4)) has a constant of 0.5 based on the average drag coefficient of the roughness characteristic of the silhouetted area of the average obstacle. By dividing the Lettau based z 0 values by the 0.5 and thus eliminating the drag coefficient from the equation, we get a new NSCE of 0.856, with scatter in the data much closer to the 1:1 line (Figure 2). The removal of the drag coefficient suggests that the geometric data generated from the lidar point cloud appears to account for spatial and temporal variability in the roughness of a snow surface.
Lidar-based snow data are becoming more readily available [19,26]. The accuracy of the scans from about 1 mm for terrestrial lidar to 10 cm for airborne lidar can account for fine-scale changes of the snowpack [26], which enables the computation of z 0 at any scale. Although anemometric data can yield reliable estimates of z 0 , meteorological towers are expensive to set up and operate. In addition, data from a single tower does not consider spatial variability as well as the geometric method [44]. Comparing the two methods does not consider the scale of the study area; the geometric measurement is taken over the entire area near the meteorological tower whereas the anemometric measurement is only influenced by the fetch area upwind of the sensors [29].
The roughness of the snowpack can vary substantially both spatially and temporally creating many implications [13,14,45]. Roughness variations can be caused by heterogeneities in land cover, vegetation, and meteorological conditions [46]; non-uniform distribution of snow cover during accumulation and melt [45,46]; snow-canopy interactions [47]; and snow redistribution by wind [48]. This was apparent in differences between the estimated z 0 for the plowed versus unplowed field ( Figure 1). Land cover varies throughout regions particularly those with a shallow snow environment, and this creates variations that depend on the underlying topography [13,14,46]. Thus, there are many different values of z 0 in the literature [7] that are broader than our observed mean range of 0.2 to 10 × 10 −3 m (Figure 1). For example, Miles et al. [31] found the z 0 of a hummocky glacier (a particularly rough underlying surface) to range between 5 to 500 × 10 −3 m, whereas Brock et al. [7] reported z 0 values for fresh snow and older snow of 0.2 × 10 −3 and 3.56 × 10 −3 m, respectively. Our results show that change in roughness between a plowed and unplowed field yielded a 20-fold difference in z 0 . The results of this study can be applied to areas of similar climate and land cover, which included flat, bare soil, and bare soil with small furrows (<1 m); and therefore, the results of this study may not scale appropriately to different land cover types. Further studies of a shallow snowpack in sagebrush steppe [49], farmland, or non-densely forested environments may be able to replicate our study results and scale from 1000 m 2 to a larger area. The z 0 values observed here had a notable change between flat soil and small furrows, so the changes in z 0 values in different environments with even minimal vegetation will have much larger effects on the z 0 values.
The inverse relation of SCA and z 0 ( Figure 3) [50] is affected by the underlying terrain and size of the roughness features. As the snow accumulation increases, the roughness elements become buried, and the topography appears to be smooth [50,51]. This relation indicates that as snow accumulates over topographic features the snow will begin to level out at a z 0 height dependent scale. A hysteresis can be noted, and it has been found that a single snowfall event on a hummocky glacier can alter the micro-topography by up to 75% due to the shallow snowpack over the small scaled features [45,52]. The CLM4 uses a z 0 value of 2.4 × 10 −3 m, a value that falls near the mean of the unplowed field, which is applicable for deep, flat snowpack surface with minimal influences from underlying terrain. However, this is not typical for shallower snowpacks or in complex terrains.
Relations between z 0 and SCA (Figure 3) can be used to improve snow-energy balance modeling by estimating the percentage of SCA via remote sensing and applying z 0 only to the portion of area it accurately describes [46,53]. Currently, most models use 100% SCA even though many areas will remain snow free due to complex terrain and can drastically change during periods of melt and accumulation [13,53]. Aerodynamic roughness length is incorporated into many climate and energy models, which require sub-grid snow distribution [54] and are still inadequate at representing SCA [46,48]. A dynamic z 0 based on SCA and land cover type can improve these on a sub-grid scale. Another complication with these models is the lack of accountability for snowpack variability throughout accumulation and melt [48,53,54].
Resolution is an important factor to consider when discussing both SCA and z 0G . The higher the resolution of the measurements (lidar, satellite, etc.), the higher the z 0-G accuracy. However, lidar datasets are often large, especially those acquired with TLS, making them difficult and time consuming to process. Lower resolution data from remote sensing or airborne lidar systems (ALS) can cause problems when scaling [53]. Quincey et al. [52] found that z 0G is typically underestimated with a small area and coarse resolution and overestimated with a large area and fine resolution when compared to anemometric data. Nonetheless, even with lower resolution, applying dynamic z 0 values may greatly improve models. Scaling can be an effective way to incorporate both an anemometric and geometric z 0 value. Based on a specific land cover type, a scaling factor can be applied to areas with the same land cover. This can help to improve modeled z 0 accuracy, once preliminary z 0 values have been established.

Conclusions
Aerodynamic roughness length within our study has shown variation spatially and temporally for a shallow accumulating snow environment. This was apparent in our results that showed differences in z 0 mean values of about 14 × 10 −3 m between the plowed and unplowed field. Thus, single-point measurements of anemometric data may not account for z 0 over a range of spatial and temporal scales. Geometrically calculated z 0 using the Lettau method has shown to be an effective and more robust form of z 0 estimation compared to the anemometric method and also producing similar, estimated values. The anemometric, single-point measurements will also not account for the snow-covered area, which changes based on its inverse relation with z 0 . However, SCA can be observed and estimated from satellite imagery or airborne lidar systems to create a more accurate estimation of z 0 .