**1. 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/) applies a single *z*0 value of 2.4 × 10−<sup>3</sup> m to all snow-covered surfaces. However, *z*0 varies both spatially [7] and temporally [8], which may result in variable estimates of turbulent heat fluxes not captured by most models [9]. Wind velocity profile measurements are often used to calculate *z*0 estimates [6], but there are a limited number of sites that measure the wind profile over a snowpack surface, making the spatio-temporal representation of *z*0 challenging.

Millimeter-scale variations in snow-surface roughness features have been estimated from a black plate pushed partially into the snow [10–12], using two-dimensional photography, digital processing, and automated post-processing software [13–16]. Snow surface elevation data are now available over large areas at the resolution (±80 mm) of airborne light detection and ranging (lidar) [17–19], terrestrial laser scanner (TLS) (resolution of ±5 mm) [20–26], and photogrammetry [25]. Although most lidar and photogrammetry efforts have only focused on snow depth [26], only a few datasets have been used to evaluate snow surface roughness at the meter-scale or sub-meter scale [27,28]. However, few of these datasets have been applied to interpolate *z*0 and create a digital elevation model of the snowpack surface for evaluating surface roughness [27]. Aerodynamic roughness length (*z*0) has been estimated from the geometry of the snow surface [2,7,29–31]. However, this method is time consuming and typically only applicable over smaller scales [13]. Also, Fassnacht et al. [27] have identified potential errors with the different methods of computing *z*0 from the geometry of the surface that result in values varying over 1–3 orders of magnitude and have suggested these methods need to be evaluated for varying scales, resolutions, and environments.

This study used TLS-derived surface geometry and vertical wind profile measurements to compare concurrent *z*0 estimates for changing snow surface features of shallow snowpacks. Here, we asked the following questions: (1) How does the aerodynamic roughness length (*z*0) vary spatially and temporally for a shallow snow environment? (2) How does *z*0 estimated from geometric measurements (*<sup>z</sup>*0*G*) compare to *z*0 estimated from anemometric measurements (*<sup>z</sup>*0*A*), and (3) How does *z*0 vary with snow-covered area based on the underlying terrain?

#### **2. 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 arrangemen<sup>t</sup> 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 (*Uz* in m/s) at height *z* (in m) above a surface is given by:

$$\mathcal{U}z = \frac{\mathcal{U}^\*}{k} \ln\left[\frac{z}{z\_0} + \psi\left(\mathbf{z}, \frac{z}{L}, L\right)\right] \tag{1}$$

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:

$$z\_0 = f\_0 z\_h \tag{2}$$

where *zh* 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:

$$(\lambda F) = \text{ l.y } z\_h \text{ } \rho el \tag{3}$$

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:

$$z\_0 = 0.5 \ z\_h \,\lambda \,\text{F} \tag{4}$$

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:

$$z\_0 = \ z\_h (1.8 \frac{A\_f}{A\_d} - 0.08) \tag{5}$$

where *Af* is the total area in square meters silhouetted by the roughness elements, and *Ad* 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 r2 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*0*A* 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<sup>2</sup> 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*0*G* (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.
