Aeolian shear stress ratio measurements within mesquite-dominated landscapes of the Chihuahuan Desert, New Mexico, USA

A field study was conducted to ascertain the amount of protection that mesquite-dominated communities provide to the surface from wind erosion. The dynamics of the locally accelerated evolution of a mesquite/coppice dune landscape and the undetermined spatial dependence of potential erosion by wind from a shear stress partition model were investigated. Sediment transport and dust emission processes are governed by the amount of protection that can be provided by roughness elements. Although shear stress partition models exist that can describe this, their accuracy has only been tested against a limited dataset because instrumentation has previously been unable to provide the necessary measurements. This study combines the use of meteorological towers and surface shear stress measurements with Irwin sensors to measure the partition of shear stress in situ. The surface shear stress within preferentially aligned vegetation (within coppice dune development) exhibited highly skewed distributions, while a more homogenous surface stress was recorded at a site with less developed coppice dunes. Above the vegetation, the logarithmic velocity profile deduced roughness length (based on 10-min averages) exhibited a distinct correlation with compass direction for the site with vegetation preferentially aligned, while the site with more homogenously distributed vegetation showed very little variation in the roughness length. This distribution in roughness length within an area, defines a distribution of a resolved shear stress partitioning model based on these measurements, ultimately providing potential closure to a previously uncorrelated model parameter. © 2006 Elsevier B.V. All rights reserved.


Introduction
Within southern New Mexico, USA, a desert shrub invasion has been occurring over the past 60 years, which is thought to be the result of the combined effects of overgrazing by cattle and soil erosion by wind and water (Schlesinger et al., 1990). The cause of the ongoing replacement of the grasslands has been linked to the preferential nutrient cycling of shrubs that separates and then eventually isolates soil nutrients into spatially distinct niches (Reynolds et al., 1999). This produces "resource islands" where small areas of shrubs are clumped together to use and trap nutrient rich soil. Okin and Gillette (2001) argued that the areas between the shrubs and the shrubs themselves have preferentially aligned themselves, as detected from remote sensing, with the dominant wind direction to preferentially grow in length, promoting soil erosion. Over time, and with relatively constant environmental conditions, these "streets" of bare soil tend to expand in both length and width as wind erosion continues to produce an environment that lacks almost any soil profile or soil nutrients (Okin and Gillette, 2001). The resulting landscape can have one of two end points: the creation of coppice dunes (which are presently developing in this region), or the complete desiccation and removal of vegetation leading to increased wind erosion until the underlying calcrete horizon is reached. This is a positive feedback dynamic from seasonally strong winds within a small directional range and little moisture availability, resulting in the unprotected soil being eroded by wind.
Vegetation or nonerodible elements on a surface aid to reduce the amount of overall soil loss by wind. By introducing vegetation in a systematic fashion, the amount of wind that reaches the surface can be reduced. Also, vegetation acts to physically block moving sediment, covers a portion of the surface, and because of its porous and flexible properties extracts momentum more efficiently than bluff barriers (Gillies et al., 2002). However, increasing the amount of vegetation on a surface increases the total shear stress and increases maximum stresses at the surface that can lead to a greater potential for wind erosion depending on the vegetation spacing and density (Musick et al., 1996).
The threshold at which a vegetated erodible surface will begin to erode by wind flow depends on many variables, including wind speed, height, width and spacing of vegetation, soil particle size, moisture content, and crust development (Wolfe and Nickling, 1993). However, roughness density, which is a measure of the size and number of roughness elements per surface area, has been shown to be a good descriptor of roughness that strongly influences sediment transport by wind (Gillies et al., 2006). By partitioning the shear stress between the erodible surface and the vegetation, Marshall (1971) expressed the sheltering capacity of vegetated surfaces as a ratio for a given roughness density. Through further development of this approach, Raupach (1992) produced a shear stress partitioning model that uses the geometric and aerodynamic properties of roughness of a surface. By utilizing a shear stress partitioning model to predict the amount of stress imposed by the wind on the soil surface, an estimation of soil loss and dust emissions is plausible.
The dynamics of surface shear stress reduction by vegetation are adequately predicted by the  drag partition model (King et al., 2005). However, the incorporation and validation for porous roughness elements in the model have not been previously satisfied because of field instrumentation limitations that result in many of the model parameters being estimated. Also, laboratory wind tunnel investigations that incorporate arrays of porous objects have primarily been restricted to using dowels or pegs instead of actual vegetation (e.g., Musick et al., 1996;Minvielle et al., 2003). Furthermore, the parameters within the model have not been fully explored to identify their spatial and temporal dependencies and their variations if the roughness elements are porous. Lastly, the dependence of the model on a roughness concentration variable originally applied to laboratory roughness in ordered configurations (e.g., staggered and solid) potentially limits its application to natural environments. This paper presents results from a field study conducted in New Mexico, USA, where direct measurements of the partitioning of stress over two types of vegetated surfaces were conducted. Through the use of instrumentation that has been extensively tested within the laboratory as well as in field studies, one of the objectives is to measure the spatial distributions of the stress budget within natural vegetation to provide insight towards developing a more comprehensive understanding of shear stress partitioning and aeolian sediment transport processes. The other main objective is to test the partitioning model of Raupach et al. (1993) for sparse vegetation by examining the spatial dependencies of its variables.

Boundary layer dynamics
The vertical wind speed profile has a characteristic form over a roughened surface that can be described by the "law of the wall": where u z (m/s) is the wind speed at a height z (m), κ is a constant of 0.4, z 0 (m) is the roughness length, d is the displacement height (m), and u ⁎ is the shear velocity (m/s). The shear velocity describes the transfer of time-averaged momentum in a turbulent, uniform, and steady flow over the surface (i.e., the surface and the roughness on it) and is related to the shear stress (τ) that the entire surface experiences as wind exerts a force over it through where ρ a is the density of the air. The force from the wind creates drag and lift forces on the sediment at the surface. If these combined forces exceed the gravity and binding energy forces acting on the surface particles, entrainment and subsequent transport will occur. In an atmospheric boundary layer, the wind is neither uniform nor steady because of the varying eddy sizes that represent the different length scales found within the internal boundary layer. This at-a-point unsteadiness and nonuniform velocity has been proposed as the dominant process that controls wind erosion and dust emissions (Gillette, 1999). The spatial and temporal maxima of shear velocity occurring on a surface can result in the threshold velocity of a surface being intermittently reached (Rice et al., 1999). Positively skewed Weibull distributions of surface velocities or shear stresses are typical of natural boundary layers even though time-averaged statistics are applied for sediment-transport-related calculations leading to an underestimated probability function of sediment transport (Schönfeldt, 2003). Research from a sediment transport aspect has explored this logistical problem and concluded that a shorter time-averaged function of shear velocity is required to capture the higher frequencies characteristic of skewed distributions (Stout, 1998;Wiggs et al., 2004). Once threshold is reached, a negative feedback or cascade process can occur, giving rise to sediment transport on the surface and possibly dust emission (Gillette et al., 1996). Therefore, for investigating the potential of wind erosion, the focus should be on the minimum shear stress that will provide enough force to entrain and maintain the transport of sand-sized particles (Rice et al., 1999), which act as the dominant driving mechanism for dust emissions (Shao et al., 1993).

Vegetation
Although vegetation acts to reduce wind erosion by simply covering a portion of the potentially erodible surface, its main effect is to extract momentum from the wind. By protruding into the air stream, vegetation increases the roughness length of the surface but also increases the zero-plane displacement height, therefore minimizing the shear stress at the surface on average, while increasing its variability. If sediment transport does occur, the vegetation can also act as a barrier to the moving sand particles. This blocking of sediment will significantly reduce sediment transport when the vegetation has a higher rigidity and increased coverage.
The coverage provided by vegetation in arid and semiarid environments varies with the supply of soil moisture and scattered storm events, creating localized patterns. Different types of vegetation also have distinct growing patterns that are sometimes circular (Larrea tridentata, creosote), dendritic (Flourensia cernua, tarbush), or linear (Prosopis glandulosa, mesquite), none of which are sim-ilar to the staggered array configuration used in laboratory experiments. The preferential growth of mesquite in linear forms in combination with the growth of coppice dunes is well documented to be a pattern that enhances wind erosion (Okin and Gillette, 2001). Complications arise when trying to parameterize this type of pattern into a model based on staggered arrays of roughness. The roughness density is one measure of describing roughness concentration and is defined as where n is the number of roughness elements of width b and height h per unit surface area (S). Although this form allows for multiple sizes of roughness, it is expressed as an average because the models that incorporate it essentially "spread" the roughness over a surface in a uniform or staggered array form because of the unknown complications associated with predicting surface stresses within random and clumped roughness configurations.

Drag partitioning
The total drag on a roughened surface from a fluid (F ) is defined as the sum of the drag on the intervening surface (F s ) and the drag on all of the individual roughness elements on that surface (F r ) (Schlichting, 1936). Marshall's (1971) laboratory study established that the partition was mostly dependent on the roughness density and only slightly dependent on the shape and arrangement of the roughness elements.
The basic understanding of the dynamic relationship between shear stress and roughness density led to the theoretical development of a model by Raupach (1992) and its final form in Raupach et al. (1993).  attempted to incorporate a spatial and temporal maximum shear stress component into the Raupach (1992) model that is based on the average shear stress ratio (R′ t ), by defining the partition of drag in terms of a ratio of maximum surface shear stress where the initiation of sediment transport could occur (τʺ s ) to the total shear stress over the surface with vegetation (τ) as the maximum shear stress ratio (Rʺ t ): where β is defined as the ratio of the drag coefficient of single surface-mounted roughness element to that of the surface without any elements, σ as the ratio of the basal to frontal area, and m is a value between 0 and 1 and is used to account for the spatial and temporal variations in the stress on the intervening surface. In this approach, β accounts for the shape of the elements and controls the partition of the drag (Raupach, 1992). The initial hypothesis of Raupach et al. (1993) was that a value of m = 0.5 is for a flat, erodible surface and a value approaching m = 1 is for an erodible surface that is quasistable with equilibrium bed topography. However, mathematically, a value of m = 1 equates to the average surface shear stress and an average shear stress ratio (R t ′), while a lesser value could be applied to represent a decrease in the sheltering ability of the roughness elements. Despite this definition for only the onset of sediment transport, measured shear stress ratios without sediment transport (R′ and Rʺ) are documented to be equivalent in value to those at threshold (Crawley and Nickling, 2003).

Study site
This study was conducted 40 km north of Las Cruces, NM in the northern Chihuahuan Desert on the Jornada Experimental Range (JER) of the USDA (http://usda-ars. nmsu.edu). The JER Agricultural Research Station is located on the Jornada del Muerto plain bordered by the San Andreas Mountains to the east and the Rio Grande valley to the west. The research conducted at JER was in an area monitored by the Long Term Ecological Research (LTER) program since 1989 to allow for historic comparisons of wind regime and sediment transport rates coupled with vegetation concentration data.
The Jornada LTER consists of 15 vegetationmonitoring sites of which the two that were observed for this study are dominated by P. glandulosa (honey mesquite). MNORT (UTM 332,360 E,3,610,376 N) and MWELL (UTM 326,307 E,3,608,925 N) were extensively monitored over the spring windstorm season (March to May). MNORT consists of developed coppice dunes typically ≥1 m in height. In contrast, MWELL had mesquite with heights ≤1 m with no coppice dune development. Scleropogon brevifolius (burrograss) and Yucca flaccida (yucca) were also present in observable quantities.
The vegetation roughness density at these sites was determined by measuring the height, longest axis, and intermediate axis (perpendicular to the longest axis) for each of the vegetation clusters within a 50 m by 50 m grid where the instruments were located (Table 1). Previously documented vegetation concentration data from Okin and Gillette (2001) and Gillette and Pitchford (2004) in addition to the measurements made for this experiment provided enough information to calculate λ.

Wind field
The experimental design at both MNORT and MWELL was similar to allow for direct comparison between the two sites and to further qualify the use of specific instrumentation. At each site, a 9.0 m anemometer tower and two 4.5 m anemometer towers with 7 and 5 R.M. Young cup anemometers, respectively, were deployed. The smaller 4.5 m towers were utilized for the calibration of the Irwin sensors as described in the next section. The anemometers (mounted on 1.0 m booms,) were spaced logarithmically with the bottom three cups of the main tower (9.0 m) matching the smaller towers (4.5 m). Therefore, the smaller tower anemometer heights were at 0.30, 0.60, 1.16, 2.28, and 4.5 m and the main tower anemometers were at heights of 0.30, 0.60, 1.16, 1.94, 3.23, 5.39, and 9.00 m. R.M. Young vanes were used for independent measurements of wind direction within the study area's flow field at 2.75 m on each of the small towers and at 9.00 m on the main tower. Positions of the towers are illustrated in Fig. 1 relative to the vegetation for both sites. All anemometers and vanes were sampled using Campbell Scientific 10X dataloggers sampling at 1 Hz and averaged at 10-min intervals.

Surface shear stress
In addition to anemometry for measuring the flow field near the surface, 16 Irwin sensors were deployed within both sites. The Irwin sensor is a differential pressure instrument mounted flush with the surface, and samples the pressure at the surface and at a height of 0.00175 m (Irwin, 1981). The Irwin sensors were sampled using 18 VDC electronic differential pressure transducers at each sensor, allowing a continuous 1 Hz sampling period with 10-min averaging procedure. The differential pressure measured with the Irwin sensor can be calibrated against surface shear stress when mounted flush with the surface. Its use within laboratory investigations is extensive (Irwin, 1981;Wu and Stathopoulos, 1994;  Monteiro and Viegas, 1996;Crawley and Nickling, 2003) and derived shear stress measurements show very good agreement with direct surface shear stress drag meters (Crawley and Nickling, 2003) and hot-wire anemometry (Irwin, 1981). These sensors were systematically positioned at the surface based on vegetation clusters to resolve the average surface flow field and also provide information on the position of the maximum and minimum shear stress areas within the area of the vegetation community being monitored ( Fig. 1). For an in situ calibration, a sensor was placed near each of the two 4.5 m towers. The surface shear stress (u ⁎S ) from the logarithmic profile, based on wind speed measurements for the bottom four anemometers, using Eqs. (1) and (2), was correlated with the differential pressure output of the Irwin sensor only when the wind direction was aligned with the maximum fetch for that tower. This assumed that the velocity profile within the sparsely vegetated area was characteristic of the sand surface roughness with a minimum 15 m fetch length. The exponential relationship between the shear stress values from the logarithmic profile and the pressure differential of the Irwin sensor at the small tower was then applied to all of the Irwin sensors at the site. This form of calibration was utilized for each storm event since the Irwin sensor's sensitivity to varying climatological conditions and follows the methodology of Gillies et al. (in press).

Vegetation
The vegetation characteristics of both MNORT and MWELL are detailed in Table 1. Mean and standard deviations of the heights of the vegetation along with the roughness densities are tabulated, distinguishing the differences between the two sites. MNORT is composed of coppice dunes with mesquite and has much larger roughness elements than MWELL, although the much larger spacing of the elements at MNORT as compared to MWELL produces a smaller difference in their vegetation concentrations (either expressed as percent cover or λ) than would be expected. Although the general spacing of the roughness elements can be deduced from these data, the spatial organization of the vegetation cannot be derived. As shown in Fig. 1, described by Gillette and Pitchford (2004), and statistically determined by Okin and Gillette (2001), the vegetation is distinctly organized at MNORT into linear forms aligned with the dominant wind direction (∼225°), while at MWELL the pattern has no directional dependence.

Wind field
Within the field study time period, two storm events were recorded at the MNORT location, while one stronger but shorter storm was recorded at MWELL. The duration of each of the storm events were N 24 h and, in the case of the first two events at MNORT, exceeded 48 h. Mean shear velocity, roughness length, and displacement height along with their standard deviations are displayed in Table 2, in addition to the maximum and minimum wind directions. Examples of typically logarithmic profiles with their bestfit linear regression are plotted in Fig. 2 for MNORT events #1 (A) and event #2 (B), and MWELL (C).
The wind profiles from the 9.0 m tower were filtered to ensure neutral conditions by selecting profiles where the mechanical turbulence would greatly exceed any convective effects. The procedure from Wolfe and Nickling (1996) was followed by removing any 10-min profiles when the lowest anemometer (0.30 m) did not exceed 3 m/s. Additionally, for the intervals when the R 2 value calculated from the logarithmic profile was b 0.98, the interval was not included in the available dataset. These steps insured that any aerodynamic variables calculated would be significant at the 0.1% level and that all profiles that were included were dominated by mechanical interactions (Sullivan and Greeley, 1993;Frank and Kocurek, 1994).
The 9.0 m tower wind profiles were then inspected graphically and were found to display upper and lower points within a transitional freestream and roughness sublayer, respectively (Fig. 2). Therefore, all calculations of the logarithmic profile were based on the 9.0 m tower using all but the top and bottom anemometers. This means that despite the anemometer tower being 9 m in height, all variables and statistics based on the wind profile from the 9.0 m tower are only based on the heights between 0.60 and 5.39 m and therefore will be referred to as the main tower.
After the above procedures were followed, an iterative method was applied to the anemometer heights to determine a displacement height for each of the 10-min averaged profiles. Eq. (1) was utilized to perform this analysis and to calculate the shear velocity and roughness length. The filtering of wind speed profiles by wind direction was applied after the limits set by the low wind speed requirements and regression coefficients of the wind velocity profiles were completed. This produced wind profiles for directions that are limited to those of the dominant wind regime.
The application of a displacement height for sparsely vegetated areas is debated because of the varying results found in previous field studies (e.g., Bottema, 1996;Wolfe and Nickling, 1996). By solving for displacement height in each of the 10-min acquisition periods the roughness length values were greatly reduced. This produced a greater agreement with the models of Lettau (1969) and Raupach (1994) for the relationship of z 0 /h and λ as the roughness lengths were reduced by up to half by the inclusion of a displacement height. Variations in the roughness length and the displacement height are plotted against the wind direction from the main tower in Fig. 3, illustrating the dependence of the wind velocity profile on the local configurations of roughness. A correlation exists between the roughness length and wind direction at MNORT that is explained by the vegetation arrangement displayed in Fig. 1, although the relationship is stronger for the higher average wind speeds of event #2. A preferential alignment of the vegetation at 225°(also the dominant wind direction) at MNORT has been previously documented through high resolution aerial photography in Okin and Gillette (2001) . Fig. 3 confirms the correlation of vegetation distribution and compass direction with increasing measured z 0 on either side of the dominant wind direction for MNORT (∼ 225°). The relationship at MWELL is weaker, dependent only on the geometry of the vegetation and shows only a slight increase in z 0 with wind direction because of the homogenous vegetation distribution at this site. The 95% confidence limits of the average z 0 has been calculated following the technique of Wilkinson (1984) resulting in ±5.3, 11.4, and 10.4% for MNORT event #1 and #2, and MWELL, respectively. This magnitude of the confidence limit is relatively low for this calculation and strengthens the correlated pattern of z 0 with the directionally aligned vegetation in Fig. 3.

Surface shear stress
The results from the Irwin sensors will be described for each event individually to examine temporal variations of the surface shear stress and then as a group to inspect spatial variations later in this section as part of the discussion of total stress. Table 3 displays, for all Irwin sensors, the number of 10-min intervals recorded, the averages, standard deviations, and maxima for events #1 and #2 (MNORT) and event #3 (MWELL), while histograms of selected individual Irwin sensors are displayed in Fig. 4 along with their higher-order statistics. The histogram was utilized to further investigate the temporal distribution of the signal record since the variations in the third-and fourth-order statistics suggested the temporal distribution of surface shear stress was not normal. The shape of this temporal distribution of surface shear stress, in particular the frequency of occurrence of the higher values is important for the initiation of sediment transport. Event #1 at MNORT recorded 9 out of the 13 sensors for more than thirty 10-min intervals individually within the event. Missing data from the time record were either from the sensor output being below its offset for that event or, in some cases, the sensor had recording problems resulting in it being omitted from that event. The results are only statistically significant for those nine sensors because of the number of observations, so the others are not included in the remaining analysis. The temporally averaged values for the surface shear stress show similar values throughout the measured area (average value of all sensors is 0.09 N/m 2 ) with coefficients of variation no more than 30%. However, the maximum values for the duration of the event vary from over double the mean value to a minimum 18% increase over the mean. The dissimilarity in the distributions of the values led to the investigation of higher order statistics as the majority of records are not normally distributed. Higher measures of skewness and kurtosis indicate that a larger amount of the distribution is contained within higher frequencies than  would be assumed for a normally distributed record. Although this will be reflected in a higher average value, associated standard deviation values will not fully describe the distribution. Therefore, a measure of skewness (standard error of skewness) that is normalized by the sample number is utilized to measure the proportion of values contained within the tails of the distribution. The standard error of skewness is defined as 2(6/i) 0.5 where i is the number of samples in the time series. When a value of skewness is twice or greater than the standard error of skewness it is statistically significant.  very small range, and the kurtosis value reflects this by being almost zero.
The surface shear stress result of event #2 at MNORT had 12 recording Irwin sensors with a total of 1905 10-min samples (Table 3B). Similar to the first event at this site, the coefficient of variation does not exceed 30% for event #2, but in contrast, the time averaged surface shear stress for all the Irwin sensors during event #2 is slightly higher (0.22 N/m 2 ). The maximum values are much higher relative to the average values for each sensor, ranging from approximately ∼20% to 130% higher. Of the 12 Irwin sensors that were in operation for this event, nine have positive skewness and kurtosis beyond 2 standard errors, and only one sensor has both negative skewness and kurtosis values. Fig. 4B displays the histograms from Irwins A, B, C, F, J, and L to further illustrate the statistical results. Irwin sensors A, B, and L have high positive skewness as well as a leptokurtic kurtosis evident by the tail of values toward the maximum. Irwin sensors C and J have similar statistical results; however, the histogram of Irwin C reveals a mode greater in value than its mean while maintaining its positive skewness. Irwin J in contrast, has many intervals of similar frequency in addition to a positively skewed tail. Lastly, Irwin F has a very symmetric distribution despite its negative skewness and platykurtic form.
Spatially, the results from the two events recorded at MNORT display a predominantly positively skewed distribution of surface shear stress throughout the field site except for a small portion in the center. If Irwin Sensors A and G are, because of their position, assumed to represent the beginning of the "streets" within the study area (see Fig. 1), a descriptive interpretation would explain the surface shear stress distributions as becoming more positively skewed with decreasing "street" width and secondarily, increasing skewness with street length. Although the pattern is not very strong in this instance, it could explain sediment transport patterns exhibited by the vegetation growth and "street" development by linking positive skewness in the shear stress signal with increased erodibility potential. Event #3 at MWELL produced a smaller number of results because of instrument errors, recording only five Irwin sensors, although a total of 479 10-min data intervals were collected (Table 3C). The average surface shear stress values are much greater for this site than MNORT (0.35 N/m 2 ), although producing maximum values only 50% larger than its average. Additionally, the coefficient of variation is lower despite the higher absolute values, to a maximum of 20% for Irwin C. Of the five sensors available, four showed positive skewness values N 2 standard errors of skewness in the surface shear stress distributions while none of them were statistically displaced from a neutral kurtosis. The histograms (Fig.  4C) for this site show that Irwin sensors B and C have similar surface shear stress distributions to the sensors at MNORT as the median is to the left of the mean and the skewed portion is well pronounced toward the higher values. Irwins I and K show surface shear stress distributions that are less skewed, with Irwin K being the only sensor without a significant positive skewness in addition to an almost leptokurtic kurtosis distribution.
Although the number of sensors at MWELL that could be used is lower than at MNORT and they are spaced further apart, they do show that the surface shear stress has a tendency toward being normally distributed in contrast to the dominantly positively skewed distribution of MNORT. Additionally, higher ratios of the maximum to the average surface shear stress values are found within the larger roughness elements associated with MNORT than those at MWELL. However, positively skewed distributions of surface shear stress are more commonly associated with higher λ values because of the reduced average surface shear stress but greater total shear stress. The difference in the organization of the vegetation between these two sites provides an explanation for the measured properties that are in opposition to those typically observed.

Shear stress ratio (SSR)
The SSR calculated with the temporally averaged and maximum surface shear stress from all the Irwin sensors and the total shear calculated from the main tower logarithmic profile are presented in Table 4 along with their associated statistics. Individual Irwin sensor SSRs indicate the spatial variation amongst the field sites, which predictability echo the values discussed in Section 4.3, because the total shear stress measurements from the main tower were utilized for all of the calculations. The second event at MNORT displayed a higher SSR value on average in comparison to event #1, although it did not produce larger maximum values the uniqueness of each wind event may be a function of varying turbulence properties. Comparing the values from MNORT and MWELL shows that the maximum and average values of MWELL are slightly higher than those at MNORT despite the 20% increase in λ.  duration of the event. The regression is forced through zero to conform with the original definition of SSR (as τʺ→0, τ′→0)  and displays a strong linear relationship slightly above the 1:1 line (R 2 values of 0.88, 0.95, and 0.94 for events #1, #2, and #3, respectively). The slope values calculated for the MWELL event are similar to values for the events recorded at MNORT, which seem to indicate no discernable trend with changing roughness density. The relationship in Fig. 5 indicates the presence of a convincing trend between the maximum and average values throughout the range of values and suggests a possible predictive tool for estimating the threshold of sediment transport.

Wind field
The variations in the wind field of z 0 and d with wind direction have been attributed to the variations in surface roughness with compass direction in Section 4.2, although other plausible explanations exist. The wind field data were calculated based on the use of wind speed profiles from a tower (i.e., the main tower) of sufficient height for the majority of desert vegetation roughened landscapes based on the fetch and vegetation height (Bottema, 1996). The shear velocities and roughness lengths reported for the surfaces in this study are within the standard deviation limits of previous measurements from 15.0 m towers also located at these sites as part of the LTER monitoring program (D. Gillette, personal communication, 2005). Furthermore, in an environment where the vegetation has been shown to be preferentially aligned with the dominant wind direction (MNORT), the surface appears rougher at an angle other than the dominant one. Fig. 3 shows that the roughness length is dependent on the compass direction when estimated within sparsely vegetated surfaces that have a directional component to them (MNORT) and even show a slight dependence for a surface that shows almost no directionality (MWELL). This indicates that single tower measurements are not satisfactory for providing an average value for surfaces of this nature and that multiple anemometer towers or multiple acquisition points need to be implemented to better assess boundary layer characteristics (Zobeck et al., 2003). The detailed comparisons of total shear stress in addition to z 0 that are presented in this paper now become site specific because of the directionality of the wind data. However, the surrounding areas that are assumed to be composed of similar vegetation development will contain the same attributes. This has implications for the roughness density calculation as well as the wind direction independent nature of the Raupach et al. (1993) model, but this is not new to the agriculture literature that incorporates direction-and height-dependent calculations into wind erosion models (Woodruff and Siddoway, 1965;Bradley and Mulhearn, 1983). Even though wind events in directions other than the dominant one are less likely to produce wind erosion at JER, the wind climate has forced the evolution of the landforms preferentially promoting desiccation of the soil surface (Gillette and Pitchford, 2004).

Surface shear stress
The in situ measurement of surface shear stress within both of these field sites allows for an in-depth examination of the effect that vegetation has on the distribution of wind-generated shear stress at the surface. Previously reported measurements of surface shear stress distributions have been restricted to laboratory wind tunnels.
The model of Rice et al. (1999) states that the highenergy skewed side of an impact energy distribution is more influential than the average in determining the transport potential. This conceptual model can be extended for the erosive potential of surface shear stress with or without sediment transport by assuming that the higher energy tail of the surface shear stress distributions is the cause of the erosion. The positively skewed distribution of surface shear stress explains the intermittency of sediment transport associated with transportlimited surfaces that appear to be non-erodible based on averaged measurements. Previous measurements of surface shear stress distributions with Irwin sensors within a staggered array by Crawley (2000) showed that although the distribution was positively skewed in most areas not all were statistically significant. The results from this study detail the statistically significant skewness values for the majority of the Irwin sensors, and in the case of MNORT an increase in skewness with fetch is suggested. When the array of cylinders used by Crawley (2000) is compared to vegetation arrays, the cylinders produce lower frequency eddies because of their bluff nature. However, their modelled size and higher wind velocities (although similar Re numbers) could explain the higher frequency eddies that create positively skewed distributions of the surface shear stress.
Through the detailed analysis of the shear stress distributions obtained with the Irwin sensors, the onset of sediment transport can be related to some function of the maximum surface shear stress. The linear relationship between the maximum and average surface shear stress above 1:1 is well established for all events (Fig. 5), and seems to have no dependency on roughness density, similar to the findings of Crawley and Nickling (2003). However, the data that form the maximum and minimum relationship (Fig. 5) in the study by Crawley and Nickling (2003) are quite variable, and upon further inspection there is a trend in their data that indicates an increase in the slope of the relationship between the maximum and average surface shear stress with increasing roughness density. This is anticipated because, as the roughness becomes denser, the amount of stress reaching the surface decreases and the relative difference in the total shear stress and the surface stress will increase (Musick et al., 1996). The difference in the roughness densities presented in this study, however, do not span the range that was reported by Crawley and Nickling (2003), which does not allow us to observe this trend.
A discrepancy exists between the values of surface shear stress in this study and that reported by Wyatt and Nickling (1997) and Crawley and Nickling (2003). The calibration process utilized in this study allowed for an in situ regulation of the measurements by the Irwin sensors because of their documented influence by climate variables, such as pressure, temperature, and humidity. In comparison, the calibration process of Wyatt and Nickling (1997) was done a priori within a laboratory wind tunnel and therefore cannot be considered for this comparison, while the investigation by Crawley and Nickling (2003) was within a wind tunnel allowing for a similar calibration to the one in this study. The absolute values of measured surface shear stress for only slighter higher roughness density values in this study are double those recorded in the wind tunnel by Crawley and Nickling (2003) even though both were using Irwin sensors. Although, the Reynolds numbers for this study are orders of magnitude larger than Crawley and Nickling (2003) or Wyatt and Nickling (1997), making the differences in surface shear stress values attributable to either the different roughness conditions or a different magnitude of wind velocities.

Shear stress ratio
The measurements made using the Irwin sensors and the main tower at both sites provide not only the measured shear stress ratio but can be utilized to predict components of the Raupach et al. (1993) model.
The SSR values measured in this study for all events are plotted against roughness density in Fig. 6, along with the field data of Musick and Gillette (1990), Wolfe and Nickling (1996), Wyatt and Nickling (1997), Baas (1998), andLuttmer (2002). The average and maximum surface shear stress-based ratios are included and fit the trend observed in previous investigations. The maximum value-based SSR does not differ from the average value SSR by a significant amount when compared to the range of values previously reported. However, the distributions in Fig. 4 support the conclusion that the average value does not represent the driving force for sediment transport.
The values chosen for parameterizing the Raupach et al. (1993) model were initially based on the geometric properties measured onsite in the case of λ and σ, while the value for β was estimated from the literature. The drag coefficient of the surface (C S ) was directly calculated by the shear velocity calculated from the smaller anemometer towers (u ⁎S ) referenced to the velocity at the average height of the roughness elements at each site (U h ) using the relationship below for the wind direction aligned with the 'streets': While the element drag coefficient (C d ) used was inferred from the value measured by Gillies et al. (2000) for an unobstructed porous shrub (Table 5), the values selected (although much higher than solid objects with this shape) are justified by the size and porous nature of the vegetation at both sites. In the case of MNORT, the vegetation itself was relatively flat and long but was raised vertically by the development of coppice dunes. This has the effect of creating a smaller shrub in a faster air stream as well as a bluff body that ultimately provides greater momentum extraction. The vegetation at MWELL was, on average, in between the sizes reported by Gillies et al. (2000) and was assigned a weighted value.
To validate the Raupach et al. (1993) model, an iterative method was used to solve for β based on the values of σ and λ presented in Table 5. The average and standard deviations of these values are also presented in Table 5 and result in β values being almost half of the values reported in the literature for both porous and solid roughness. For MNORT, the two events produced different results with the second event predicting a higher value of β in addition to a much larger distribution, possibly signifying that the drag coefficient of either the surface or vegetation is not constant. The value of β calculated for MWELL is closer to the assigned value with the standard deviation considered, although it was also under predicted. This suggests that element drag coefficients of 0.54 ± 0.22 and 0.69 ± 0.40 should represent the roughness at MNORT, while an element drag coefficient of 0.44 ± 0.15 represents the average roughness at MWELL. These values are less than the proposed values, but further calculation using Eq. (4) to solve for m produces errors associated with the model's limits. For this reason, the continuing analysis utilizes the initially proposed β values.
The Raupach et al. (1993) model can be utilized for predicting the m parameter when using Eq. (4) with the maximum measured surface shear stress. Previously reported values for this parameter are 0.14 to 0.18 from Wyatt and Nickling (1997) and 0.53 to 0.58 by Crawley and Nickling (2003). In this study, the estimated values of m are 0.19 and 0.28 for MNORT and 0.28 for MWELL with standard deviations of 0.13 to 0.21 and 0.12, respectively. These values were calculated based on each 10-min interval of matching Irwin sensors and 9.0 m tower data providing an extensive range of SSR values. The larger distribution of the values at MNORT in comparison to MWELL could be attributed to the larger amount of 10-min intervals recorded, but upon further inspection, the m parameter is highly correlated with the roughness configuration.
In Fig. 7A, B, C, the calculated m values are plotted against the wind direction for the 10-min interval of the SSR, on which they were based, and additionally the roughness length is also plotted against wind direction on a secondary y-axis. Fig. 7D, E, F and the correlation coefficients between m and z 0 of 0.48, 0.74, and 0.86 for events #1, 2, and 3, respectively, show that the calculated value of m is highly influenced by the geometric properties of roughness. This solution is a result of the variation in the measured shear stress ratio values showing a  prior dependency on z 0 . However, this pattern transforms the modelled shear stress ratio from values independent of the vegetation arrangement to values with an m that describes the spatial variation with wind direction. This is similar to the original definition by Raupach et al. (1993) that suggested the value of m accounted for the spatial distribution of surface stress characteristics but displayed a reduced complexity to the solution of m if a function could be solely associated with z 0 . Although a high degree of scatter in the relationship of z 0 and m at MNORT yields regression coefficients of 0.23 and 0.57 for events #1 and #2 (Fig. 7D, E), and a linear pattern is observed with slopes of 7.1 and 9.0, respectively. At MWELL, the relationship of z 0 and m also displays a linear relationship, but with a coefficient of 0.75 and a slightly greater slope value of 7.3 (Fig. 7F). This result could have positive implications for improving the practical application of the Raupach et al. (1993) model by reducing the uncertainty previously associated with assigning a value to m for different surfaces. This relationship seems to hold reasonably well for these two sites and, although they are within a small range of λ, they have different vegetative characteristics and represent two types of roughness density distributions possibly suggesting a general relationship between z 0 and m.

Conclusions
Traditional meteorological instrumentation along with Irwin sensors were used to resolve the surface shear stress ratio at two sparsely vegetated sites at the Jornada Experimental Range in New Mexico in two Fig. 7. The model-derived m (primary y-axis) and calculated roughness length (z 0 ) (secondary y-axis) are plotted against the wind direction for event #1 (A) and event #2 (B) at MNORT, and event #3 (C) at MWELL. Calculated values of m plotted against roughness length (z 0 ) for event #1 (D) and event #2 (E) at MNORT and for event #3 (F) at MWELL. different spatial distributions of mesquite communities. The measurements indicated that the aerodynamic properties of one of the surfaces (MNORT) were attributable to the alignment of the vegetation with the dominant wind direction. Although a large range of roughness lengths were recorded for that site, it shows good agreement with previous models, suggesting that the values obtained for z 0 were appropriate for the given roughness density.
The strong directional dependence of z 0 associated with vegetation alignment and correlated with the dominant wind direction at MNORT could be exaggerated because only one recording station was implemented. Although single tower measurements are commonly used in the logarithmic-profile calculation over vegetated surfaces, the results of this study show that over a site with heterogeneously distributed vegetation, this approach may be insufficient to capture the aerodynamic properties of an entire surface with similar properties. This suggests further experiments should involving a network of anemometer towers, comparing at-a-point measurements to surface averaged measurements, to strengthen or dismiss this hypothesis.
The measurements of surface shear stress showed highly positively skewed temporal distributions and displayed a slight increase in skewness when the distance between vegetation perpendicular to the wind direction was reduced at MNORT, while no such correlation existed at MWELL. The more regularly roughened site of MWELL, however, did show a greater homogeneity of surface shear stresses at each sensor in comparison to those measured at MNORT with its more ordered roughness.
The maximum surface shear stress values were utilized in the SSR analysis and displayed a very good agreement with the SSR based on average values although no roughness density-based trend was apparent. This relationship is also shown by Crawley and Nickling (2003); however, a discrepancy in the absolute values of the surface shear stress values exists. The difference in surface shear stress values is attributed to the scale and porous nature of the roughness in addition to the magnitude of wind speeds measured in this study, but is of little consequence as the SSRs are similar for both studies for the same roughness density. This may have implications on modelling vegetation with wind tunnels, although some disagreement may be attributed to the Irwin sensor calibration method utilized.
The modelled values of the SSR derived from the Raupach et al. (1993) model produced average values of 0.19, 0.28, and 0.28 for the three recorded events for the solution of m. However, a relationship was found between the roughness length of the surface and the solution of m for a given wind direction. This linear relationship showed an increase in m with an increase in roughness length and was shown to be in good agreement at MWELL, although a weaker correlation was found with the data from MNORT. This result contributes insight into the credibility of using a parameter to quantify the maximum surface shear stress and possibly a way of estimating its value upon further measurements with a larger variation in surface roughness density.