Spatial Distribution of Shrubs Impacts Relationships among Saltation, Roughness, and Vegetation Structure in an East Asian Rangeland

Vegetation influences the occurrence of saltation through various mechanisms. Most previous studies have focused on the effects of vegetation on saltation occurrence under spatially homogeneous vegetation, whereas few field studies have examined how spatially heterogeneous cover affects saltation. To examine how spatial heterogeneity of vegetation influences saltation, we surveyed the vegetation and spatial distribution of shrubs and conducted roughness measurements at 11 sites at Tsogt-Ovoo, Gobi steppe of Mongolia, which are dominated by the shrubs Salsola passerina and Anabasis brevifolia. Saltation and meteorological observations were used to calculate the saltation flux, threshold friction velocity, and roughness length. The spatial distribution of shrubs was estimated from the intershrub distance obtained by calculating a semivariogram. Threshold friction velocity was well explained by roughness length. The relationships among roughness, saltation flux, and vegetation cover depended on the spatial distribution of shrubs. When the vegetation was distributed heterogeneously, roughness length increased as the vegetation cover decreased, and the saltation flux increased because the wake interference flow became dominant. When the vegetation was spatially homogeneous, however, the saltation flux was suppressed even when the vegetation cover was small. These field experiments show the importance of considering the spatial distribution of vegetation in evaluating saltation occurrence.


Introduction
Dust is a serious environmental problem in East Asia, affecting not only human and livestock health, but also agriculture, grazing, land degradation, climate change, and biochemical cycles [1][2][3][4]. Dust processes can be divided into three stages: emission, transport, and deposition. Among these, the emission stage is complex and poorly understood because of the unique landscape features of arid regions [3,5]. Dust emission has been divided into erosivity, which depends on wind shear stress, and erodibility, which is affected by multiple elements of the land surface [6,7], including soil particle size, soil moisture, salt accumulation, snow cover, and vegetation. Dust emission processes from a vegetated ground surface, however, are not yet fully understood.
Dust emission is predominantly driven by saltation bombardment and aggregate breakup during transport events [8,9]. Vegetation affects saltation by three mechanisms [10]: (1) it covers the ground surface, protecting it from the wind; (2) it extracts momentum from the wind acting on a surface, and (3) it traps sediment. Prior studies have examined the influence of vegetation cover on dust emission with a focus on the first mechanism [11], and various vegetation cover thresholds for wind erosion to occur have been proposed [12][13][14][15]. The mitigating effect of vegetation on wind erosion has generally been described by the drag partition theory, in which plants are treated as roughness factors (aerodynamic Land 2021, 10, 1224 2 of 18 roughness length z 0 or lateral cover λ) on the land surface. In models applying drag partitioning schemes [11,16,17] in which plants are treated as lateral cover λ [18], wind momentum is divided into two parts, one acting on the ground surface and the other acting on roughness elements, including plants [19]. The third mechanism, the trapping effect of plants, depends mainly on the height of the vegetation [20].
Although the effects of vegetation on sediment transport have traditionally been explained in terms of roughness factors, vegetation cover, and the height of the vegetation community, it has recently been pointed out that the spatial heterogeneity of vegetation must also be considered [3,[21][22][23][24][25]. For example, in arid lands, even if the vegetation cover exceeds a certain threshold, sediment transport can still occur when the wind is strong [26,27]; thus, failure to take account of the effect of the spatial heterogeneity of vegetation may lead to underestimation of wind erosion [3,23]. Heterogeneity of dryland vegetation is characterized by patchy vegetation, with bare ground between vegetation patches. This structure enhances sediment transport, and connectivity has been proposed as a concept to explain these structures and functions [28,29]. Okin and Gillette (2001) [21] suggested that in a mesquite shrubland where the shrub distribution is anisotropic, sediment transport can be larger than that estimated from the vegetation cover when intershrub gaps (which they called streets) are large. Moreover, physically based models generally assume that roughness elements, expressed in terms of z 0 or λ, are distributed homogeneously across the ground surface and thus do not take into account the spatial heterogeneity of vegetation (e.g., the telephone pole problem [22]). Webb et al. (2021) [25] argued that it is important to consider the spatial arrangement of vegetation in assessing the susceptibility of landscapes to sediment transport. Okin (2008) [22] has proposed a model that incorporates intershrub distance to account for spatial heterogeneity of vegetation, and this model has outperformed previous models [30], but few field studies evaluating the effects of vegetation on sediment transport have considered spatial heterogeneity. Therefore, field studies that take account of the spatial heterogeneity of vegetation are needed to further elucidate the dust emission process, thereby contributing to the development of sustainable countermeasures against land degradation and the improvement of the accuracy of numerical models. In this study, we quantitatively evaluated the effect of various vegetation factors, including the spatial distribution of shrubs, on the occurrence of saltation in a shrub-steppe in Tsogt-Ovoo, Mongolia, a dust emission hotspot. Vegetation factors include vegetation cover and the length of gaps between shrubs (intershrub distance). Vegetation cover is the proportion of the ground area that is directly sheltered by plants and dead leaves, but this metric assumes that vegetation is spatially homogeneous. The intershrub distance normalized by the vegetation height is thus used as an indicator of the spatial distribution of vegetation. The intershrub distance is an important predictor of vegetation sheltering effects [25]. Although saltation is mitigated when vegetation cover exceeds a certain threshold, differences in the intershrub distance lead to different sediment transport rates among sites with similar vegetation cover [25]. Therefore, we hypothesized that the interaction of these two factors influences sediment transport. In order to focus on the effects of vegetation cover and the spatial distribution of shrubs on saltation occurrence, this study was conducted at sites with similar community species composition and community height. Vegetation height has a significant effect on the size of the sheltered area downwind of the vegetation [22,25]. Although vegetation shape, porosity, and flexibility are also factors that influence the protective effect of vegetation, because we considered intershrub distance to be more important [3,25], we excluded the effects of those factors from the analysis by selecting sites with similar species composition.
We examined the relationship between the occurrence of saltation and vegetation factors by observing multiple sites in the field, taking into account the spatial distribution of shrubs. We found that the threshold friction velocity was linearly related to the roughness length z 0 /h. When the intershrub distance was large, wake interference flow was dominant, and roughness length tended to increase with decreasing vegetation cover. When the intershrub distance was small, vegetation was relatively homogeneous, the roughness Land 2021, 10, 1224 3 of 18 length did not vary much with changes in vegetation cover, and the occurrence of saltation was suppressed even when the vegetation cover was small. These results suggest that it is important to consider the spatial distribution of shrubs when evaluating the effect of vegetation on wind erosion in an arid rangeland dominated by shrubs.

Study Area
Our study area was in Tsogt-Ovoo sum (district), Mongolia, in northeast Asia. Tsogt-Ovoo is located in the northern part of the Gobi Desert, an arid region known for frequent dust storms [1,31] (Figure 1). According to SYNOP data at Tsogt-Ovoo station from 2000 to 2020, the air temperature ranged from −10.0 to 21.2 • C, and the annual average precipitation was 189.8 mm. The study area is located in a valley, and the vegetation community changes stepwise from upland to lowland sites. To examine the occurrence of saltation in the shrub-steppe, we selected an observation area located on two pediments (north and south pediments), where the vegetation community is dominated by Salsola passerina Bge. (S. passerina) and Anabasis brevifolia C.A. Mey (A. brevifolia) ( Figure 2). These two species, which are both small subshrubs of the Amaranthaceae family, are distributed on clayey-pebbly soils of submountain plains. These shrublands provide excellent grazing for camels, sheep, and goats, the main livestock in the study area [32]. We conducted vegetation surveys and carried out saltation and meteorological observations at 11 sites in the study area from late April to early May 2019 ( Figure 1). dominant, and roughness length tended to increase with decreasing vegetation cover. When the intershrub distance was small, vegetation was relatively homogeneous, the roughness length did not vary much with changes in vegetation cover, and the occurrence of saltation was suppressed even when the vegetation cover was small. These results suggest that it is important to consider the spatial distribution of shrubs when evaluating the effect of vegetation on wind erosion in an arid rangeland dominated by shrubs.

Study Area
Our study area was in Tsogt-Ovoo sum (district), Mongolia, in northeast Asia. Tsogt-Ovoo is located in the northern part of the Gobi Desert, an arid region known for frequent dust storms [1,31] (Figure 1). According to SYNOP data at Tsogt-Ovoo station from 2000 to 2020, the air temperature ranged from −10.0 to 21.2 °C, and the annual average precipitation was 189.8 mm. The study area is located in a valley, and the vegetation community changes stepwise from upland to lowland sites. To examine the occurrence of saltation in the shrub-steppe, we selected an observation area located on two pediments (north and south pediments), where the vegetation community is dominated by Salsola passerina Bge. (S. passerina) and Anabasis brevifolia C.A. Mey (A. brevifolia) ( Figure 2). These two species, which are both small subshrubs of the Amaranthaceae family, are distributed on clayeypebbly soils of submountain plains. These shrublands provide excellent grazing for camels, sheep, and goats, the main livestock in the study area [32]. We conducted vegetation surveys and carried out saltation and meteorological observations at 11 sites in the study area from late April to early May 2019 ( Figure 1).

Field Survey
We established 11 sites (each 20 × 20 m). Insofar as possible, we tried to select observation sites that would be representative of the vegetation upwind of dust source areas, which could be several to tens of kilometers from the observation sites [33]. At each site, three piezoelectric particle counters (PCs; ud-101, Chuo Kosoku, Tokyo, Japan) were used to count saltation particles at 10 s sampling intervals. The PCs were installed at the same height (approximately 0.1 m) but at different, randomly selected points, and the results from the three PCs at each site were averaged. Although the ud-101 is not directional, we mounted each PC on a rotating platform equipped with wind vanes so that it would always face upwind ( Figure 3). Therefore, we observed saltation whenever there was wind, regardless of its direction. To obtain profiles of wind speed and direction and air temperature, three portable weather meters (Kestrel Weather Meter 5500, Nielsen-Kellerman, Boothwyn, PA, USA) were installed at different heights (approximately 2.1, 1.6, and 1.1 m; sampling interval, 20 or 30 s). PCs and weather meters were installed on the southeast side of each site because we expected northwestern winds to predominate. We moved the instruments from one set of sites to another set of sites after a dust event passed because of the limitation of the number of instruments; we observed dust events at 3 or 4 sites at the same time.

Field Survey
We established 11 sites (each 20 × 20 m). Insofar as possible, we tried to select observation sites that would be representative of the vegetation upwind of dust source areas, which could be several to tens of kilometers from the observation sites [33]. At each site, three piezoelectric particle counters (PCs; ud-101, Chuo Kosoku, Tokyo, Japan) were used to count saltation particles at 10 s sampling intervals. The PCs were installed at the same height (approximately 0.1 m) but at different, randomly selected points, and the results from the three PCs at each site were averaged. Although the ud-101 is not directional, we mounted each PC on a rotating platform equipped with wind vanes so that it would always face upwind ( Figure 3). Therefore, we observed saltation whenever there was wind, regardless of its direction. To obtain profiles of wind speed and direction and air temperature, three portable weather meters (Kestrel Weather Meter 5500, Nielsen-Kellerman, Boothwyn, PA, USA) were installed at different heights (approximately 2.1, 1.6, and 1.1 m; sampling interval, 20 or 30 s). PCs and weather meters were installed on the southeast side of each site because we expected northwestern winds to predominate. We moved the instruments from one set of sites to another set of sites after a dust event passed because of the limitation of the number of instruments; we observed dust events at 3 or 4 sites at the same time. For the vegetation survey, we established nine permanent survey quadrats (each 1 × 1 m), spaced at regular 5 m intervals, at each site. In each quadrat, the total vegetation cover and the height of the vegetation community were recorded. In addition to the dominant species S. passerina and A. brevifolia, shrubs such as Potaninia mongolica Maxim, perennial forbs such as Allium polyrhizum Turcz. ex Regel., and annual forbs such as Artemisia pectinata Pall. were observed. High precipitation in summer 2018 led to high growth of annual forbs, and many dead plants and litter remained at each site were trapped by shrubs such as S. passerina and A. brevifolia and formed mounds. The vegetation cover in each quadrat was calculated as the proportion of the ground surface covered by all living and dead plants, and the total vegetation cover was calculated as the average vegetation cover in all quadrats. This remaining dead vegetation cover was included in vegetation cover because dead plant detritus also contributes to the roughness of the ground surface and affects saltation occurrence [34,35]. In addition, we calculated community height as the average of the natural height of representative vegetation mounds in each quadrat. by shrubs such as S. passerina and A. brevifolia and formed mounds. The vegetation cover in each quadrat was calculated as the proportion of the ground surface covered by all living and dead plants, and the total vegetation cover was calculated as the average vegetation cover in all quadrats. This remaining dead vegetation cover was included in vegetation cover because dead plant detritus also contributes to the roughness of the ground surface and affects saltation occurrence [34,35]. In addition, we calculated community height as the average of the natural height of representative vegetation mounds in each quadrat. Figure 3. Image of a piezoelectric particle counter (ud-101) on a platform equipped with wind vanes. The instrument rotated depending on the wind direction so that it always faced upwind.

Estimating the Spatial Arrangement of Shrubs with a Semivariogram
We used a semivariogram for estimating the spatial arrangement of shrubs. Each site was divided into 4 blocks (each 10 m square), and we placed markers at the corners of each block. We took photographs of each block from the corners of sites and the center of the sites. The photographs were taken at about 7.5 m away from the vertex of each block, and at an angle from a height of 2.7 m. These images of each block were converted into images from above by using a marker-based projective transformation method in the OpenCV library [36]. The distributions of shrubs, litter, and dead forbs on the images were estimated by using Support Vector Machine classification in ArcGIS (ESRI Japan Inc., Tokyo, Japan, v.10.5) and by manual interpretation of transformed images taken from two directions. The estimated spatial distribution of the shrubs was overlaid by a 0.5 m grid and the cover of shrubs in each grid was calculated.
A unidirectional semivariogram was calculated for each site from the shrub cover data on a 0.5 m grid. The semivariogram ( ) is expressed by the following formula. Figure 3. Image of a piezoelectric particle counter (ud-101) on a platform equipped with wind vanes. The instrument rotated depending on the wind direction so that it always faced upwind.

Estimating the Spatial Arrangement of Shrubs with a Semivariogram
We used a semivariogram for estimating the spatial arrangement of shrubs. Each site was divided into 4 blocks (each 10 m square), and we placed markers at the corners of each block. We took photographs of each block from the corners of sites and the center of the sites. The photographs were taken at about 7.5 m away from the vertex of each block, and at an angle from a height of 2.7 m. These images of each block were converted into images from above by using a marker-based projective transformation method in the OpenCV library [36]. The distributions of shrubs, litter, and dead forbs on the images were estimated by using Support Vector Machine classification in ArcGIS (ESRI Japan Inc., Tokyo, Japan, v.10.5) and by manual interpretation of transformed images taken from two directions. The estimated spatial distribution of the shrubs was overlaid by a 0.5 m grid and the cover of shrubs in each grid was calculated.
A unidirectional semivariogram was calculated for each site from the shrub cover data on a 0.5 m grid. The semivariogram γ(x) is expressed by the following formula. where E{x} is the mean of x, and z i and z j are shrub cover of pairwise grid cells for which the lag distance is x (m). For each empirical variogram, the theoretical model was estimated with spherical, exponential, and linear models, and the range was calculated. Range is the lag distance x at which the semivariogram reaches a maximum and levels off. The average distance between shrubs D (m), which is twice the range value [21], is the length of nonvegetated ground surface between shrubs. Wind flow across a vegetated surface depends on the density of roughness elements. When intershrub distances are short, the wind can flow over the top of the vegetation canopy (skimming flow), whereas when the intershrub distances are large, the wind can enter the gaps between the plants and its momentum can act on the ground surface (wake interference flow) [10]. In this study, to estimate the average of unidirectional intershrub distance, we calculated a unidirectional semivariogram according to the vector mean wind direction per event (as described below) at each site in a range of 22.5 • and determined their ranges r uv . The average intershrub distance normalized by community height D/h (m/m) is thus given by the following equation.

Estimation of Saltation Flux
A total of six events were observed during the survey period. To estimate threshold friction velocity and aerodynamic roughness length z 0 during each event, an event period was identified by referring to the time-series variations of saltation flux, wind speed, wind direction, and temperature. In particular, we identified each event period such that the wind direction would be similar throughout each individual event. In the following analysis, each event at each site (site event) is treated as one dataset, and the following calculations were carried out for 20 site events.
For each site event, we calculated the 60 s average of the number of saltation particles. The saltation flux q PC was estimated as follows [37]: where d 50 is the median particle diameter of sand, ρ s is the density of sand (= 2.5 g/cm 3 , Kimura and Shinoda (2010) [13]), d PS is the sensor diameter (=0.012 m), and t 0 is the measurement time. In this calculation, we assumed, referring to Abulaiti et al. (2013) [38], that the median sand particle diameter was 0.108 mm. The saltation flux at each site was calculated as the average of flux values measured by the three PCs installed at that site.

Calculations of z 0 and d 0
The aerodynamic roughness length z 0 and the zero-plane displacement d 0 are calculated from the 10 min average of wind speed by using the following log profile equation.
where u, u * , and k are wind speed m s −1 at a height z (m), wind friction velocity m s −1 , and the von Karman constant (=0.4), respectively. A subscript of 0 indicates a surface value, at roughness length z 0 (m) and zero-plane displacement d 0 (m). Ψ M is the universal profile function of the stability parameter ζ. If the atmosphere is neutral, Ψ M is zero.
To determine d 0 , the calculations were performed iteratively by using the maximum correlation method to maximize the multiple correlation coefficient R [39,40]. where is the observed wind speed, and U ij (m/s) is the wind speed calculated by the log profile equation. Subscripts i and j refer to the profile number and the three heights that the wind was measured (1.1, 1.6, and 2.1 m, numbered 1-3 from the ground up), respectively. The value R has a maximum when the values of z 0 , d 0 , and u * i are reasonably determined. The values z 0 and u * i are temporary values determined for a given value of d 0 by the method of least squares.
To ensure neutral conditions, the data used for analysis were selected by using the following criteria, with reference to Ishizuka et al. (2009Ishizuka et al. ( , 2012 [41,42], Mikami (1997) [40], and Mikami et al. (1996) [39]: where θ j (K) is air temperature. The subscripts to u and θ represent the measurement heights, where 1 is 1.1 m, 2 is 1.6 m, and 3 is 2.1 m. After d 0 and z 0 are calculated with Equations (4) and (5) using the 10 min average data, u * per 60 s is recalculated using the log profile equation with the 60 s average of wind speed.

Calculation of u * t
The threshold friction velocity u * t is calculated using the 60 s average friction velocity data and the saltation flux. The threshold friction velocity is calculated by the following method with reference to Ishizuka et al. (2009) [41]. First, the squared difference e of saltation flux between observed data and modeled data is calculated for each 60 s: where q PC is the saltation flux observed by PCs at a reference height and q m is that obtained by modeling, using the equation of Owen (1964) [43] where c is a coefficient, ρ is the density of the atmosphere (=1293 g m −3 ), g is gravitational acceleration (=9.8 m s −2 ), and u * t is the threshold friction velocity m s −1 . E is the summation of e during one site event for a certain u * t and c: The combination of u * t and c that minimizes E is determined by successively substituting u * t from 0 to the maximum friction velocity in increments of 0.1 and c in two significant figures.
In Equation (7), the coefficient c is a fitting parameter. q PC represents the actual saltation flux, which depends on the wind strength. The coefficient c in Equation (7) is used as an index of the saltation flux in this study to allow comparisons among events with different average wind speeds. We consider the coefficient c to represent the ratio of the potential saltation flux with the vegetation elements and ground surface conditions at each site to the bare ground.

Analysis
We performed our analysis in two steps. First, we used Ward's method with cover and community height as explanatory variables to hierarchically cluster the 11 sites into groups. We then compared saltation flux, distance between shrubs D/h, vegetation cover, and community height among the groups. Cover and community height were compared among sites, and saltation flux and distance between shrubs D/h were compared among site events. Welch's t test was used for these comparisons, and the Benjamini-Hochberg method [44] was used for multiple comparisons. The Benjamini-Hochberg method corrects for multiple comparisons by controlling the false discovery rate (FDR). In this study, we used this method to estimate the q value, which represents the significance level of the FDR, from the p value estimated by Welch's t test.
Second, we conducted a regression analysis using the two indicators of saltation u * t and c. Some site event data were excluded from this analysis because of the assumption of neutral conditions. In addition, some site event data were excluded on the basis of the calculation results for z 0 and d 0 . When the multiple correlation coefficient R did not converge in the range of d 0 < (bottom weather meter height), we excluded that site from the analysis. We also excluded sites with large saltation fluxes even when d 0 > 0, because the wind speed near the ground surface, which is very low when d 0 > 0, would be expected to result in small saltation fluxes [45]. In addition, referring to Daniels (1973) [46], we excluded sites with too small z 0 values for the shrubland from the analysis. Finally, data from 8 sites, for a net of 10 site events, were used in subsequent analyses (Table 1, Figure 4). Table 1. Conditions of the soil surface and vegetation at each site. The initial letter of the Site ID (S or N) indicates that the site is located on the south or north pediment. The soil surface conditions are described qualitatively in relation to gravel and crust conditions in the whole study area. Under "Gravel size", "medium" means 1-2 cm and "small" means 0.5 to 1 cm in size. Under "Gravel amount", "large" means that about 70% of the surface is covered by gravel, and "middle" means that about 40% of the surface is covered by gravel. We also recorded the presence of crusts and cracks: a "weak" crust means a slightly aggregated surface, and "small" cracks means the crust is slightly dry and starting to crack. For cover and community height, the mean and standard error (S.E.) in the survey quadrats are shown (n = 9). We estimated three relationships by multiple linear regression. Firstly, we examined the relationship between threshold friction velocity u * t and roughness length normalized by community height z 0 /h. Secondly, we verified the relationship between the coefficient c and vegetation elements (vegetation cover and distance between shrubs (D/h)) and their interaction. Together, these analyses represent the relationship between saltation occurrence and vegetation frequency and quantity. Finally, we examined the relationship between z 0 /h and vegetation elements and their interaction. The vegetation cover was We estimated three relationships by multiple linear regression. Firstly, we examined the relationship between threshold friction velocity * and roughness length normalized by community height /ℎ. Secondly, we verified the relationship between the coefficient c and vegetation elements (vegetation cover and distance between shrubs ( /ℎ)) and their

Wind Conditions and Dust Events
During the entire study period, six saltation events were observed. For each event, we compared the vector mean of wind across sites. In events 1, 2, 4, and 5, western to southwestern winds dominated (mean wind direction (S.E.): event 1, 254.8 (4.7); event 2,

Wind Conditions and Dust Events
During the entire study period, six saltation events were observed. For each event, we compared the vector mean of wind across sites. In events 1, 2, 4, and 5, western to southwestern winds dominated (mean wind direction (S.E.): event 1, 254.8 (4.7); event 2, 250.3 (7.6); event 4, 235.9 (3.8); event 5, 267.8 (3.5)). In events 3 and 6, northern to northeastern winds dominated (10.9 (10.0) and 47.1 (7.2), respectively). Wind speeds were highest in event 1 (9.3 m/s (0.1)) and low in events 3 and 5 (4.4 m/s (0.1) and 4.5 m/s (0.1), respectively). For the site events used in the regression analyses, the wind directions during the same event were approximately the same (Figure 4). Event 5 at any site was not used in the regression analysis.

Ground Surface Conditions
Although we conducted no quantitative gravel or crust surveys during this study, soil surface conditions at all sites were very similar. At all sites, the ground surface was covered by a weakly aggregated soil crust (Table 1). At most sites, approximately 70% of the ground surface was also covered by gravels 0.5-2 cm in size. Soil moisture was measured in each quadrat at each site once during the study period. The volumetric water content was almost the same within and among groups (mean soil moisture (S.E.): group 1, 0.083% (0.003); group 2, 0.082% (0.002); group 3, 0.081% (0.007).

Linear Regression Analysis
u * t was well explained by z 0 /h.; the derived regression equation was u * t = 0.82 + 0.08 ln(z 0 /h) (R 2 = 0.87, p < 0.001) ((a) in Table 2, Figure 6a). In contrast, u * t did not show a clear trend in relation to vegetation cover or community height. The multiple regression analysis results for coefficient c against cover and D/h, however, showed that cover and the interaction term (cover * D/h) had significant effects on coefficient c ((b) in Table 2). The regression equation was c × 10 4 = 1.98 − 0.167 + 0.00548 D/h * cover − 0.0216D/h (R 2 = 0.56, p = 0.049; D/h and cover were centered at their respective means of 73.1 and 23.6). Coefficient c tended to decrease as cover increased, and the rate of decrease depended on D/h. When D/h was large, coefficient c decreased drastically as cover increased, whereas when D/h was small, coefficient c did not change significantly as the cover changed ( Figure 6).  The multiple regression analysis between /ℎ and vegetation factors (cover and intershrub distance) showed that /ℎ was significantly influenced by vegetation cover and the interaction term (cover * /ℎ ) ((c) in Table 2). The regression equation was ℎ ⁄ × 10 = 4.79 − 0.0603 ℎ ⁄ − (0.298 + 0.0104 ℎ ⁄ ) * cover ( = 0.63, = 0.029 , with /ℎ and vegetation cover centered at their respective means of 73.1 and 23.6). /ℎ decreased as cover increased, especially when /ℎ was large (Figure 7). In contrast, when /ℎ was small, the change in /ℎ as cover increased was small. The multiple regression analysis between z 0 /h and vegetation factors (cover and intershrub distance) showed that z 0 /h was significantly influenced by vegetation cover and the interaction term (cover * D/h) ((c) in Table 2). The regression equation was z 0 /h × 10 2 = 4.79 − 0.0603 D/h − 0.298 + 0.0104D/h * cover (R 2 = 0.63, p = 0.029, with D/h and vegetation cover centered at their respective means of 73.1 and 23.6). z 0 /h decreased as cover increased, especially when D/h was large (Figure 7). In contrast, when D/h was small, the change in z 0 /h as cover increased was small. and the interaction term (cover * /ℎ ) ((c) in Table 2). The regression equation was ℎ ⁄ × 10 = 4.79 − 0.0603 ℎ ⁄ − (0.298 + 0.0104 ℎ ⁄ ) * cover ( = 0.63, = 0.029 , with /ℎ and vegetation cover centered at their respective means of 73.1 and 23.6). /ℎ decreased as cover increased, especially when /ℎ was large (Figure 7). In contrast, when /ℎ was small, the change in /ℎ as cover increased was small.

Discussion
The threshold friction velocity is affected by various vegetation factors. The roughness length z 0 /h varies depending on the distribution of roughness elements on the ground surface. King et al. (2006) [47] pointed out that the roughness length varies with wind direction when the vegetation arrangement is heterogeneous. Although the wind direction was different for each event (Figure 4), we estimated the roughness length for each site event and succeeded in capturing the changes in roughness length due to different wind directions under heterogeneous vegetation. In this study, the intershrub distance D/h was not linearly related to z 0 /h, but the relationship between z 0 /h and vegetation cover varied with D/h; this result suggests that roughness length reflects the spatial arrangement of vegetation. The threshold friction velocity u * t was linearly related to z 0 /h, which expresses the spatial distribution of vegetation. Wind tunnel experiments have also shown that the threshold friction velocity is explained by roughness length [17], and our field observations confirmed this experimental result. According to Shao (2008) [48], the actual threshold friction velocity u * t at a site can be determined by multiplying the ideal threshold friction velocity on dry bare ground u * t0 by correction functions that take account of various ground surface conditions, including vegetation.
where f λ (λ), f θ (θ s ), and f cr (cr) are the correction functions for lateral cover λ due to vegetation and gravel, soil moisture θ s , and crust cr. This study showed that threshold friction velocity is proportional to the natural logarithm of roughness length z 0 /h (Section 3.4). By converting z 0 /h to λ with Equation (10) [18], we can calculate f λ (λ) by Equation (11) [48] z 0 h = 0.96 λ 1.07 (10) f λ (λ) = (1 − m r σ r λ) 1 2 (1 + m r β r λ) where m r is a correction parameter, σ r is the ratio of basal area to frontal area of roughness elements, and β r is the ratio of the wake coefficient on bare ground to that on a vegetated ground surface. We assumed m r = 0.5, σ r = 1.0, and β r = 200, values based on those used in WRF-chem [49], a widely used numerical weather prediction model. As a result, we could convert the regression equation obtained in this study (Section 3.4) to Equation (12).
At the study sites, the amounts of crust and gravel and gravel size were approximately the same among sites (Table 1). In addition, the ground surface was very dry (Section 3.3; [50]), and no rain fell during the survey period. If the effect of soil moisture is ignored, this result satisfies Equation (9). The threshold friction velocity u * t0 was 0.34 for the surface conditions of the study site without consideration of the influence of vegetation.
The value of f λ (λ) was greater than one, and the threshold friction velocity increased as roughness increased. Vegetation acting as surface roughness elements extracts momentum from the wind acting on the ground surface [26,51]. Increasing roughness causes the threshold friction velocity to increase and the drag of the ground surface to decrease.
The value of the roughness length standardized by the community height z 0 /h was in the range of 10 −3 − 10 −1 ; this value is reasonable when compared to the results of previous studies (e.g., [18,52]) in which the relationship between roughness length and vegetation cover was examined for homogeneous vegetation distributions. Vegetation, as a roughness factor on the ground surface, influences wind speed and wind momentum [26,51]. According to Wolfe and Nickling (1993) [10] three flow regimes describe the relationship between vegetation and boundary layer wind flow: (1) when the cover is small (<16%), each individual roughness element affects the wind flow separately (isolated roughness flow). In this regime, an increase in vegetation cover causes the roughness length z 0 to increase [3,52]; (2) when the vegetation cover is moderate (16-40%), the wakes of individual roughness elements interfere with each other (wake interference flow); (3) at high coverage (>40%), skimming wind flow occurs across the top surface of the roughness elements. In wake interference flow or skimming flow, roughness reaches a maximum value when the vegetation cover reaches a certain threshold; then, as the vegetation cover increases further, the roughness surface becomes aerodynamically smooth and roughness decreases and then levels off (e.g., [3,18,52]). The threshold cover of these three regimes and the cover where the roughness reaches a maximum can vary depending on the vegetation structure. For example, Liu et al. (2021) [52] reported that the threshold cover in the three regimes varies depending on plant shape. When the roughness elements are cylindrical, wind flow shifted to wake interference flow at 6.6% cover, and 21.1% cover causes the wakes to completely overlap and the regime to shift to skimming flow. Our results suggest that the threshold cover of these three regimes can also vary depending on the spatial distribution of the shrubs. A large intershrub distance indicates that the shrubs constitute large patches. Large gaps separate individual shrubs or patches of shrubs. As vegetation cover increases, additional shrubs fill in these gaps. When the shrub cover exceeds a certain level, the wakes of the individual roughness elements (shrubs or patches) begin to interfere with each other, the regime shifts to wake interference flow, and the roughness length decreases. In contrast, when the intershrub distance is small, individual shrubs are isolated or distributed in small patches; this distribution pattern is relatively homogeneous. In this study, because there were only a few sites with small intershrub distances and moderate vegetation cover, additional research is needed. However, when the vegetation cover was small, the roughness length increased as the vegetation cover increased, and when the vegetation cover was large, the roughness length decreased under a skimming flow regime such that it was similar to that when the vegetation cover was small.
In this study, the sites were selected so that the community height would be similar among sites but the vegetation cover and intershrub distance would be different. Many previous studies have shown that the saltation flux is small above a threshold vegetation cover [12][13][14][15], and our results support those results (Section 3.1). Coefficient c decreased with increased vegetation cover, indicating that an increase in vegetation cover reduces the occurrence of sediment transport. However, the multiple regression analysis results suggest that sediment transport was affected by the intershrub distance (Section 3.4). When the intershrub distance was large, the saltation flux increased rapidly as the vegetation cover decreased. In contrast, when the intershrub distance was small, which indicates a homogeneous distribution of shrubs, the saltation flux hardly increased even when the vegetation cover decreased. The reason for this result is that a homogeneous distribution of shrubs reduces the amount of continuous bare ground between shrubs, where wind momentum can act on the ground surface and cause particle saltation. In Okin's model (2008) [22], wind speed is sharply reduced downwind of shrubs and recovers with distance from the shrubs. Total saltation flux in a site is summation of saltation flux at the microsite where friction velocity exceeds the threshold [22]. Because wind speed is lower in the area sheltered by shrubs, the wind speed does not exceed the threshold wind speed of bare ground conditions at the microsite, and thus the occurrence of saltation is suppressed. When the intershrub distance is large, the downwind areas sheltered by the plants are isolated. Therefore, an increase in vegetation cover directly increases the sheltered area, and the saltation flux significantly decreases as vegetation cover increases. In contrast, if the intershrub distance is small, then the shrubs are close to other shrubs of the upwind direction. Therefore, areas protected by shrubs may overlap. Thus, even when the vegetation cover is low, the saltation flux is relatively small because the measurable area is already protected, and the saltation flux does not decrease significantly even when the vegetation cover increases because of the overlap of the protected areas.
Among the sites selected in this study, the intershrub distance varied by site and azimuth. In between-group comparisons, intershrub distance D/h was greatest and variability was relatively high in group 2, which had moderate vegetation cover (Section 3.1). Therefore, this result suggests that vegetation canopy gaps may be largest, depending on the vegetation distribution pattern, when the vegetation cover is moderate. Many previous studies have reported that the spatial distribution of shrubs in arid or semiarid lands changes from random or aggregated to regular as the shrubs grow [53][54][55][56]. As shrubs grow, their distribution shifts from the random distribution caused by seed dispersal to a patchy distribution that reflects the spatially heterogeneous distribution of soil nutrients and moisture. In the patchy distribution, although the intershrub distance depends on the distribution patterns of nutrients and moisture (Figure 8a,b), as vegetative cover and patch size increase because of shrub growth and litter trapping, the patches become more likely to encroach on the intervening bare ground, thereby causing the intershrub distance to decrease (Figure 8b,c). Although the vegetation coverage in group 2 exceeded the threshold cover of 12-15% for sediment transport in field experiments, wake interference flow dominates when the distance between shrubs is large, and drag transfers to the ground surface in gaps between vegetation may cause sediment transport to be large [3,23]. Therefore, it is important to consider not only vegetation cover thresholds but also the spatial distribution of vegetation in rangeland management.
The limitations of this study include the limited sample size and the lack of high accuracy in estimating the spatial distribution of shrubs. In particular, there were few sites with large vegetation cover and large intershrub distance, and there were also few sites with small intershrub distance and moderate vegetation cover. Therefore, a non-linear relationship between z 0 /h and vegetation cover might not be detected by our method of analysis. However, when the vegetation cover is large, the intershrub distance is expected to be small, because the bare ground surface area is small, and arid land vegetation is characterized by a heterogeneous distribution. As a result, it is difficult to find many such sites in the field. The use of unmanned vehicles (drones) to obtain aerial images for determining the spatial distribution of shrubs has recently been proposed [3,56,57]. At our study sites, the shrub-steppe community was dominated by S. passerina and A. brevifolia, but different vegetation communities occur in a bottom part of the valley that is slightly higher than the other bottom part (yellow region in Figure 1) where Nitraria sibirica Pall. is dominant. Compared with N. sibirica, S. passerina and A. brevifolia are relatively small shrubs, and they tend to be relatively homogenously distributed when the vegetation cover is similar. In contrast, the relatively larger N. sibirica shrubs form mounds (so-called fertility islands); [58,59] that tend to be distributed heterogeneously. The spatial heterogeneity of shrubs affected the occurrence of saltation even in relatively homogeneous S. passerina and A. brevifolia communities; thus, the spatial distribution of shrubs may have a greater effect on the occurrence of saltation in communities of larger shrubs such as N. sibirica. It is also possible that plants in vegetation communities with different species compositions may lead to different roughness factor values because of differences in their porosity and flexibility. Thus, further investigation of communities dominated by other species is required. the vegetation distribution pattern, when the vegetation cover is moderate. Many previous studies have reported that the spatial distribution of shrubs in arid or semiarid lands changes from random or aggregated to regular as the shrubs grow [53][54][55][56]. As shrubs grow, their distribution shifts from the random distribution caused by seed dispersal to a patchy distribution that reflects the spatially heterogeneous distribution of soil nutrients and moisture. In the patchy distribution, although the intershrub distance depends on the distribution patterns of nutrients and moisture (Figure 8a,b), as vegetative cover and patch size increase because of shrub growth and litter trapping, the patches become more likely to encroach on the intervening bare ground, thereby causing the intershrub distance to decrease (Figure 8b,c). Although the vegetation coverage in group 2 exceeded the threshold cover of 12-15% for sediment transport in field experiments, wake interference flow dominates when the distance between shrubs is large, and drag transfers to the ground surface in gaps between vegetation may cause sediment transport to be large [3,23]. Therefore, it is important to consider not only vegetation cover thresholds but also the spatial distribution of vegetation in rangeland management. Figure 8. Schematic diagram of changes in intershrub distance (double-headed arrows) with increasing vegetation cover. Shrubs are randomly distributed in a landscape with small shrubs and low vegetation cover (a), but their distribution becomes patchy (b) as the shrubs grow. When the small cover regime (a) shifts to a moderate cover regime (b), only some of the shrubs remain because of the spatial heterogeneity of soil nutrients and intra-and interspecific competition among shrubs. As the shrubs continue to grow, the patch size increases (c) such that the patches begin to encroach on the intervening bare ground, with the result that the intershrub distance decreases.
The limitations of this study include the limited sample size and the lack of high accuracy in estimating the spatial distribution of shrubs. In particular, there were few sites with large vegetation cover and large intershrub distance, and there were also few sites with small intershrub distance and moderate vegetation cover. Therefore, a non-linear relationship between /ℎ and vegetation cover might not be detected by our method of analysis. However, when the vegetation cover is large, the intershrub distance is expected to be small, because the bare ground surface area is small, and arid land vegetation is characterized by a heterogeneous distribution. As a result, it is difficult to find many such sites in the field. The use of unmanned vehicles (drones) to obtain aerial images for determining the spatial distribution of shrubs has recently been proposed [3,56,57]. At our study sites, the shrub-steppe community was dominated by S. passerina and A. brevifolia, but different vegetation communities occur in a bottom part of the valley that is slightly higher than the other bottom part (yellow region in Figure 1) where Nitraria sibirica Pall. is dominant. Compared with N. sibirica, S. passerina and A. brevifolia are relatively small Figure 8. Schematic diagram of changes in intershrub distance (double-headed arrows) with increasing vegetation cover. Shrubs are randomly distributed in a landscape with small shrubs and low vegetation cover (a), but their distribution becomes patchy (b) as the shrubs grow. When the small cover regime (a) shifts to a moderate cover regime (b), only some of the shrubs remain because of the spatial heterogeneity of soil nutrients and intra-and interspecific competition among shrubs. As the shrubs continue to grow, the patch size increases (c) such that the patches begin to encroach on the intervening bare ground, with the result that the intershrub distance decreases.