Comparisons suggest more efforts are required to parameterize wind flow around shrub vegetation elements for predicting aeolian flux

Upon interacting with the atmosphere, vegetation could alter the wind distribution and consequently the erodibility of nearby region. The parameterization of wind distribution around vegetation is crucial for the prediction of surface aeolian flux. This paper compared the performances of existing empirical distribution models in the estimation of aeolian flux for shrub vegetation, focusing on distribution pattern and vegetation porosity (main parameter of distribution function). Predicted dust fluxes directly entrained by air flow show weak sensitivity to both distribution pattern and porosity in the case of low vegetation density, which suggests some aspects in dust forecast models might be simplified. However, both distribution pattern and porosity show significant effect on sand saltation transport rate in the lee of vegetation element and, consequently, on the formation and evolution of surface aeolian landforms. The contribution of dust fluxes released in wind increase zone to the total emission by using current parameterizations increases with both the decrease of wind speed and the increase of vegetation density. Nevertheless, the parameterization of wind increase zone needs to be validated and improved by further experimental and numerical investigations.

Vegetation plays a very important role in the change of global climate and the sustainability of ecological system [1][2][3][4] . Firstly, as its basic physiological function, vegetation could fix the carbon dioxide in the atmosphere 5 and hold moisture in the soil 6 . Secondly, as its physical function -being regarded as ground obstacles, it could hinder or reduce the surface aeolian flux through direct or indirect manners [7][8][9][10][11] , which is particularly important for arid or semi-arid regions 1,4 . For example, previous studies suggest that aeolian dust aerosol has direct and semi-direct effect on the climate of Asian arid and semiarid regions [12][13][14] . Vegetation could directly prevent surface from erosion because of bestrow. It also could form a speed reduction region in the leeward side through interacting with air flow to indirectly decrease the erosion of surface soils 7,8,10,11 . Therefore, it is the key to find a reasonable way to combine the two manners together for quantitatively determining the effects of vegetation on wind erosion 15,16 . Shear stress partitioning model 1,15,16 is the most widely used model for qualifying the effects of vegetation on wind erosion. Its core concept is assuming that the dynamical effect of vegetation to an erodible surface is to increase threshold shear velocity (u *ft ) by absorbing momentum from air flow. Thus, the threshold shear velocity of soil surface covered by vegetation (u *ft1 ) is corporately determined by threshold on bare soil, vegetation basal cover (C), lateral cover (λ), the ratio of the drag coefficient of roughness element to that of unvegetated soil (β) and tuning parameter m 15 . The implement of shear stress partitioning model into wind erosion forecast model have been indeed improved the prediction of dust emission 1 . However, there are still uncertainties existing in estimations of the magnitude of dust events 17 or of the exact location of dust sources 18 . These uncertainties aren't induced by model precision, but are more likely to originate from the shear stress partitioning model 16,18 . This is because that shear stress partitioning model assumed the same threshold everywhere in vegetated surface (uniform shear stress distribution). In fact, recent experiment 10 and simulations 11 all reveal that the shear stress in vegetated surface distributes heterogeneously (non-uniform shear stress distribution). The presence of vegetation reduces aeolian flux by altering (mainly decreasing in the leeward side) the shear stress distribution of local region but not raising the threshold velocity of the whole surface 8 . And a wind erosion model including the distribution of shear stress performs better than that employing stress partitioning model 8,19 . Hence, the parameterization of shear stress distribution around vegetation becomes the key of all issues.
The parameterization of shear stress distribution includes two aspects -distribution pattern and distribution function. Typically, the shear velocity is roughly proportional to wind speed at a reference height above surface 10 , hereafter, wind speed distribution is used instead of shear stress distribution. For convenience, a vegetation element is usually simplified as a cylinder 8 . There are three main distribution patterns. The first pattern is the triangle shape proposed by Raupach 7 . He supposed that the wind speed doesn't vary continuously but is all zero within the zone, while equals the coming wind speed outside the zone. Nevertheless, this assumption seems to be more applicable for solid cylinder cases but less effective for a vegetation element with porous structures 8 . The second one is the rectangular shape proposed by Okin 8 (Fig. 1a,c). He supposed that the wind speed recovers gradually from the lowest value back to the coming wind speed within the zone. The last one is proposed by Leenders et al. 9 (Fig. 1b,d). Different from the two formers, they proposed that two zones, wind speed reduction zone and increase zone, exist around vegetation. The wind speed reduction zone is assumed to be half-ellipse, while the speed increase zone to be full-ellipse. There are many studies on the wind speed distribution in the leeward side of roughness elements [8][9][10][11][20][21][22][23][24][25][26] . But systematic and comprehensive description of wind speed around porous vegetation elements is still scare. Here, three typical distribution functions of wind speed reduction for vegetation element are focused. The first function is proposed by Okin 8 through fitting the experimental data behind a porous windbreak from Bradley and Mulhearn 20 . The second one is proposed by Leenders et al. 9 through fitting their measuring data around shrub vegetation elements (Hyphaene thebaica and Commiphora africana). But, they only modified the parameters of Hagen's distribution function 8 , remaining the form unchanged. The last one is recently proposed by Mayaud et al. 10 through fitting their field observed data around three shrub types (Stipagrostis amabilis, Rhigozum trichotomm, and Zygophyllum stapfii). However, they chose the distinguishing form of distribution function in comparison to the first two types. Furthermore, to be more reasonable, Leenders et al. 9 and Mayaud et al. 10 selected both plant height and porosity as the main factors controlling the variation of wind speed in their distribution functions, while Okin 8 only employed plant height. In contrast, although wind increase zone has also been observed in other studies 10,11 , the description of this zone for porous vegetation was only conducted by Leenders et al. 9 .
After a review of literature, it could be found that there is still a lack of studies on comparing the performances of these proposed parameterizations in predicting aeolian flux. Firstly, although the difference in wind speed reduction among different distribution functions was studied 10 , nevertheless, the sensitivity of dust flux release in a vegetated surface to the parameterizations of wind speed distribution (including both pattern and function) isn't clear yet. Secondly, the importance of the wind increase zone to the estimation of total dust flux in study area isn't www.nature.com/scientificreports www.nature.com/scientificreports/ fully understood, and it is also required to evaluate how the sensitivity of increased dust fluxes to parameters that are used to describe the wind distribution in the zone. Thirdly, the sensitivity of saltation flux in the leeward side, which is crucial for the development and evolution of coppice dunes 27,28 , to the parameterizations of wind speed distribution, should be tested. Therefore, this paper aims to investigate these issues.

Methods
In nature, vegetation elements (or patches) distribute randomly 1,8 . Here, for simplification, the cylinder vegetation elements are assumed to be uniformly distributed in a square area (defined as l × l) 25,26 , as shown in Fig. 1. The vegetation height is denoted as H, and the diameter is denoted as D. The basal cover and lateral cover are defined as C = πD 2 /(4l 2 ) and λ = DH/l 2 , respectively. The coming wind is supposed to be unidirectional and statistically stable in time and space, which means that no turbulence is considered.
Patterns and functions of wind speed distribution around vegetation elements. As stated above, two patterns (proposed by Okin 8 and Leenders et al. 9 ) are selected here for porous vegetation elements. For rectangular shape pattern (Fig. 1a), the width of the zone equals the diameter of vegetation elements, and the length is the distance, L x , from the central location of vegetation to the location at which the wind speed recovers to the coming wind speed. For the half-ellipse reduction zone (Fig. 1b), the semi-minor axis and semi-major axis are 0.5D and L x , respectively. For the full-ellipse increase zone, the semi-minor axis and semi-major axis are 0.25D and 0.5D, respectively.
Here, the central location of a vegetation element is set as the origin of coordinate, and the wind speed around the element is denoted by * u x y ( , ). The wind distribution function proposed by Okin 8 for wind reduction zone is shown in Eq. (1), where C O1 = 0.32, C O2 = 4.8, u *ref is the referring shear speed of incoming wind, and x is the horizontal coordinate. The distribution function proposed by Leenders et al. 9 for wind reduction zone is expressed in Eq.
. The distribution function for wind increase zone is described in Eq. (4). Ø and C p are the wind-increase factor and area-increase factor, respectively. x 0 and y 0 are the coordinates of the central location of wind increase zone.
Erosion-related quantities in this paper. Some abbreviations are given at first to help compare performances of published parameterizations of wind speed distributions. In order to understand the differences better, four groups of combinations (pattern + function) are employed. Group 1st is abbreviated as "O2008" which means that both pattern and function are from Okin 8 . Group 2nd is abbreviated as "OM2017" which means that pattern is from Okin 8 and function is from Mayaud et al. 10 . Group 3rd is abbreviated as "L2011" which means that both pattern and function are from Leenders et al. 9 . Group 4th is abbreviated as "LM2017" which means that pattern is from Leenders et al. 9 and function is from Mayaud et al. 10 . Five erosion-related quantities are focused in this work. The first quantity is the protecting efficiency (P r ) that describes the proportion of the non-erodible area to the total study area. This non-erodible area, within which the wind speed is always lower than threshold value, includes the basal area and the area caused by the sheltering role of vegetation. Mathematically, are the total study area (l × l), the erodible area, and the non-erodible area, respectively. The second one is the averaged release rate (F vt ) of fluid-entrained PM10 dust in study area. The averaged release rate contains three parts: the release in wind increase zone (F in ), the release in wind reduction zone (F re ), and the release in normal wind zone (F nor ), as shown in Eq. (5). The fluid-entrained PM10 flux could be estimated by the Eq. (6) that is proposed by Zhang et al. 29 on the basis of wind tunnel experimental data. The third one is the surplus PM10 flux (TF) defined as the differential of fluxes between considering and without considering the increase of wind speed in wind increase zone, as shown in Eq. (7). The fourth one is the proportion P in of the TF versus the total flux that doesn't consider the effect of wind increase, as shown in Eq. (8). The last one is the reduction of transport rate (Q r ) defined as the differential of rates between considering and without considering the decrease of wind speed within wind decrease zone, as shown in Eq. (9). It's employed to statistically evaluate the integrated reduction level of transport rate caused by www.nature.com/scientificreports www.nature.com/scientificreports/ the sheltering role of a single vegetation element. The saltation transport rate q x is estimated by Eq. (10) 30 , where u *it is the impact entrainment threshold, g the gravitational acceleration and ρ a the air density. q x,0 tells the estimated rate by using the referring shear speed, while q x,re by using the reduced shear speed.
Other settings. Shrub

Results
At the beginning, the proposed parameterizations of leeward wind distribution are compared (Fig. 2). For convenience, OE2008, LE2011 and ME2017 are employed to represent the wind distribution functions proposed by Okin 8 (Eq. 1), Leenders et al. 9 (Eq. 2) and Mayaud et al. 10  A vegetation element could extend its effective protecting area, where no wind erosion occurs, by the sheltering effect. But, the degree of extension depends on wind speed, as shown in Fig. 3a. With the increase of incoming wind speed, the protecting efficiency P r decreases towards basal cover over all groups. Nevertheless, there are still some differences. At the lowest wind speed, the P r of rectangular distribution pattern (O2008 and OM2017) reaches around 1.2 times higher than that of half-ellipse pattern (L2011 and LM2017). Different from LM2017 that P r decreases gradually, a sharp decrease of P r occurs between wind speed 8~9 m/s in L2011. This difference should be caused by the wind speed distribution function. In contrast, the averaged release rate of PM10 (F vt ) increases with the incoming wind speed, but shows very weak sensitivity to both distribution pattern and function (solid scatters in Fig. 3a). It also reveals weak sensitivity to vegetation porosity (scatters in Fig. 3b), no matter www.nature.com/scientificreports www.nature.com/scientificreports/ what the wind speed is. Besides, Fig. 3b shows two remarkable differences in prediction of dust release between by considering the distribution of wind speed (scatters) and by without considering the distribution of wind speed (i.e., using shear stress partitioning model) (solid line). Firstly, the estimated PM10 fluxes by considering wind speed distribution are generally higher than those by without considering wind speed distribution. Secondly, it can be seen that, there is no emission of PM10 when the wind speed is below 8 m/s for without considering wind speed distribution, while the emission of PM10 still happens by considering wind speed distribution.  www.nature.com/scientificreports www.nature.com/scientificreports/ It reveals that, at a fixed lateral cover, the relative proportion (P in ) of TF reduces with the increase of wind speed for both distribution functions (Fig. 4a). And, the difference in the relative proportion between two functions is negligible. Also, as shown in the figure, the lateral cover affects the relative proportion significantly. At a fixed wind speed, the relative proportion increases with lateral cover. For example, the value of P r in the case of λ = 0.125 is 6~7 times larger than that in the case of λ = 0.03125, which seems to be invariant with wind speed. The tests of the sensitivity of TF to main parameters (φ and C p ) are shown in Fig. 4b. It reveals that, the TF increases with both φ and C p . However, in comparison to area-increase factor C p , TF shows higher sensitivity to wind-increase factor φ. Exactly, curve fit of data (R 2 > 0.98) suggests TF increases with C p linearly, while increases with φ exponentially. Thus, the ratio of the largest value to the lowest value for TF is about 1.5 within the given range of C p , but that for TF is higher than 5 within the given range of φ.
The sand transport rate focused here is calculated on the basis of the hypothesis that only one vegetation element exists in the whole flat land surface. The wind distribution function proposed by Okin 8 shows no strict constraint on the end point of sheltering zone. The calculating results in Fig. S1 in supporting information suggest that the reduction of transport rate (Qr) reaches saturation when the leeward distance is about 20 H. So, the sheltering length is approximately set to be 20 H when the transport rates are calculated by using Okin's distribution function. As shown in Fig. 5a, Qr increases with wind speed over all groups. The transport rates calculated by rectangular pattern (O2008 and OM2017) are far higher than those by half-ellipse pattern (L2008 and LM2017). Under the same pattern, the difference in transport rates owing to distribution functions is small. Figure 5b shows the change of Qr with vegetation porosity by using two different wind speed distribution functions. Although Qr decreases with vegetation porosity (because the leeward wind speed would increase with vegetation porosity according to proposed parameterizations), the response of transport rate to porosity is both wind-dependent and function-dependent. It could be found that, the increase of incoming wind speed could enhance the variation of Qr caused by vegetation porosity (by comparing the square scatters with the triangle scatters). Generally, Qr calculated by the distribution function of Leenders et al. 9 is less sensitive than that by the distribution function of Mayaud et al. 10 . For example, in the case of U z = 15 m/s, Qr decreases about 15% (referring to the maximum value) for L2011 within the given range of porosity, but decreases about 80% for LM2017 within the same range.  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion and Conclusion
Prediction of dust events is not only important for human being's health and safety but also for our understanding of the climate change of our planet 3,12-14 . It is thus crucial for us to accurately estimate the dust load in vegetated land surface, particularly in arid and semi-arid region 3,4 . As the most important driving factor for aeolian dust emission, three selected parameterizations (OE2008, LE2011 and ME2017) of leeward wind speeds were compared firstly. Significant differences in the lowest wind speed as well as the corresponding location exist among the three parameterizations, which should originate from the difference in vegetation canopy shape [8][9][10] . However, the total release rate of dust flux, which is directly entrained by air flow, shows weak sensitivity to both distribution pattern and porosity in the case of low vegetation density (low coverage). Although the contribution of dust release in wind increase zone is considered, the proportion of dust release in wind increase zone accounts for about 5% of total dust emission in the case of low vegetation density or high wind speed, which wouldn't alter the weak sensitivity. This thus suggests that some aspects in dust forecast models might be simplified. For instance, it seems that the difference in canopy shape owing to vegetation types as well as their growth stages couldn't need to be considered in low vegetation density locations. Nevertheless, it doesn't mean that the dust release predicted by shear stress partitioning model (uniform distribution of shear stress 7,15 ) would agree with that by the method here (non-uniform distribution of shear stress 8,19 ), in the case of low density. Calculated results reveal that dust release could occur in much lower wind speed by considering non-uniform distribution of shear stress than that by uniform distribution of shear stress. Therefore, employment of uniform distribution of shear stress is one of possible reasons that why current dust forecast models yield uncertainties in the location and magnitude of dust events 16 .
The importance of wind increase zone in the lateral side of vegetation elements is tested in this work. Under the settings proposed by Leenders et al. 9 , the contribution of wind increase zone in total dust release gradually increases with both the increase of vegetation density and the decrease of incoming wind speed. When the average wind speed over the whole study area is below fluid threshold (no matter this is caused by either the low incoming wind speed or the high vegetation cover), the increase of wind speed around vegetation could make the local speed raise up to a value above fluid threshold and consequently lead to dust emission. So, the increase of wind speed around vegetation element may be another reason that causes the uncertainties by current dust forecast models. Calculating results also reveal the sensitivity of increased dust flux to the two main parameters (wind-increase factor and area-increase factor). It is thus important and valuable to make sure whether the proposal on the parameterization of wind increase zone is reasonable or not. Mayaud et al. 10 also detected the wind increase in the lateral side, but they didn't quantitatively describe the extent and magnitude of wind increase zone. Their study further suggested that the ratio of height to diameter (or width) and the porosity could affect the acceleration of wind speed. Hence, it is required more data (including experimental, numerical and theoretical data) to confirm the parameterization of wind increase zone by including more related physical factors.
The change of sand transport rate is the key to the formation and development of coppice (or nebkha) dunes in the lee of vegetation element 27,28 . Although total dust release rate is weakly sensitive to both distribution pattern and function, the reduction of transport rate is sensitive to them. For a single element, the half-ellipse pattern seems to be more reasonable than rectangular pattern, according to recent measurement 10 and numerical simulation 26 . The difference in sensitivity of reduction of transport rate to vegetation porosity, which may be caused by the different canopy shape between LE2011 9 and ME2011 10 , suggests that more efforts are required to carefully determine the effect of porosity on the recovery of leeward wind speed.
Furthermore, two additional points needs to be noticed. Firstly, the possible change of distribution pattern or function due to the increase of vegetation coverage should be investigated. It should be reminded that these results obtained in this work may be only applicable in low density cases. The interaction among vegetation elements could greatly change the flow pattern (e.g., skimming flow) in moderate or high density cases 34 , consequently, the distribution pattern or function couldn't be suitable any more. Besides, experiments 35,36 suggest the canopy shape has great effect on vegetation drag coefficient; and numerical simulation 11 shows the effect of canopy shape on leeward wind speed distribution. The wind speed at the lowest grid in dust forecast models might thus be estimated inaccurately with the increase of vegetation density because of the difference in canopy shape, which could also be able to rise up the uncertainties in prediction of dust events. Also, the sensitivity of total fluid-entrained dust flux to both wind distribution pattern and porosity in moderate or high vegetation density may be different from that in low vegetation density as shown above. Secondly, the effect of height on the sheltering role of vegetation should be constrained. All distribution functions introduced here didn't show any constraint on the effect of height, indicating that the sheltering role of vegetation is proportional to vegetation height. In another word, if the height is infinite, the sheltering role would be infinite. This is not the truth. For a fixed diameter (or width), the effect of height on the sheltering role is limited. Raupach 7 proposed an empirical expression to conceptually describe the variation of sheltering with the ratio of height to width. However, the parameters in the expression are difficult to obtain. Sadique et al. 26 recently proposed a simple way (dividing the effect at the location of roughness height) to limit the effect of element height, based on the numerical data of rectangular solid roughness elements. Thus, their try affords us lesson and could be altered to apply in the study of porous vegetation element.
In conclusion, this paper compared the performances of existing empirical distribution models in predicting aeolian flux for shrub vegetation, focusing on distribution pattern and vegetation porosity (main parameter of distribution function). Canopy shape shows significant influence on wind distribution around vegetation element. Results reveal that the predicted protecting areas would vary, to some extent, among different models; however, the predicted dust fluxes directly entrained by air flow are weakly sensitive to both distribution pattern and porosity in the case of low vegetation density, regardless of canopy shape. The predicted dust fluxes in wind increase zone mightn't be overlooked, particularly in the cases of low incoming wind and high density conditions. The dust fluxes in wind increase zone show the sensitivity to both area-increase factor and wind-increase factor, but are likely to be more sensitive to the latter. The predicted sand saltation transport rate in the lee of a vegetation element is sensitive to both distribution pattern and porosity. The difference in the reduction of transport rate www.nature.com/scientificreports www.nature.com/scientificreports/ owing to distribution pattern of wind reduction zone increases with wind speed, while the variation of transport rate with porosity is function-dependent. Rectangular pattern produces more reduction of transport rate than half-ellipse pattern. The Qr predicted by the distribution function of Mayaud et al. 10 is more sensitive to vegetation porosity than that by the distribution function of Leenders et al. 9 . This difference in sensitivity comes directly from the different leeward wind speed variations. These variations are controlled by the role of porosity in distribution functions that is affected by canopy shape (or shrub type). It thus implies the necessity of considering the canopy shape (or shrub type) in further understanding of the form and evolution of leeward surface landforms. These results shown above suggest that, we need more experimental measurements or numerical simulations to confirm or improve the wind distribution model around shrub vegetation, in order to exactly determine aeolian flux in the presence of vegetation.