Predicting the roughness length of turbulent flows over 1 landscapes with multi-scale microtopography 2

8 The fully rough form of the law of the wall is commonly used to quantify velocity 9 profiles and associated bed shear stresses in fluvial, aeolian, and coastal environments. A key 10 parameter in this law is the roughness length, z 0 . Here we propose a predictive formula for z 0 that 11 uses the amplitude and slope of each wavelength of microtopography within a discrete-Fourier-12 transform-based approach. Computational fluid dynamics (CFD) modeling is used to quantify 13 the effective z 0 value of sinusoidal microtopography as a function of the amplitude and slope. 14 The effective z 0 value of landscapes with multi-scale roughness is then given by the sum of 15 contributions from each Fourier mode of the microtopography. Predictions of the equation are 16 tested against z 0 values measured in ~10 5 wind velocity profiles from southwestern U.S. playa 17 surfaces. Our equation is capable of predicting z 0 values to 50% accuracy, on average, across a 18 four order-of-magnitude range. We also use our results to provide a simpler alternative formula 19 that, while somewhat less accurate than the one obtained from a full multi-scale analysis, has an 20 advantage of being simpler and easier to apply. 21


Problem statement
The velocity profiles of turbulent boundary-layer flows are often quantified using the fully rough form of the law of the wall: where u(z) is the wind velocity (averaged over some time interval) at a height z above the bed, u * is the shear velocity, κ is the von Kármán constant (0.41), and z 0 is an effective roughness length that includes the effects of grain-scale roughness and microtopography (e.g., Bauer et al., 1992;Dong et al., 2001).Velocity profiles measured in the field are commonly fit to Eq. ( 1) to estimate u * and/or τ b for input into empirical sediment transport models (often after a decomposition of the bed shear stress into skin and form drag components) (e.g., Gomez and Church, 1989;Nakato, 1990).Fits of wind-velocity profiles to Eq. ( 1) also provide measurements of z 0 .Given a value for z 0 , a time series of u * and/or τ b can be calculated from Eq. (1) using measurements of velocity from just a single height above the ground.This approach is widely used because flow velocity data are often limited to a single height.Equation (1) only applies to z ≥ z 0 , and may be further limited in its accuracy within the roughness sublayer, i.e., the range of heights above the ground comparable to the height of the largest roughness elements.The roughness sublayer is the layer where the mean velocity profile deviates from the law of the wall as the flow interacts with individual roughness elements.This layer is typically considered to extend from the ground surface to a height of approximately twice the height of the tallest roughness elements.
Values of z 0 depend on microtopography/land cover (quantifying this dependence in unvegetated landscapes is a key goal of this paper) and are typically in the range of 10 −2 -10 1 mm for wind flow over arid regions (Prigent et al., 2005).Most existing methods for estimating z 0 using metrics of surface roughness or microtopography rely on the concept of a dominant roughness element, the size and density of Published by Copernicus Publications on behalf of the European Geosciences Union.which the user must specify a priori (e.g., Lettau, 1969;Arya, 1975;Smith and McLean, 1977;Jacobs, 1989;Taylor et al., 1989;Raupach, 1992Raupach, , 1994;;Kean and Smith, 2006a).Procedures are available for estimating z 0 in landscapes with multi-scale roughness, but they often rely on idealizations such as treating the microtopography as a sequence of Gaussian bumps (e.g., Kean and Smith, 2006b).Nearly all natural landscapes have microtopographic variability over a wide range of spatial scales.Identifying the dominant scale objectively and uniquely can be difficult.For example, the top plot in Fig. 1 shows a hypothetical case of a landscape composed of two superposed sine waves.The effective roughness length of a landscape is related to the presence/absence or extent of flow separation, and flow separation is primarily controlled by the derivatives of topography (slope and curvature) rather than the amplitude of the bedforms/roughness elements (Simpson, 1989;Lamballais et al., 2010).As such, roughness elements of smaller amplitude but steeper slopes may exert greater control on z 0 values compared with roughness elements that are larger in amplitude but gentler in slope.
Given a landscape with multi-scale roughness in which each scale has distinct amplitudes and slopes, it can be difficult to identify the dominant scale or scales of roughness for the purposes of estimating z 0 .
Figure 1 illustrates two examples of microtopography from playa surfaces in the southwestern US.The middle plot shows a transect through the Devil's Golf Course in Death Valley, California, and the bottom plot shows a transect through a relatively smooth section of Lordsburg Playa, New Mexico.These plots are presented using different vertical scales because the amplitude of the microtopography at the Death Valley site is approximately 100 times greater than that of the Lordsburg Playa site.Both landscapes have no vegetation cover, no loose sand available for transport, and are flat or locally planar at scales larger than ∼ 1 m.As such, they are among the simplest possible natural landscapes in terms of their roughness characteristics.Nevertheless, as Fig. 1 demonstrates, they are characterized by significant roughness over all spatial scales from the resolution of the data (1 cm) up to spatial scales of ∼ 1 m.To our knowledge, there is no procedure for predicting z 0 in a way that honors the multi-scale nature of microtopography in real cases such as these.To meet this need, we have developed and tested a discrete-Fourier-transform-based approach to quantifying the effects of microtopographic variations on z 0 values.The method simultaneously provides an objective measure of the spatial scales of microtopography/roughness that most strongly control z 0 .
In a recent paper similar in spirit to this one, Nield et al. (2014) quantified the z 0 values of wind-velocity profiles over playas as a function of various microtopographic metrics.Nield et al. (2014) proposed an empirical, powerlaw relationship between z 0 and the root-mean-squared variations of microtopography, H RMSE : where the coefficient c is equal to ln(−1.43) or 0.239 m −0.34 .Equation ( 2) is one example of several predictive formulae that Nield et al. (2014) proposed for different surface types.Equation ( 2) applies to surfaces with large roughness elements or that exhibit mixed homogenous patches of large and small roughness elements.Nield et al. (2014) concluded that "the spacing of morphological elements is far less powerful in explaining variations in z 0 than metrics based on surface roughness height".In this paper we build upon the results of Nield et al. (2014) to show that z 0 can be most accurately predicted using a combination of the amplitudes and slopes of microtopographic variations.
The presence of multi-scale roughness in nearly all landscapes complicates attempts to quantify effective z 0 values for input into regional and global atmospheric and Earthsystem models.In such models, topographic variations are resolved at scales larger than a single grid cell (10-100 km at present, but steadily decreasing through time as computational power increases) but the aerodynamic effects of topographic variations on wind-velocity profiles at smaller scales are not resolved in these models and must be represented by an effective z 0 value (sometimes in combination with an additional parameter, the displacement height, which shifts the location of maximum shear stress to a location close to the top of the roughness sublayer; Jackson, 1981).Topographic variations at spatial scales below 10-100 km are typically on the order of tens to hundreds of meters.Currently available maps of z 0 values do not incorporate the aerodynamic effects of topography at such scales.For example, Prigent et al. (2005) developed a global map of z 0 in deserts by correlating radar-derived measurements of decimeter-scale roughness with z 0 values inferred from wind-velocity profiles.This approach assumes that the dominant roughness elements that control the effective z 0 value over scales of 10-100 km occur at the decimeter scale captured by radar.It is possible that, in some landscapes, the roughness that controls z 0 occurs at scales that are larger or smaller than those measured by radar.Therefore, a procedure is needed that predicts z 0 values using data for topographic variations over a wide range of scales, including but not limited to decimeter scales.This study aims to fill that gap.

Study sites
We collected wind-velocity profiles and high-resolution topographic data using terrestrial laser scanning (TLS) from 10 playa sites in the southwestern US (Fig. 2) during the spring of 2015.These sites were selected based on the range of microtopographic roughness they exhibit (Table 1).Roughness can be quantified in multiple ways, but H RMSE , the root-mean-squared deviation of elevation values measured at a sampling interval of 0.01 m, provides one appropriate metric (Nield et al., 2014).The 10 sites range in H RMSE from approximately 0.55 to 36 mm (see Sect. 2.1).In addition to the H RMSE we computed S av , the average slope computed at 0.01 m scale, for each site.Values of S av range from 0.01 to 0.159 m/m (Table 1).
Each study site was an area of at least 30 m×30 m with relatively uniform roughness, as judged visually and by analysis of the TLS data.The minimum fetch required for an equilibrium boundary-layer flow is typically assumed to be 1000 times the height of the dominant roughness elements (Counehan, 1971).Based on this criterion, 30 m was adequate fetch for 7 of the 10 sites, i.e., all except for the three Death Valley sites, where roughness elements were up to 300 mm; hence, the area of homogeneous roughness was verified to a distance of only ∼ 100 times the height of the dominant roughness elements.However, the required fetch must also depend on the maximum height above the ground where velocities are measured to compute a z 0 value locally, since any roughness transition triggers an internal boundary layer that grows indefinitely in height with increasing distance downwind of the transition.Using the Elliot (1958) formula for the height of The playa surfaces at our study sites were predominantly crusted and ranged from flat, recently formed crust to wellformed polygons with deflated and broken crust ridges.All of the sites were completely devoid of vegetation.Sand blows episodically across some portions of the playas we studied, but we chose study areas in which we observed no sediment transport during fast winds.We considered only landscapes without vegetation and loose, erodible sand because such cases must be understood first before the additional complications of flexible roughness elements and saltation-induced roughness can be tackled.That said, we anticipate that concepts from this paper may be relevant to quantifying z 0 over vegetated landscapes also.
Our goal is to understand the controls on boundary-layer flows over rough terrain generally, not playa surfaces specifically or exclusively.As such, we use playa surfaces as "model" landscapes.Playas are useful for this purpose because they are macroscopically flat but exhibit a wide range of microtopographic roughness at small scales.The relative flatness of playas at scales larger than ∼ 1 m makes it possible to characterize their boundary-layer flows using relatively short anemometer towers.Of course, playas are also of special interest to aeolian geomorphologists because they can be major dust sources when sand from playa margins is transported across the playa surface, disturbing crusted surfaces and liberating large volumes of silt-and clay-rich sediments.
The questions addressed in this paper could, in principle, be addressed using wind tunnel experiments.Wind tunnels certainly have the advantage of user control over wind velocities.However, Sherman and Farrell (2008) documented that z 0 values in wind tunnels are, on average, approximately an order of magnitude lower than those measured in the field for otherwise similar conditions (e.g., grain size).One interpretation of the Sherman and Farrell (2008) results is that the confined nature of wind tunnel flows and/or their limited fetch can limit the development of boundary layers in equilibrium with bed roughness.For this reason, we focused on measuring wind flow over natural surfaces with homogeneous roughness characteristics over distances of at least 30 m surrounding our measurement locations.

Terrestrial laser scanning and analyses of playa surface microtopography
A Leica C10 terrestrial laser scanner was used to acquire point clouds of the central 10 m × 10 m ground surface upwind of the anemometers at each of the 10 study sites.The areas surrounding each 10 m × 10 m area were also surveyed to check for approximate homogeneity in the roughness metrics out to areas of 30 m × 30 m, but the central 10 m × 10 m areas were the focus of the subsequent data analysis.Each area was scanned from four stations surrounding the 10 m × 10 m area and merged into a single point cloud using a Leica disk target system.Registration errors were a maximum of 2 mm in all cases.The Leica C10 has an inherent surface-model accuracy of 2 mm, but this value decreases as the number of overlapping scans increases (Hodge, 2010), resulting in a value of approximately 1 mm in the case of four overlapping scans.
The scanner was mounted on a 3.5 m tripod to maximize the angle of incidence (low angles of incidence elongate the "shadows" or occlusions behind microtopographic highs; Brown and Hugenholtz, 2013).All of the returns within each 1 cm 2 domain were averaged to create a digital elevation model (DEM) with point spacing of 0.01 m.Voids were filled using natural-neighbor interpolation.Voids requiring interpolation were limited to < 1 % of the area at the smoothest five sites (Lordsburg and Willcox playas), between 1 and 3 % at the two Soda Lake sites, and between 10 and 20 % at the three Death Valley sites.
In addition to the calculation of basic topographic metrics such as H RMSE and S av (the latter being the average slope computed at 0.01 m scales) (Table 1), we also computed the average amplitude spectrum of all 1-D topographic transects at each study site.The amplitude spectrum is equal to 2 times the absolute value of the complex discrete Fourier transform (DFT).The average amplitude spectrum refers to the fact that the one thousand amplitude spectra of each 1-D transect computed along the east-west direction were averaged to obtain a single average spectrum for each study site.We used the DFT implemented in the IDL programming language.The DFT coefficients were also used as input to a filter that uses the amplitude and slope of each Fourier mode to compute its contribution to the z 0 value.We created "mirror" images of each transect before application of the DFT.This approach has been shown to work as well or better than windowing for minimizing truncation error (i.e., incomplete sampling) in data sets characterized by the broadband/multiscale variability characteristic of many environmental data series (Smigelski, 2013).

Measurement and analyses of wind profiles
Wind speeds were measured at 1 s intervals and at seven heights above the surface (0.01, 0.035, 0.076, 0.16, 0.52, 1.22, and 2.80 m) using four Inspeed Vortex rotating cup anemometers and four AccuSense hot-wire anemometers (F900 series) (the latter calibrated to work over the 0.15-10 m s −1 range of wind velocities) (Fig. 3).The hot-wire sensors were secured to an L-shaped steel frame and placed above the surface such that the small opening in the sensor head was oriented as perpendicular to the wind direction as possible (Fig. 3).The 10 m s −1 range of the hot-wire sensors was not a limiting factor because all of the hot-wire sensors were located close to the ground, i.e., within 0.16 m from the surface, where velocities were lower than 10 m s −1 during our deployments.We collected data at each of the 10 sites for 10 to 30 h spanning multiple days, times of day, and a wide range of wind velocities.
The lowest cup and the highest hot-wire anemometers were positioned at the same height (0.16 m) above the surface to standardize measurements between the two types of wind sensors.When positioned at the same height, the hotwire sensors measured wind speeds (based on the factory calibration) that were approximately 10 % lower than the values obtained from the cup anemometers.We used the ratio of the wind velocities measured by the bottom cup anemometer to the wind velocities measured by the top hot-wire sensor to standardize the hot-wire measurements to the cup anemometer measurements in real time.This scaling-factor approach also serves a second purpose, which is to minimize the effects of wind-direction variability on the velocities measured by the hot-wire sensors.The cup sensors measure wind speeds effectively from nearly any direction, but the hot-wire sensors are required to be oriented within 20 • perpendicular to the wind for greatest accuracy.The hot-wires were manually repositioned following large and sustained changes in wind direction, but short-duration changes may have resulted in oblique incidence angles with a bias towards lower velocities.Continually rescaling the velocities measured by the highest hot-wire sensor to the lowest cup sensor mitigated this potential problem.
Scaled values from the bottom three (0.01, 0.035, and 0.076 m) hot-wire sensors were combined with the four cup anemometers to calculate shear velocities, u * , and aerodynamic roughness lengths, z 0 , based on the average velocities measured in each 12 s interval via least-squares fitting of the wind velocities to the natural logarithm of the distance www.earth-surf-dynam.net/4/391/2016/Earth Surf.Dynam., 4, 391-405, 2016 above the ground.To extract a z 0 value from the velocity profile data, we followed the procedure of Bergeron and Abrahams (1992), who emphasized the need to regress u on ln z rather than ln z on u.The shear velocity is equal to the slope of the regression of u on ln z multiplied by κ (Eq.6 of Bergeron and Abrahams, 1992) and the roughness length is equal to the exponential of the following: minus the intercept divided by the slope (Eq.7 of Bergeron and Abrahams, 1992).
The shear velocity is equal to the slope of u vs. ln z divided by κ.The roughness length is equal to the exponential of the following: minus the intercept divided by the slope.The 12 s interval was chosen based on the results of Namikas et al. (2003), who found that time intervals greater than 10 s resulted in the most accurate results, while those obtained from smaller averaging intervals were less reliable.Values of z 0 can be influenced by deviations from neutral stability.A common way to address this issue is to remove profiles from the analysis in which the velocity at a given height is below some threshold value (e.g., Nield et al., 2014).In this study we repeated our analysis using only those profiles with a wind velocity of at least 3 m s −1 at a height of 0.16 m.The mean and standard deviations of z 0 were nearly identical to those obtained using all of the data, likely reflecting the fact that we targeted time periods of fast winds for measurement.
During the data collection, the hot-wire sensors were moved to approximately 25-50 random locations within each site.We moved the hot-wire sensors to numerous locations within each site because wind velocities measured close to the ground are sensitive to the microtopography of the specific spot above which they are measured -i.e., the z 0 value measured on the stoss side of a microtopographic high tends to be smaller than the z 0 value measured on the lee due to the convergence/divergence of flow lines.Since our goal was to characterize the average or representative z 0 value over each surface, it is appropriate to move the hot-wire sensors around the surface to ensure that the results are not specific to one location but instead represent a statistical "sample" of the flow above the surface at multiple locations.This approach is also consistent with how the computational fluid dynamics (CFD) model output was analyzed (see Sect. 2.3).
Velocity profiles can deviate from Eq. ( 1) close to the ground over rough terrain.As such, it is important to identify which sensors, if any, are located within the roughness sublayer prior to computing u * and z 0 values by fitting windvelocity data to Eq. (1).To do this, we plotted the average of all wind-velocity measurements at each site as a function of ln z.The results (described in Sect.3.2) show that the lowest two (hot-wire) sensors (located 0.10 and 0.035 m above the ground) at the three Death Valley sites and the rough Soda Lake site deviated noticeably from Eq. (1).The fact that these sensors were within the roughness sublayer is consistent with the fact that the height of the largest roughness elements at these sites is greater than or comparable to 0.035 m (the height of the second lowest sensor).Data from the lowest sensor at the next four smoothest sites (i.e., smooth Soda Lake, the two Willcox Playa sites, and the rough Lordsburg Playa site) also deviate noticeably from Eq. (1).Data from these sensors were not used in the calculation of u * and z 0 at those sites.In addition, we verified in all cases that the removal of these sensors deemed to be within the roughness sublayer improved the mean correlation coefficients, R 2 , at each site.Only profiles with R 2 values greater than 0.95 were retained.

Computational fluid dynamics
CFD modeling was used to quantify the effects of the amplitude and slope of sinusoidal microtopography on z 0 .We used the 2013 version of the PHOENICS CFD model (Ludwig, 2011) to estimate the time-averaged wind velocities associated with neutrally stratified turbulent flow over sinusoidal topography with a range of amplitudes and slopes.PHOEN-ICS uses a finite-volume scheme to solve simultaneously for the time-averaged pressure and flow velocity.PHOENICS solves the flow equations using the iterative SIMPLEST algorithm of Spalding (1980), which is a variant of the SIM-PLE algorithm of Patankar and Spalding (1972).The solution was considered converged when the state variables changed by less than 0.001 % from one iteration to the next.We used the renormalization group variant of the k-ε closure scheme first proposed by Yakhot and Orszag (1986) and later updated by Yakhot et al. (1992), which is widely used for sheared/separated boundary-layer flows.
Inputs to our model runs include a topographic profile (in these cases, a sinusoid of a prescribed amplitude and maximum slope), a grain-scale roughness length, z 0g (set to be 0.002 mm for all runs), and a prescribed horizontal velocity at a reference height far from the bed (u r = 10 m s −1 at z r = 10 m was used for all of the model runs presented).The value of z 0g was chosen based on the measured value of z 0 at the two flattest sites (Lordsburg smooth and intermediate), both of which yield z 0 = 0.002 mm as described in Sect.3.2.This value is also consistent with the grain-scale roughness expected at a site with a median grain size of very fine sand if the Bagnold (1938) relation z 0g = d 50 /30 is used.The ground surface is prescribed to be a fully rough boundary, i.e., one that results in a law of the wall velocity profile characterized by a roughness length equal to z 0g (0.002 mm) and a shear velocity equal to κu r / ln(z r /z 0g ) (0.26 m s −1 ) in the absence of topography.In the CFD model the ground surface is treated using a wall-function approach, i.e., the velocity profile within the first cell is assumed to be logarithmic with a microscale roughness length equal to z 0g if the flow is turbulent; otherwise, a laminar profile is used based on the viscosity of air.At the upwind boundary of the model domain an "inlet" law of the wall velocity profile is prescribed with a roughness length equal to z 0g .At the downwind boundary (i.e., the "outlet") a fixed-pressure boundary condition is used., (c, d), (e, f), and (g, h).
The computational grids we used consisted of 2-D terrainfollowing coordinate systems.Thirty logarithmically spaced grid points were used in the vertical direction, ranging from 0.1 mm to 10 m above the bed.We used 2000 grid points in the horizontal direction.The absolute size of the horizontal domain varied depending on the slope of the bedforms.That is, the topographic profile was identical for all of the runs (except for the fact that an amplitude of 0.05 m used for half of the runs and an amplitude of 0.1 m was used for the other half).Steeper slopes were obtained by decreasing the horizontal grid spacing or "compressing" the input topographic profile horizontally.The minimum length/fetch of the model domain was 30 m.Our analysis of the wind profiles output by the model was restricted to the last 20 % of the model domain, i.e., the portion farthest downwind.This was necessary because the upwind boundary of the model is a roughness transition triggered by the interaction of the input velocity profile (characterized by roughness length z 0g ) with the microtopography.This roughness transition generates an internal boundary layer that grows with distance from the upwind end of the domain.Within the internal boundary layer, the velocity profile is characterized by an effective roughness length z 0 set by the amplitude and slope of the bedforms.To properly compute the value of z 0 based on velocity profiles from the top of the roughness sublayer to a height of 3 m, it is necessary to restrict the analysis of the wind profiles to the downwind end of a model domain that is at least 30 m in length as described in Sect.1.2.
Model runs were performed using two different amplitudes (0.05 and 0.1 m) and a range of maximum slopes from 0.001 to 2.0.Each of the 400 vertical velocity profiles of the last 20 % of the model domain were fit to Eq. ( 1) from the top of the roughness sublayer (assumed to be equal to twice the maximum height of the bedforms) to a height of 3 m (to match the maximum height measured in the field).The 400 z 0 values were then averaged to obtain an effective z 0 value for each value of the sinusoidal amplitude and slope.Values of z 0 were fit to the expression where a is the amplitude (in units of m) of the sinusoid, S is the maximum slope (the slope at the point of inflection of the sinusoid in units of m/m), and c 1 , c 2 , and c 3 are unitless coefficients.
2.4 Fourier analysis of topography and a multi-scale approach to quantifying z 0 The results of the CFD modeling (described in Sect.3.3) suggest that the slope and amplitude of microtopographic variations control z 0 values via the sigmoidal relation of Eq. ( 3).This result provides a basis for quantifying the multi-scale controls on z 0 within a discrete-Fourier-transform-based approach that treats each Fourier mode as a sinusoid, uses Eq. ( 3) to quantify the effective roughness associated with each sinusoid, and then sums the contributions of each sinusoid to determine the total effective z 0 value, fully taking into account microtopographic variations across a wide range of scales.Within the implementation of the DFT in the IDL programming language, the amplitude of each Fourier mode is equal to 2 times the amplitude of the complex Fourier coefficient, i.e., a n = 2|f n |, and the maximum slope is given by S n = 2π ka, where k is the natural wavenumber.As such, the generalization of Eq. ( 3) to multi-scale topography as quantified using the DFT is where c 4 is a unitless coefficient analogous to c 1 , but with a potentially different value, and k is the natural wavenumber defined as the inverse of the wavelength.We verified that Eq. ( 4) returns the same value of z 0 as predicted by Eq. ( 3) for the case of a sinusoidal bed if c 4 = c 1 .We also verified that the z 0 values predicted by Eq. ( 4) were independent of the total number of data points and the sampling interval of the input data (provided that the dominant scales of roughness were represented and resolved).The best-fit value of c 4 was obtained by a trial and error minimization of the leastsquared difference between the predictions of Eq. ( 4) and the mean z 0 values measured at the 10 sites.An alternative approach to Eq. ( 4) that is easier to apply and does not rely on the Fourier transform is where c 5 and c 6 are unitless coefficients.

TLS surveying
Figure 4 presents color maps of the topography of the roughest and smoothest sites at each playa.Table 1 presents summary statistics for the 10 sites, including the topographic metrics H RMSE and S av .
Figure 5 plots the average amplitude spectrum of all 1-D topographic transects for each site.These spectra demonstrate that significant topographic variability exists at all spatial scales of measurement, i.e., from 0.02 to 10 m (note that two samples are required for Fourier analysis; hence, the smallest wavelength captured in our analysis is 2 x or 0.02 m).A surface with a single scale of roughness, such as wind ripples, would have power concentrated at a narrow range of wavelengths, unlike the "broadband" spectra of Fig. 5. Also, note that the different shapes of the spectra reflect the different spatial scales that dominate topographic variability at each site.At Willcox Playa, for example, the largest roughness elements occur at horizontal spatial scales ∼ 1-3 m (Fig. 4d and e).As a result, the power spectra for the Willcox sites exhibit a "bend" at wavelengths of approximately 1-3 m, indicating that the amplitude of the microtopography drops off substantially at wavenumbers larger than 0.3-1, i.e., wavelengths smaller than 1-3 m.A similar bend occurs in the Lordsburg rough site but at a higher wavenumber corresponding to a wavelength of ∼ 0.03-0.1 m.The color map of the Lordsburg rough site is consistent with this -i.e., it shows a "dimpled" surface with large roughness elements ∼ 0.03-0.1 m in size.

Measurement and analyses of wind profiles
Figure 6 plots the relationship between the average wind velocity (normalized to the value measured at 2.8 m above the ground) and the natural logarithm of height above the ground for all sites.The data have been normalized to emphasize how z 0 and deviations from Eq. ( 1) vary among the sites (neither of which depend on absolute velocity values).Note that the three Death Valley sites have been shifted to the left along the x axis by 0.1 m s −1 to help differentiate the plots.
The law of the wall predicts a constant slope when u is plotted vs. ln z.When the velocities are normalized as in Fig. 6, a steeper slope corresponds to a smaller z 0 value.The slopes of the lines in Fig. 6 systematically decrease (hence mean z 0 values increase) from the smoothest playa (Lordsburg) to the roughest (Death Valley).Within each playa, the slopes also systematically decrease from relatively smooth sites to rough sites (Table 1).The plots in Fig. 6 suggest that the lowest two sensors (located 0.01 and 0.035 m above the ground) at the Death Valley sites and the rough Soda Lake site reside within the roughness sublayer and hence should not be used to obtain z 0 values via least-squares fitting of data to Eq. ( 1).The same is true for the lowest sensor at the four smoother sites (all but the two smoothest sites at Lordsburg Playa).Histograms of z 0 values measured at each site are presented in Fig. 7a and c.As noted in Sect.2.2, a z 0 value was calculated for each 12 s interval for which a least-squares fit of u to ln z yielded a R 2 value of greater than 0.95. Figure 7 shows that z 0 values are approximately lognormally distributed.Sites that have higher-amplitude microtopographic variations at the 0.01 m scale (as measured by the average amplitude spectra in Fig. 5) have higher z 0 values.Aside from measurement error/uncertainty, there are two reasons for variance in measured z 0 values.The first is the fluctuating nature of turbulence itself.This source of variance can be reduced by averaging the wind velocities over longer time intervals before fitting to Eq. ( 1).The second source of variance comes from moving the hot-wire sensors to different locations around each site, thereby "sampling" different patches of microtopography.We found this second source to be the dominant source of variation based on the fact that z 0 values exhibit much greater variability over timescales of ∼ 1 h, i.e., the timescale over which the hot-wire sensors were moved around the landscape.
Values of mean z 0 for each site have a power-law dependence on H RMSE (Fig. 8a), i.e., where c = 6 ± 1 m −1 and b = 2.0 ± 0.1 and the uncertainty values represent 1σ standard deviations.Equation ( 6) is broadly consistent with the results of Nield et al. (2014) (Eq.2).The value of the exponent b we obtained is slightly higher than that of Nield et al. (2014), but such a difference is not unexpected considering that we are studying different playas.
There are several limitations with using H RMSE as the sole or primary predictive variable for z 0 .First, a nonlinear relationship between z 0 and H RMSE yields unrealistic values when applied outside the range of spatial scales considered here and in Nield et al. (2014).For example, using Eq. ( 6) with H RMSE values in the range of predicts z 0 values in the range of 6-54 m, i.e., values larger than any value ever measured.Playa surfaces rarely, if ever, have H RMSE values of 1-3 m, but many other landscapes (e.g., alluvial fans) do.Since the goal of this work is to use playas as model landscapes for understanding the multi-scale controls on z 0 above landscapes in general (not playas specifically), it is necessary for any empirical equation to predict reasonable results for a broad range of landscape types and a range of spatial scales beyond the specific range considered in the model calibration.Second, H RMSE values are problematic to use as the sole or dominant variable for predicting z 0 values because they contain no information about terrain slope.A topographic transect with a point spacing of 0.01 m can be "stretched" to obtain any slope value, with important consequences for flow separation and z 0 values.
Figure 8b plots the relationship between mean z 0 and S av , the mean slope computed between adjacent points at the 0.01 m scale, for the 10 study sites.This figure documents a systematic nonlinear relationship between z 0 and S av , suggesting that the nonlinearity between z 0 and H RMSE in Eq. ( 6) may reflect a dependence of z 0 on S av in addition to a dependence of z 0 on H RMSE .This hypothesis is consistent with Fig. 8c, which demonstrates that H RMSE values are highly correlated with S av values, i.e., that, in the playa surfaces we studied, playas with larger microtopographic amplitudes are systematically steeper.We would not expect such a correlation between amplitude and steepness to apply to all landform types because, as microtopography transitions into mesotopography and H RMSE increase from 0.1 to 1 and higher, slope gradients do not continually steepen without bound.If our goal is to understand the controls on z 0 values in landscapes generally, the data in Fig. 8 suggest that it is necessary to quantify the separate controls of amplitude and slope on z 0 values.This was the purpose of the CFD modeling described in the next section.

Computational fluid dynamics
To demonstrate the suitability of PHOENICS for modeling atmospheric boundary-layer flows and to establish that the effective roughness length depends on the microtopographic variability at multiple scales, we performed a numerical experiment using the central microtopographic profile measured at the Soda Lake smooth site as input (plotted in Fig. 9a).We measured a mean z 0 value of 4.6 mm from ve- locity profiles at this site.Figure 9b presents the velocity profiles predicted by the PHOENICS model for 2-D flow over the profile, following the procedures detailed in the Methods section.PHOENICS predicts an effective roughness length of 2.4 mm based on a least-squares fit of the velocity to the logarithms of distance above the ground from a height equal to twice the height of the dominant roughness elements to the top of the model domain.As such, the PHOENICS model predicts a z 0 value similar to the value we measured in the field (relative to the 4 order of magnitude variation in z 0 values we measured across the study sites).
To demonstrate that the z 0 value depends on microtopographic variability at multiple scales, we filtered the Soda Lake smooth profile diffusively to remove some of the smallscale (high-wavenumber) variability while maintaining the large-scale variability (i.e., the root-mean-squared variability of the filtered and unfiltered profiles is identical).Figure 9 plots the original profile, the filtered profile, and their amplitude and z 0n spectra.The z 0 values for the unfiltered and filtered cases are 2.4 and 0.15 mm, respectively, based on fitting the velocity profiles predicted by PHOENICS.That is, the filtered profile has a z 0 value more than an order of magnitude smaller than the original profile despite the fact that the amplitude of the large-scale microtopographic variations is the same as the original profile.Eq. (3) predicts 2.8 and 0.25 mm, respectively, for the z 0 values.The z 0 value decreases in the filtered case because steep slopes that trigger flow separation are significantly reduced at a wide range of scales by filtering, lowering the z 0 value.The results of this numerical experiment demonstrate that z 0 values depend on variability microtopographic variability at multiple scales.There is also a general theoretical argument that supports this conclusion.If one accepts that both the amplitude and slope of the microtopography influence the effective roughness length (which we will demonstrate below for the case of a sinusoid), it follows that there is no single Fourier mode that controls the effective roughness length, unless the topography is a perfect sinusoid.This is because the slope is a high-pass filter of the topography (i.e., the slope is proportional to k * a n , where a n is the Fourier coefficient) and hence is more sensitive to high-wavenumber components of the topography than the amplitude is.
Figure 10 presents color maps that illustrate the output of the CFD model for an example case (a = 0.05 m and S = 0.79 m/m).Figure 10a, which shows a color map of the turbulent kinetic energy, illustrates the growth of the internal boundary layer with increasing distance from the upwind boundary of the domain as the input velocity profile interacts with and adjusts to the presence of the microtopography.Figure 10b and c  zones of flow separation that occur in this example.These figures also illustrate the terrain-following and logarithmically spaced nature of the computational grid in the vertical direction.
Figure 11 plots the z 0 values computed from an analysis of the CFD-predicted wind profiles over sinusoidal topography for two different values of the sine-wave amplitude (a = 0.05 and 0.1 m) and for a range of values of the maximum slope S from approximately 0.001 to 2.0.For maximum slope values less than approximately 0.004, the z 0 value is equal to z 0g , as we would expect (the topography is effectively flat).As the slope of the microtopography increases, the wind field is increasingly perturbed by the roughness of the terrain.Eventually, flow separation is triggered and flow recirculation zones are created in the wakes of each bedform, further increasing z 0 values.For very steep slopes, i.e., S ∼ 0.4-1, z 0 values still increase with increasing slope but at a slower rate than for gentler slopes since the flow as already separated and additional steepening has only a modest effect on the spatial extent of flow separation and z 0 values.The nonlinear dependence of z 0 on S is well fit by a sigmoidal relationship of the form given by Eq. (3).Best-fit values are c 1 = 0.1, c 2 = 0.4, and c 3 = 2.0.
3.4 Fourier analysis of topography and a multi-scale approach to quantifying z 0 Using a minimization of the squared difference between the mean measured values of z 0 and the values predicted by Eq. ( 4) for all study sites, we found the optimal value of c 4 to be 1.5.Figure 12 plots z 0n values computed by Eq. ( 4) as a function of the natural wavenumber, k.The sum of all the z 0n values is the predicted value of z 0 for each surface.There is also value, however, in examining the dependence of z 0n on the wavenumber.The plot in Fig. 12 shows which spatial scales are most dominant in controlling the value of z 0 for a given landscape (see arrows in Fig. 12).On Lordsburg Playa, the only spatial scales that have non-negligible slope gradients are those at 0.01-0.3m.At the rougher sites, the dominant roughness elements are found at different scales, from 0.1-1 m (Soda Lake) to 1-10 m (Willcox Playa) to 0.3- 3 m (Death Valley).This plot also shows that in some cases there is one dominant scale of roughness elements (e.g., Soda Lake and Death Valley), while in others there are two or more scales that are equally dominant (e.g., Willcox Playa).
Figure 13 plots the z 0 values predicted by Eq. ( 4) vs. the mean measured values for the 10 study sites.Note that there appears to be only nine points plotted in Fig. 13 because two of the points (for Lordsburg smooth and Lordsburg intermediate) are nearly indistinguishable.The correlation between the logarithms of the predicted and measured mean z 0 values is quite good (R 2 = 0.991).Equation ( 4) is capable of predicting z 0 values to 50 % accuracy, on average, across a 4 order of magnitude range.
An alternative approach is to use the values of H RMSE and S av to estimate z 0 using Eq. ( 5).We found c 5 = 16 and c 6 = 2.0 to yield the highest R 2 value (0.978).Eq. ( 5) is thus a useful formula with an advantage of simplicity, but it is somewhat inferior to the multi-scale analysis of Eq. ( 4) based on its lower R 2 value.

Discussion and conclusions
The values of c 3 and c 2 respectively reflect the magnitude of the nonlinear increase in z 0 values as slope increases and the slope value where back-pressure effects begin to limit the rate of increase in z 0 with increasing slope.The values of c 3  3) and ( 4) do not utilize curvature, and hence neither equation can be the basis of a perfect method for predicting z 0 .It is likely that the only way to precisely estimate z 0 is to compute the actual flow field over the topography using a CFD model.Any other approach will likely involve some type of approximation.We propose that Eq. 4, while imperfect, yields a good approximation for z 0 values in real-world terrain (i.e., those with many Fourier coefficients contributing to z 0 ), based on the R 2 value of 0.991 we obtained.Eq. ( 5) provides an alternative for users who prefer its simplicity.Eq. ( 5) is not accurate for all possible S av values, since z 0 cannot increase without bound as S av increases.As such, Eq. ( 5) should only be considered applicable for microtopography with S av values less than approximately 0.15 m/m.
We developed and tested a new empirical formula for the roughness length, z 0 , of the fully rough form of the law of the wall that uses the amplitude and slope of microtopographic variations across multiple scales within a discrete-Fouriertransform-based approach.A sigmoidal relationship between z 0 and the amplitude and slope of sinusoidal topography developed from CFD model results was used to quantify the effects of each scale of microtopography on z 0 .The model was developed and tested using approximately 60 000 z 0 values from the southwestern US obtained over 2.5 orders of magnitude in distance above the bed.The proposed method is capable of predicting z 0 values to 50 % accuracy, on average, across a 4 order of magnitude range.This approach adds to our understanding of and ability to predict the characteristics of turbulent boundary flows over landscapes with multi-scale roughness.

Figure 1 .
Figure 1.Plots of synthetic (top) and real (middle and bottom) topographic transects illustrating the multi-scale nature of topography using natural playa surfaces as examples.

Figure 2 .
Figure 2. Aerial images of the study sites.

Figure 3 .
Figure 3. Photographs of the equipment used for measuring wind profiles.(a) Mast holding four hot-wire anemometers (left) and four cup anemometers (right; note that only the lowest three are visible) at the Animas intermediate study site.(b) Close-up photograph of the hot-wire sensors at the Soda Lake smooth site.For scale, note that the top hot-wire sensor is located 0.16 m above the surface in both photographs.

Figure 5 .
Figure5.Plots of the average amplitude spectrum, A, of 1-D transects of the microtopography of each site as a function of the natural wavenumber, k.The colors red, green, blue, and black are used to represent the Death Valley, Soda Lake, Willcox, and Lordsburg sites, respectively.Thicker curves represent rougher sites within each playa.

Figure 7 .
Figure 7. (a, b) Normalized histograms of z 0 values measured at each site and (c, d) probability distributions for each site, assuming z 0 values are lognormally distributed.

Figure 8 .
Figure 8. Plots of mean z 0 at each site vs.(a) H RMSE and (b) S av .(c) Plot of H RMSE vs. S av .

Figure 9 .
Figure9.Demonstration of the dependence of z 0 values on the multi-scale nature of microtopography.(a) Plot of a profile through the Soda Lake smooth site (thin curve).Also shown is the same plot with diffusive smoothing (thicker curve).Smoothing maintains the amplitude of microtopographic variations at large spatial scales (i.e., the amplitude spectrum is unchanged at large scales) but removes some of the small-scale (high-wavenumber) variability.(b) Plots of the mean velocity profiles predicted by PHOENICS over the original and filtered profile.(c) Amplitude spectra of the two plots in (a).(d) Contributions of each Fourier mode to the z 0 values for the two plots in (a).

Figure 10 .Figure 11 .
Figure 10.Illustrations of the output of the PHOENICS CFD model for the example case (with amplitude a = 0.05 m and maximum slope S = 0.79 m/m) of flow over a sinusoidal bed.(a) Color map of turbulent kinetic energy, KE.This map illustrates the growth of the internal boundary layer triggered by the effective roughness change as the input velocity profile (characterized by a grain-scale roughness z 0g ) interacts with and adjusts to the microtopography.The color vector maps in (b) and (c) illustrate the zones of flow recirculation that occur in the lee side of each bedform.

Figure 13 .
Figure 13.Plot of mean measured z 0 values vs. predicted values (using Eq. 4) for the 10 study sites.Error bars denote 1σ variations in the measured z 0 values.
Plots of mean wind velocity (normalized by the velocity measured at the highest sensor, located 2.8 m above the ground) (x axis) as a function of the natural logarithm of height above the ground (y axis).The colors red, green, blue, and black are used to represent the Death Valley, Soda Lake, Willcox, and Lordsburg sites, respectively.Within each playa, thicker lines are used to represent the rougher sites.Open circles indicate stations located within the roughness sublayer.These sensors were not used to calculate z 0 .