**4. 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 ge<sup>t</sup> 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−<sup>3</sup> 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−<sup>3</sup> m, whereas Brock et al. [7] reported *z*0 values for fresh snow and older snow of 0.2 × 10−<sup>3</sup> and 3.56 × 10−<sup>3</sup> 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<sup>2</sup> 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−<sup>3</sup> 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*0*G*. 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*0*G* 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.
