Spatial vegetation structure and its effect on wind erosion of Alxa dryland ecosystem

In wind erosion models, previous parameters related to vegetation morphology and density are limited in describing the spatial distribution of vegetation that influences surface heterogeneity. Thus, it is not fully understood how spatial vegetation patterns affect wind erosion on a field-scale. Based on an investigation of 36 plots of vegetation in Alxa Plateau, northwestern China, we established a multivariate linear model for temporally and spatially averaged aerodynamic roughness length (Z 0) incorporating the height, roughness density, regularity of vegetation patches (curvature) and spacing between patches (connectivity). The curvature positively interacted with the connectivity in affecting the mean Z 0, while it was the most important factor affecting the standard deviation of Z 0. The connectivity modulated the roughness density in affecting the standard deviation of Z 0. The spatial-related terms contributed 37% and 62% to the model variance of the mean and standard deviation of Z 0, respectively. Our results validate the importance of spatial vegetation structure in the vegetation-airflow interactions, with a suggestion of estimating the heterogeneity of surface erodibility by intuitive spatial parameters. Based on that spatial vegetation patterns reflect the ecosystem states, a strengthened linkage between wind erosion and vegetation stability may be useful in erosion regulation in drylands.


Introduction
In arid and semi-arid ecosystems, wind erosion is a critical component of land degradation and desertification (Kok et al 2018, McLaughlin 2019. Vegetation in these areas plays an important role in reducing wind erosion by providing a sheltered area, reducing wind momentum, and impeding sand transport (Wolfe and Nickling 1993). Due to water scarcity and variability (Kéfi et al 2007, Wang et al 2017, natural rain-fed vegetation in drylands often exhibits spatial patterns with discontinuous patches surrounded by bare soil, such as spots, bands, labyrinth, and other irregular patterns (Aguiar andSala 1999, Maestre et al 2021). As a result, land surface erodibility is strongly influenced by the complex interplay between airflow and characteristics of vegetation patches (Lee and Soliman 1977, Mayaud et al 2016, Trautz et al 2017.
The ability to quantify surface roughness with spatial vegetation structure can provide guidelines for wind erosion management practices.
Airflow through sparse vegetation is often characterized by aerodynamic roughness length (Z 0 ) in meteorology (Raupach 1992). Based on the 'law of the wall' , horizontal wind velocity decreases logarithmically with decreasing height above the surface (Bagnold 1974). As a length scale, Z 0 is the theoretical height at which the horizontal wind velocity equals zero (Wolfe and Nickling 1993). Using Z 0 to calculate the shear velocity, which represents the surface shear stress exerted by the airflow, and comparing the shear velocity to the wind velocity threshold to initiate mass transport is the basis for informing wind erosion (King et al 2005, Foroutan et al 2017. While Z 0 was commonly treated as a constant for simplification in many wind erosion models (Okin 2008, Munson Seth Figure 1. Aerodynamic roughness length on the rough surfaces including (a) the wind velocity recovery behind an isolated roughness element, and (b) as the spacing between roughness elements decreased, the sheltered areas overlapped and the wind velocity recovery is interfered, aerodynamic roughness length increases as a result. Dashed lines indicate the wind velocity recover in the downwind. Shaded areas are the sheltered areas provided by roughness elements and are related to morphological characteristics. Adapted from Webb et al (2021). CC BY 4.0.

Figure 2.
The telephone pole problem. Suppose the lateral area is constant in the three scenarios, (a) a single telephone pole, (b) it is cut uniformly into four shorter poles, and (c) it is cut and reshaped into short poles of different sizes. In the three scenarios, while the model estimations based on the lateral area are the same, the wind erosion changes in physical reality. Adapted from Mayaud and Webb (2017). CC BY 4.0. et al 2011, in practice it varies across heterogeneous surfaces that are mainly affected by roughness elements (figure 1).
An alternative method is to estimate Z 0 by the morphological characteristics of roughness elements related to their density (Marticorena et al 1997, Pi et al 2020, called lateral cover (λ), front area index, or roughness density (Lettau 1969): where N is the number of roughness elements, H and B are the height and breadth of each roughness element, respectively (figure 2(a)), and S is the observed area. However, models that include only morphological parameters (height and lateral cover) have difficulty fitting the variability of Z 0 in the natural field with inhomogeneous vegetation patches (King et al 2005). Because λ characterizes the amount of roughness elements per unit area (Raupach et al 1993), but does not quantify their spatial distribution. The weakness of estimating wind erosion without the spatial structure of roughness elements is highlighted by the telephone pole problem (Okin 2008). For example, cutting a single telephone pole into four shorter poles without changing the lateral cover leaves the estimation of horizontal sand fluxes unchanged (figures 2(a) and (b)), which does not correspond to the physical reality that the horizontal fluxes are decreased.
To complement the deficiency in estimating wind erosion with λ, the connectivity of canopy gaps or bare soil was developed to depict spatial vegetation structure by characterizing the spacing between vegetation (McGlynn and Okin 2006, Zhang et al 2021. Nevertheless, it is unable to quantify the variability of the vegetation patch sizes. Also, assuming that the telephone pole is plastic and can be cut into shorter poles of different sizes (figure 2(c)), as opposed to uniform poles, this may lead to the fluctuations in Z 0 as demonstrated in previous studies (King et al 2006, Kono andOkuro 2021). However, a quantitative description of the regularity for vegetation patches (regular or irregular) is lacking.
The patch size distributions of vegetation have been found to follow power-laws in a wide range of drylands, indicating irregular vegetation patterns with many small patches and rare large patches (Kéfi et al 2007, Scanlon et al 2007. It has been demonstrated recently that the patch size distributions deviate from power-laws when the patterns changed from irregular to regular (Berdugo et al 2019). A parameter called curvature is used to depict the shape of patch size distributions and determine the fitness of the distribution to power-laws (Schneider andKéfi 2016, Berdugo et al 2019). While the shape of patch size distributions is believed to reflect the interplay between vegetation and its environment (Kéfi et al 2011, this parameter has not been applied to simulate abiotic processes that related to the spatial distribution of vegetation. Here, based on the investigation of different spatial structure in desert vegetation in Alxa Plateau, northwestern China, we introduced the curvature to describe the regularity of the vegetation patches, and the connectivity to characterize the spacing between patches. We performed a multi-model inference procedure to examine the standardized effects of wind speed, mean plant height, roughness density and spatial structure of vegetation on the temporally and spatially measured Z 0 on field-scale. We hypothesized that (a) the rough surfaces with heterogeneously distributed vegetation are better described by the combination of both the parameters characterizing bare soil and vegetation patches; (b) these two spatial parameters are effective in estimating Z 0 by interacting with each other or modulating morphological characteristics and roughness density, thus improving the evaluation of surface erodibility.

Study area
This study was conducted along the southeastern margin of the Tengger desert (figure S1) in the southeastern of Alxa Plateau (104 • 12 ′ -105 • 30 ′ E, 37 • 21 ′ -38 • 18 ′ N). The area has a typical continental climate. The average annual temperature is 8.3 • C, with a potential evaporation range from 2800 to 3000 mm. The annual precipitation falls between 150 and 200 mm, of which 80% is concentrated between June and September (Chen 2010, Ma andWang 2020). The average annual wind speed is 3.4-4.7 m s −1 . The prevailing wind direction is northnorthwest. Dust and sand storm occurs centrally between March and May, and days with wind greater than 17 m s −1 were approximately 10 d per year (Li et al 2022). Basic soils mainly are aeolian sandy soils, light brown soils and gray desert soils. The vegetation type is temperate desert, with average coverage approximately 18%. Shrubs are the dominant life form, with Reaumuria soongorica, Salsola passerine, Ammopiptanthus mongolicus and Nitraria tangutorum being the predominant species.

Field measurement
We investigated the spatial patterns of vegetation in 36 plots (size: 30 m × 30 m) at relatively flat topography. The plots were more than 5 km apart from each other. Orthophotos with a resolution of 0.3 cm for each plot were snapped by a UAV (DJI Phantom 4 equipped with a 4K camera) at the height of 15 m (Prošek andŠímová 2019). Then mosaic images with a resolution of less than 0.1 cm were obtained by the Agisoft Photoscan (version 1.3). To characterize the roughness elements in vertical space, we measured the heights of all alive shrubs and the remains of dead ones. The plant height of perennial grasses identifiable on the images, such as Allium mongolicum and Cynanchum mongolicum, were also measured. We summed the height of sand mounds and plants as the height of nabkhas. For each plot we collected nine soil samples randomly on the intershrub bare ground to determine their particle size distribution by a laser granularity instrument (Malvern MS 2000, Gao et al 2014). Surface wind profiles were observed during strong wind events by temporally and spatially averaged measurements (Zobeck et al 2003). We measured wind velocity at 3 m, 2 m, 1 m and 0.3 m above the ground by a portable micro weather station with four anemometers (Onset S-WSA-M003) at 5-10 randomly selected intershrub positions (five positions for regular and highly uniform vegetation and more for complicated situations) within each plot. Wind velocity was measured at 1 s intervals for 30 min, and a 10 s mean was recorded. We first measured wind velocity profiles on a bare ground in our study area to correct the equipment.

Data analysis
We calculated the average height of all vegetation patches in each plot. The mosaic images were first processed by color filtering via the eCognition (version 9.3) to distinguish vegetation patches from bare soil. Then the vegetation layers were imported to Arc-GIS (version 10.3) and modified visually. The areas of vegetation patches were calculated, and the breadth of each patch was estimated as the diameter of a circle. The lateral cover (λ) was calculated according to equation (1), and the fractional cover of vegetation was also calculated.
We used the shape of patch size distribution to characterize whether the spatial patterns of vegetation are regular or irregular (Berdugo et al 2019). The patch size distribution of each plot was described as the cumulative frequency of vegetation patch sizes (X), i.e. the areas of vegetation patches, greater than or equal to a value (x), P(X ⩾ x). We fitted a quadratic regression for each patch size distribution (Schneider andKéfi 2016, Berdugo et al 2019): where the coefficient of the quadratic item (a) was used to characterize the degree of curvature of the regression. For an irregular pattern, the patch size distribution follows a power-law distribution, which shows a straight line in the logarithmically transformed axes, and the coefficient is close to 0. As the coefficient deviates from 0, the degree of curvature increased, and the patch size distribution follows a lognormal or exponential distribution, which exhibits a more regular pattern (Berdugo et al 2019). The curvature ranges from −1 to 1, we took the absolute value of the coefficients because regression curves are generally convex, which fit negative coefficients. The connectivity of canopy gaps was obtained by converting the classified images of vegetation patches and bare soil to zero-one matrices (1 for canopy gaps and 0 for vegetation). The images were previously coarsened to a resolution of 0.1 m because we considered this to be the minimum size of canopy gaps which have no potential for wind erosion (Herrick et al 2005). Each row of the matrices was set as a transect to measure canopy gap sizes. The connectivity of one direction can be calculated as (McGlynn and Okin 2006): where⃗ r is a lag vector with gap size |r|, n is the number of successive sets of gaps along⃗ r in a matrix, and I i is a class indicator, the value equals 1 for canopy gaps and 0 for vegetation. As |r| increases, the connectivity decays according to an exponential function: where C 0 is the connectivity for |r| = 0, which equals to the fractional cover of canopy gaps, and α may be interpreted as the range of traditional variogram with a value equal to |r| when the connectivity is C 0 /e. Therefore, α is the gap size where the connectivity become stable. Based on the direction of ⃗ r changes from 0 • to 360 • with a 2 • intervals (Zhang et al 2021), we obtained the connectivity of all directions and calculated the average α to characterize the connectivity of canopy gaps for each plot (McGlynn and Okin 2006). To estimate Z 0 , the wind profiles were chosen as the wind speed at 3 m above the ground greater than 6 m s −1 for the purpose of approximate neutral air conditions. The equation is as following (Jackson 1981): where U is the wind velocity at height Z, and u * is the shear velocity, k is the Karman constant, d is the zero-plane displacement, which refers to the upward displacement of the effective surface on a rough surface. This parameter was estimated as a correction factor because the height of the weather station was determined from the actual surface. To determine the value of Z 0 and d, we performed iterations to maximizing correlation coefficient of the least square correlation for U and ln (Z − d) (Mikami et al 1996, Liu et al 2021. Then Z 0 equals to the anti-log of the quotient of the intercept and slope. We calculated the mean and standard deviation of Z 0 with all the positions measured in each plot, and the wind speed greater than 6 m s −1 at 3 m above the ground was also averaged. We developed a multiple regression model (García-Palacios et al 2018) to explore the factors affecting the mean and standard deviation of Z 0 in the natural environment, including wind velocity, weighted mean soil particle size, the height, fractional cover, lateral cover, and spatial structure of vegetation. We first considered the correlations between these factors and excluded those with strong correlations (coefficients greater than 0.5) to avoid multicollinearity in the regression model. The quadratic items for all the factors were also included, and the two-way interactions of wind velocity, mean height, lateral cover, and spatial structure were involved. The best model was selected based on the minimized corrected Akaike information criterion (AICc) by the package MuMIn (Bartoń 2022). After that an averaged model was obtained by the models with ∆AICc < 2 from the optimal model (table S1). To compare the relative strengths of the factors, the data were standardized before regression, and the contribution of each factor to the explanatory power of the models was represented by the ratio of its corresponding estimated coefficient to the sum of all the coefficients . Response variables were log-transformed to meet the assumptions of normality.
Based on the results of the parameter estimation, the interaction terms affecting the mean and standard deviation of Z 0 were illustrated through the predictions of the mean model, while the other factors in the model were fixed to the mean value, which was exactly zero because the data were normalized (Le Bagousse-Pinguet et al 2017). All the data analyses were performed in R (version 4.1.1, R Core Team 2021), except the calculation of connectivity was coded by Matlab (version R2020b).

Relationships among environmental and vegetation indices
Of all the environmental and vegetation indices, wind velocity at 3 m above the ground was only weakly correlated to the mean particle size and not with vegetation characteristics and structure (figure 3). The height and breadth of vegetation patches were strongly correlated. Apparently, the lateral cover was positively correlated to the height and fractional cover of vegetation. The mean size of soil particles was strongly correlated to the spatial structure of vegetation. Therefore, we excluded the mean particle size, breadth and fractional cover of vegetation in the multiple linear models.
To further understand the spatial structure parameters, three examples were shown in figure 4. Intuitively, the plot b (figure 4(Ab)) showed a spots pattern and was more regular than the others, with the ratio between the maximal and minimal patches smaller than 10 3 (figure 4(B)), and the curvature of the patch size distribution was 0.96. The patch size distributions of the other two plots (figures 4(Aa, b)) showed approximately straight lines on the log-transformed axes, which were better fitted to the power-laws, with the curvature were 0.02 and 0.07, respectively. Although the shape of the patch size distributions was quite different, the connectivity averaged from all the directions for plot a and b were similar (figure 4(C)), with a value of 2.51 m and 3.11 m, respectively. In our study area, correlation between the curvature and the connectivity was negative ( figure 4(D)). Large patches of bare soil were visible in the plot when the connectivity was large (9.62 m) in irregular patterns (figure 4(Ac)).

Effect of environmental and vegetation parameters on the aerodynamic roughness length
The explanatory power of our model was similar for the mean and standard deviation of Z 0 (figure 5), with R 2 s were 0.62 and 0.64, respectively. The wind velocity negatively affected the mean whereas the average height of patches had positive effects on both the mean and standard deviation (P < 0.05). The linear regression model did not show a significant relationship between the lateral cover and the mean of Z 0 . However, it was positively correlated to the standard deviation, and interacted with the connectivity. The curvature had a negative effect on the standard deviation, indicating that the variation of Z 0 was smaller in regular patterns. The interaction of the two spatial parameters significantly affected the mean and accounted for 18.75% of the model variance. For the mean Z 0 model, structure-related terms contributed 37.04%, while the terms related to height and roughness density contributed 45.04% of the variance. In addition, structure-related terms explained up to 62.30% of the model variance for standard deviation, indicating that the spatial structure of vegetation played an important role in estimating Z 0 .

Predictions of the interactions with spatial structure parameters
The regularity of vegetation patches, the connectivity of canopy gaps, and roughness density were interdependent in influencing Z 0 . As predicted by the interaction between the structure parameters, the mean Z 0 of regular patterns was smaller than that of irregular patterns when gaps were weakly connected ( figure 6(A)). In contrast, the regular pattern was related to the higher mean value of Z 0 when the connectivity was much larger ( figure 6(A)). With small connectivity, the standard deviation of Z 0 varied gradually as the roughness density increased ( figure 6(B)). When the connectivity was large, the standard deviation increased rapidly as roughness density increased.

Discussions
Here we described the spatial structure of heterogeneously distributed vegetation in the Alxa ecosystem by characterizing whether the patch patterns were regular or not, combined with the spacing between vegetation patches. By adding the spatial structure parameters in linear regression multi-variable models, we showed that compared to the morphological characteristic and roughness density that used to estimate Z 0 , the interaction of the connectivity and curvature were positively correlated to the mean Z 0 at the field scale, while the curvature was mainly correlated to the standard deviation of Z 0 .
The ability to characterize vegetation heterogeneity is vital to evaluate mass transport in drylands, such as water run-off and wind erosion (Trautz et al 2017, Rodríguez et al 2018. With the development of aerial mapping technology, it is a common idea to extract the spatial information of vegetation from digital images (Scanlon et al 2007, Xu et al 2015. Aiming at the limitation of roughness density to describe the spatial distribution of vegetation in previous wind erosion models, McGlynn and Okin   (2006) quantified vegetation spacing using the connectivity of canopy gaps based on image processing. It is convenient and intuitive to use the connectivity to characterize the connecting pathways in bare soil (figure 4), where resources could be easily transported (Okin andGillette 2001, Okin et al 2015). However, in our investigation, the regularity of vegetation patches was contrastingly different for canopy gaps with similar connectivity. Vegetation patches of various sizes not only affect the development of connecting pathways, but also provide different sheltered areas for the surface, thereby influencing downwind wind velocity recovery (Wolfe andNickling 1993, King et al 2006). Thus, we added the curvature of the patch size distributions as a compliment to the connectivity to further illustrate the spatial variability of wind erosion.
In our averaged models, λ had no effect on the mean Z 0 , whereas the interaction between the spatial parameters had a significant effect on the mean Z 0 . Due to the harsh environment (Maestre et al 2021), the diverse morphology of different shrub species in our study area, such as the evergreen A. mongolicus, Caragana korshinskii with an inverted cone canopy shape, and cushion Oxytropis aciphylla, may influence the calculation of λ and thus undermine its effectiveness in estimating Z 0 in the field (Liu et al 2021). Based on the uniform or staggered distribution of roughness elements (Wooding et al 1973), previous models have pointed out that the spacing of roughness elements has an important effect on Z 0 (Wolfe and Nickling 1993). In our model, the curvature modulated the effect of the connectivity on Z 0 . Irregular patterns showed higher Z 0 when the connectivity was small, probably because large patches interfere more with wind speed recovery (King et al 2006). In contrast, when the connectivity was large, the mean Z 0 decreased as vegetation patterns shifted from regular to irregular. In our investigation, large patches of bare soil emerged in irregular patterns with large connectivity, and the difference in Z 0 between bare soil and vegetated surface reached a magnitude of 10 4 (Wieringa 1992), resulting in a significant decrease in the mean value. The positive effect of the interaction resulted in the largest Z 0 in regular patterns with the largest connectivity. Although this is consistent with the effect of airflow regimes on Z 0 (Mayaud et al 2016), in the absence of observational evidence, we expect to validate the reliability of model extrapolation in future wind tunnel experiments.
While acknowledging that in reality Z 0 varies across heterogeneous surfaces, it is often treated as a constant in models that evaluate wind erosion in natural fields (Munson Seth et al 2011, Foroutan et al 2017. Our results showed that the standard deviation of Z 0 was significantly affected by the characteristics of vegetation. In accordance with the telephone pole problem (Okin 2008, Webb et al 2021, the connectivity modulated λ in affecting the standard deviation of Z 0 . When the connectivity increased, the increase of λ mainly take place in vertical space, and the standard deviation of Z 0 increased as a result. Moreover, the variability in standard deviation of Z 0 was strongly affected by the curvature. Although the fluctuation of Z 0 due to irregular vegetation patches has been validated in the field (King et al 2006, Kono andOkuro 2021), to our knowledge, this is the first time that the variability of Z 0 has been captured by an intuitive parameter derived from aerial images. This also illustrated that a constant assumption is reasonable for regular patterns where the variability of Z 0 is small, whereas the assumption is not tenable for irregular patterns.
The spatial structure of vegetation has long been linked to the stability or stable states of ecosystems. Based on a positive feedback mechanism, Okin et al (2015) pointed out that the connectivity is related to the shift of ecosystem states. That is, the larger the erodible bare soil patches are, the more easily resources are lost, vegetation degradation exacerbates as a result, and bare soil patches are further developed. Meanwhile, the patch size distribution deviates from a power-law has been hypothesized as an indicator of imminent degradation (Kéfi et al 2014). In our study area, aggregated soil particle size was closely correlated to the spatial distribution of vegetation. Under wind erosion conditions, vegetation patterns influence the deposition of nutrient-rich soil particles (Suter-Burri et al 2013) and wind desiccation associated with evapotranspiration (Trautz et al 2017). In turn, the redistribution and loss of resources affects plant growth (Rietkerk et al 2004). Therefore, for vegetation management in arid areas subject to wind erosion, the optimization of vegetation patterns is of great importance. Although we did not find an optimal combination of the spatial parameters in the context of complex interplay between vegetation and the environment, on the basis of the model presented here, the link between wind erosion and vegetation stability is strengthened. According to a field study, the placement of barriers to adjust the connectivity of desert vegetation played a significant role in the recovery of degradation (Peters et al 2020). Future studies should consider both the effect of the regularity of vegetation patches and the connectivity of bare soil as a reference for desert vegetation restoration.

Conclusions
The heterogeneity of the rough surface is better quantified by using connectivity and curvature to describe the spatial distribution of both bare soil and vegetation patches. Spatial vegetation structure affected the temporal and spatially averaged Z 0 significantly at the field scale. The accuracy of models with constant assumption of Z 0 in evaluating wind erosion depend on the regular of vegetation patches. For irregular vegetation patterns, including the connectivity and curvature can improve the estimation of Z 0 variability. Our findings strengthen the link between vegetation stability and wind erosion, which holds promise to support the wind erosion management practices in drylands.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.6084/m9.figshare.20438658.v1. Data will be available from 5 December 2022.

Authors' contributions
Conceptualization: X W and T C Data curation: T C and Y P Formal analysis: T C, Y P and Y L Methodology: X W and T C Writing original draft: X W and T C All authors contributed critically to the drafts and gave final approval for publication.