The Impacts of Vegetation and Meteorological Factors on Aerodynamic Roughness Length at Different Time Scales

: The aerodynamic roughness length (z 0m ) is a crucial parameter for reliably simulating turbulent exchanges between the land surface and the atmosphere. Due to the large number of input variables related to vegetation growth and aerodynamic conditions near the surface, estimating z 0m precisely is difﬁcult and, to date, no universal model has been established. Understanding the z 0m changes in time series data and the relative contributions of vegetation indices and meteorological factors is important to providing a basis for modelling z 0m . In this paper, the main meteorological factors that inﬂuence z 0m in different seasons are presented based on data from three automatic weather stations (AWSs) that represent various land surface patterns in the Heihe river basin. A correlation analysis identiﬁed the dominant factors that inﬂuence z 0m changes at half-hour and daily scales; then, a factor analysis was performed to identify the different contributions of vegetation indices and meteorological factors to z 0m at different time scales. The results show that meteorological factors (wind speed, wind direction and atmospheric stability) are the main driving factors for z 0m at the Arou and Guantan sites, which are situated in grassland and forest mountain areas, respectively, and that the vegetation indices have no impact on the z 0m variations in these areas. In contrast, for the Daman site, situated in ﬂat farmland, the vegetation indices are the primary driving factors, while meteorological factors such as wind speed and atmospheric stability are secondary factors, and wind direction has no signiﬁcant inﬂuence. Finally, a detailed analysis was conducted to detect the relationships between half-hourly z 0m measurements and three dominant meteorological factors.


Introduction
The aerodynamic roughness length (z 0m ) is the height above the ground where the wind speed is zero under neutral conditions [1]. Aerodynamic roughness length is a key variable in land surface flux simulations and atmospheric boundary studies [2,3]. The z 0m metric is a virtual height under ideal atmospheric conditions that cannot be measured and quantified directly; therefore, it is important to develop a reliable parameterized method for estimating z 0m at different spatial and temporal scales.
Micrometeorological wind profiles measurements have generally been used to estimate site-specific z 0m based on the Monin-Obukhov similarity theory [4,5]. However, because modern studies involve land surface processes at a greater spatial scale, in recent years, an increasing number of models have been established to estimate regional z 0m . Raupach developed an approach for estimating z 0m in terms of height, width, spacing and wind-drag properties of surface roughness elements [6,7], Jasinsk simulated the parameters in the Raupach model by classifying four vegetation types [8]. With the development of remote sensing techniques, links have been built between regional z 0m and vegetation indices such as the normalized different vegetation index (NDVI) and leaf area index (LAI) [9,10]. In addition, the signatures of bidirectional reflectance distribution function (BRDF) derived from multi-angle remote sensing data have been used to describe the canopy structural components of roughness length [11]. However, previous studies of modelling z 0m mainly considered topographic and vegetation features on the surface and neglected the effects of aerodynamic factors, which results in significant errors [12]. Schmid found that z 0m is not only related to the vegetation geometry but is also affected by the near-surface wind speed [13]. Zhou found that z 0m changes with wind direction and atmospheric stability at a given location [14,15]. However, few quantitative studies have been performed to evaluate the proportions of all the factors that may influence the instantaneous or daily average z 0m estimates over different underlying surfaces.
Therefore, to fill this gap in the field and provide a reference for modelling regional z 0m , this paper investigates the driving factors of z 0m using onsite time-series data over three typical vegetation types. Depending on research and application needs, z 0m can be summarized to hourly or daily data or even fixed within a specified season; we select half-hourly z 0m and daily z 0m as the main objectives in this study. Several statistical methods are available to quantify the effects of various factors on z 0m , including regression analysis, principal component analysis (PCA), factor analysis and artificial neural network analysis. Considering that all observed variables mainly reflect changes of z 0m in a certain place in two ways-growth of vegetation and variations of atmospheric conditions-we adopt factor analysis in this study to search for joint factors and evaluate the contributions of factors related to z 0m .
In recent decades, intensive ground observation facilities supported by continuous research projects have been established in the Heihe River basin, and these facilities have collected multiple meteorological data time-series [16][17][18] which are convenient for analysing the differences between the factors that influence z 0m over different landform and vegetation types in a single climate zone. For this study, we selected the Guantan, Arou and Daman stations, which respectively represent underlying surfaces of forest, grassland and cropland. Vegetation indices calculated from remotely sensed images are used to describe the vegetation growth. Finally, a detailed analysis of the variation characteristics of half-hourly z 0m is presented, along with each main driving meteorological factor.

Site Description
In this study, data were collected from field automatic weather stations (AWSs) at three sites in the Heihe River basin distributed in the upper and middle reaches ( Figure 1). Arou station is located in an east-west oriented valley with a maximum width of 3 km from north to south. The site is surrounded by natural alpine meadow with grass reaching a height of approximately 20-30 cm during the growing season. The Guantan site is surrounded by continuous mountains with elevations ranging from 2700-3400 m, and the weather station is located on the shady slope of the mountain in the Dayekou sub-catchment. The overstory tree species at the study site is Qinghai spruce, and the forest floor is nearly covered by a layer of moss with a depth of 10 cm [19]. Daman station is located in a typical oasis on very flat terrain near Zhangye City. The underlying surface is mainly irrigated maize that reaches a maximum height of approximately 180 cm [20]. The elevation map was produced using the ASTER GDEM product of METI and NASA (available online: http://www.gdem.aster.ersdac.or.jp).

Onsite Data
The AWSs at these three sites collected 144 records per day including wind speed, wind direction, air temperature and air humidity at 10-min intervals at different heights. A description of the AWSs and the time periods covered by the data used in this study are listed in Table 1. Remotely sensed vegetation indices including the normalized difference vegetation index (NDVI) and leaf area index (LAI) were selected as the key indicators to measure the characteristic of vegetation growth. The time-series cloud-masked NDVI values of the sites were calculated using the MODIS red and infrared band reflectance after atmospheric correction: 2010 for the Guantan station and 2012 for the Arou station with 1 km resolution. The Savitzky-Golay (SG) filter was adopted to fill gaps to acquire the daily NDVI [21]. The 1-km/5-day composite LAI product in the Heihe River basin was acquired from the China science data centre (available online: http://card.westgis.ac.cn/), which consists of a normalization algorithm to integrate Terra-MODIS, Aqua-MODIS, FY3A-MERSI and FY3B-MERSI data. The time-continuous LAI were retrieved using neural networks and a maximum-value compositing algorithm [22][23][24].

Calculation of the Obukhov Length L
The Obukhov length L near the surface can be used to characterize the effect of atmospheric stability on the turbulence [25]. Therefore, this study investigates the implications of vertical measurement profiles in relation to the Monin-Obukhov similarity, using the gradient Richardson number (Ri) to infer L with the following equations : where z = (z 1 − z 2 )/ln(z 1 /z 2 ) is the height where the calculated Richardson number is valid, z 1 is fixed at 10 m height, while z 2 is 2 m for Arou and Guantan site and 3 m for Daman site, respectively. Compared to L, Ri has the advantage that it contains only meteorological gradients that can be determined experimentally. Ri is expressed as the non-dimensional ratio between the buoyancy forces and the inertial ones due to the wind [26]. The sign of this number is determined by the temperature gradient, the negative values are an index of instability and the positive ones are an index of stability.
where T is a reference atmospheric temperature, ∂θ v /∂z is the difference in virtual potential temperature (VPT) between two measurements at heights z 1 and z 2 , g = 9.8 m/s 2 is the gravitational acceleration, and ∂u/∂z is the wind gradient between z 1 and z 2 .

Calculating the Aerodynamic Roughness Length from AWS Data
To guarantee the reliability of the wind profile results, the raw datasets were pre-processed and quality control was performed to obtain aerodynamic roughness values based on the following criteria: (1) the minimum wind velocity of all layers exceeded 2 m/s; (2) data obtained at rainy times were discarded; (3) missing or invalid data were interpolated based on neighbouring data [18,27].
The aerodynamic surface roughness length can be determined iteratively based on wind profile data. Using the Monin-Obukhov similarity, the z 0m values can be calculated via the logarithmic wind profile equation [6,7]: where u and θ denote the wind speed and potential air temperature, respectively, at the height z above ground level, u * is the friction velocity, θ * is the friction temperature; k is von Karman's constant (k = 0.4); d is the zero-plane displacement; z 0m is the aerodynamic roughness length; z 0h is the thermal roughness length and θ 0 is the potential temperature near the surface. The expressions of the stability functions Ψ m and Ψ h depend on the stability conditions in the surface layer, which are described by the stability parameter Z/L [28]. When Z/L < 0 (unstable conditions) [29], When Z/L > 0 (stable conditions) [30], In this study, we use the least-squares method to determine z 0m [31,32]. Equation (3) can be written as the following simplified form: Commonly, d is related to the vegetation height (h) such that d equals 0.67 h over surfaces covered by dense and homogeneous vegetation [33]. Therefore, d is determined here by iterating over 0.1 m intervals of from 0.1 m to 3 m. For each given d, a correlation coefficient exists for the fitting Equation (9), and the optimal d value is chosen when the correlation coefficient reaches its maximum, the invalid value will be returned in loop iterations if the wind speed gradient is abnormal. Finally, the 10-min z 0m values can be calculated according to the optimal d value. To maximally avoid effects from missing data, the z 0m data are averaged over 30-min intervals and daily average z 0m values are also calculated.

Correlation Analysis
Correlation coefficients measure the strength and direction of a relationship between two variables. Here, Pearson product-moment correlation coefficients were calculated to investigate the relationship between the latent factors and z 0m . We selected the NDVI and the LAI as the most likely factors that represent the effects of vegetation growth on z 0m and adopted wind speed (WS), wind direction (θ), relative humidity (RH), air temperature (T), Obukhov length (L) and air pressure (P) as the meteorological factors. Because Obukhov length is an instantaneous parameter that describes atmospheric conditions, calculating an average daily level of Obukhov length is meaningless and excluded from the correlation analysis of daily z 0m . To match the half-hourly z 0m data, we used daily NDVI and LAI values under the assumption that changes in vegetation indices can be ignored over shorter periods. Because wind direction is a categorical variable, we proposed a definition of wind deflection angle (∆θ) between the prevailing wind direction (θ p ) and the instantaneous wind direction, described as follows: The prevailing wind direction usually depends on the local topography, latitude and the locations of land and sea; therefore, ∆θ is categorized as a topographic factor in the subsequent factor analysis. The prevailing wind directions at the Arou, Guantan and Daman sites are 162 • , 121 • and 318 • relative to due north, respectively, according to all observation data.

Factor Analysis
Factor analysis is a statistical method for decomposing input data as a linear combination of common factors and specific factors in a lower dimensional space. As a feature extraction approach, it seeks to acquire factor scores for input data in a lower dimensional space via a linear transformation. Those factor scores can be viewed as extracted features of the original data in another, lower-dimensional space. In this study, we first used correlation analysis to determine the main driving factors that influence the changes in z 0m at half-hourly and daily intervals. Then, the main driving factors were normalized to eliminate the differences in dimensions between the variables: where X j is the indicator value of the j factor, and X j and σ j are the average and standard deviation of all the sample data of the j factor, respectively. Then, factor analysis was performed to determine the contributions of the vegetation indices and meteorological factors that affect z 0m . The mathematical model of factor analysis is as follows [34]: where X = (X 1 , X 2 , . . . , X j , . . . , X k ) is a k-dimensional stochastic vector composed of k samples, F = (F 1 , F 2 , . . . , F n ) is the common factor of X (n < k), ε k×1 is an unobserved stochastic error factor with zero mean and finite variance and A k×n is the correlation coefficient matrix of the top n factors: where λ 1 , λ 2 , . . . λ n are the first n characteristic values of the covariance matrix X k×1 , and u 1 , u 2 , . . . u n are the corresponding standard orthogonal eigenvectors. Based on the number of eigenvalues that exceed 1, the amounts of the common factors, n, can be determined. Then, the variance contribution of the common factor F n×1 is obtained by calculating the sum of the squares of the elements in the j-th column of A k×n : which represents the sum of the variance contributions of the same common factor Fj to the variable. The cumulative contributions of the common factors can then be calculated accordingly. The factor analysis was performed using SPSS version 25 software (IBM, New York, NY, USA), using PCA as the extraction method and varimax as the rotation method.

Driving Factors for Changes in Half-Hourly and Daily z 0m
The results of bivariate correlation analyses between half-hourly z 0m , daily z 0m and each possible driving factor are shown in Table 2, which indicates that the main driving factors for z 0m differ among the three stations. For the Arou and Guantan sites, the half-hourly and daily z 0m values are positively correlated with the wind deflection angle and negatively correlated with the wind speed, while no significant correlation exists with the NDVI and the LAI, which represent the vegetative growth conditions. In contrast, at the Daman site, the half-hourly and daily z 0m values are positively correlated with the NDVI and the LAI and negatively correlated with wind speed, while the wind deflection angle has no correlation to z 0m . Similarly, half-hourly z 0m is significantly correlated with the Obukhov length; however, neither half-hourly z 0m nor daily z 0m have any significant correlation with relative humidity, air temperature or air pressure at any of the sites. The factors that are significantly correlated with half-hourly z 0m and those that are significantly correlated with daily z 0m were selected to conduct a subsequent factor analysis. As Table 2 shows, only two factors have significant correlations with daily z 0m at the Arou and Guantan sites; therefore, it is unnecessary to conduct a factor analysis. Consequently, we conducted the factor analysis for daily z 0m only at the Daman site and for half-hourly z 0m at all three sites.

Contributions of Driving Factors to Half-Hourly z 0m
Before conducting the factor analysis, Kaiser-Meyer-Olkin (KMO) [35] and Bartlett's sphericity tests [36] were applied to examine the appropriateness of the data for factor analysis. The KMO test determines the degree of correlation between variables with partial correlation coefficients: the KMO values for the half-hourly z 0m values at the Arou, Guantan and Daman sites were 0.625, 0.591 and 0.703, respectively, all of which satisfy the minimal requirement; that is, they are above 0.5. Bartlett's sphericity test determines whether variables are independent based on the variables' correlation coefficient matrixes. The Bartlett's sphericity test results were 67.222, 58.761 and 121.928 at the Arou, Guantan and Daman sites, respectively, and the p-value for each site was 0: all of these satisfy the requirement that p should be below 0.001. Based on the results of these two tests, the data are suitable for factor analysis. Tables 3 and 4 show the factor analysis results. For half-hourly z 0m at the Arou and Guantan sites, the first common factors are wind speed and Obukhov length, which are aerodynamic factors (A_factor): their contributions for Arou and Guantan are 44.03% and 42.87%, respectively. The second common factor (the topographic factor (T_factor)) is mainly expressed by wind deflection angle, and the contributions at Arou and Guantan are 38.33% and 38.66%, respectively. For the Daman site, the first common factor is the vegetation factor (V_factor) consisting of NDVI and LAI: the factor loads for NDVI and LAI are 0.959 and 0.941, respectively. The contribution of the vegetation factor is 53.14%. The aerodynamic factor, composed of wind speed and Obukhov length, is a secondary factor that contributes of 30.29%.

Contributions of Driving Factors to Daily z 0m
The KMO value for daily z 0m at the Daman site is 0.731, which satisfies the minimal requirement (>0.5). The Bartlett's sphericity test value is 333.777 and the p value is 0. Both the KMO test and Bartlett's sphericity test demonstrate the feasibility of the data for factor analysis. Similar to the results of half-hourly z 0m , the first common factor for daily z 0m at the Daman site is the vegetation factor, which is composed of NDVI and LAI with factor loads of 0.969 and 0.960, respectively. The contribution of the vegetation factor is 54.46%. The contribution of the aerodynamic factor is 30.67%.

The Applicability of Monin-Obukhov Similarity
In this research, we adopted a traditional method of calculating observed z 0m by fitting the profile observations using Monin-Obukhov similarity (MOS). As MOS was developed under homogeneous and flat underlying conditions, the applicability of MOS at the Arou and Guantan sites, which are located in mountainous areas, needs to be examined further. We filter the profiles with minimum wind velocity over 1 m/s and separate the wind speed data of one certain layer (20 m at the Arou site and 2 m at the Guantan site) as the target for validation. The remainder of the profiles are used to calculate the friction velocity, the displacement height and the aerodynamic roughness length in the MOS framework. Then, the wind speed of the target layer can be calculated accordingly. Figure 2 shows the linear fitting relationships between the measured wind speed and the calculated wind speed of the target layer: the coefficients of determination (R 2 ) values of the fitted equations are 0.95 and 0.68 for the Arou site and the Guantan site, respectively, indicating significant relationships between the measured wind speed and the calculated wind speed at both sites, confirming the applicability of MOS.

The Wind Speed Factor
To formulate the detailed change rule between z 0m and wind speed, we select the measured AWS data in August, when the vegetation reaches its peak-growth stage. The averaged z 0m values and standard deviations are calculated at 0.1 m/s wind speed intervals. Figure 3 reveals that the z 0m values decrease with increasing wind speed. The reduction in z 0m as the wind speed increases is most significant at the Arou site, and to a certain extent, the z 0m values remain relatively low and stable when the wind speeds become high (Figure 3a). Because of the barrier effects of the densely distributed stands of spruce around the Guantan site, the maximum wind speed in August is just over 4 m/s; thus, the descending z 0m trend is not pronounced (Figure 3b). At the Daman site, z 0m has a slowly descending tendency that is significantly correlated with increased wind speeds (Figure 3c). At slow wind speeds, the standard deviations of z 0m at all the stations are generally larger, possibly because the weak wind creates irregular turbulence inside the vegetation, resulting in large variations of z 0m under similar wind speed conditions. The most likely explanations for the decrease in aerodynamic roughness with increasing wind speed are as follows: as the wind speed increases, the wind changes the geometry of nonrigid bodies (such as vegetation) above the ground surface to some extent. In other words, as the stalks or blades of a plant lean in the same direction under wind action, the airflow drag due to the rough elements decreases; as the wind speed continues to increase, there is a critical point at which the plants reach equilibrium and their shapes no longer change with the wind speed. When the observed wind speed is small, the z 0m measured at the Guantan and Daman sites never reach this stable stage, whereas, at Arou, the z 0m measurements reflect this equilibrium.

The Wind Direction Factor
The influence of wind direction under unstable conditions was also explored to assess the homogeneity of these three sites. We plotted the distributions of half-hourly wind speed under different wind directions on polar coordinates; the frequency percentages of z 0m were classified at eight common wind directions (N-NE-E-SE-S-SW-W-NW). As Figure 4a shows, the dominant prevailing winds are west and southeast at the Arou site. Although the wind speed values are similar, significant differences occur in the z 0m between these two wind directions (Figure 4b): under a west wind, 98% of the z 0m values are less than 0.04 m, while under a southeast wind, 58% of the z 0m values under exceed 0.04 m. Therefore, easterly winds tend to result in larger z 0m values at the Arou site. At the Guantan site, wind is concentrated in two opposite directions, which can be partially explained by the mountains around the site. The z 0m values under the different wind directions also vary significantly (Figure 5b). Under south and southeast winds, approximately 80% of the z 0m values are less than 0.4 m. In contrast, under a northwest wind, 90% of the z 0m values are greater than 0.4 m. These results occur because the Guantan site is located on the shady slope of the mountain; thus, wind blowing directly against the slope results in larger z 0m values. The Daman site is located in an irrigated flat farm area where the topography has little influence on the wind. The prevailing wind directions in August are southeast and northwest, and there is no significant difference in the frequency distribution of z 0m between the different wind directions, as shown in Figure 6b.    Table 5 shows the mean values of z 0m in August calculated at the two prevailing wind directions at each site and shows that the z 0m values at Arou and Guantan site differ significantly with wind direction, while those at the Daman site are similar even under two opposite wind directions. Considering that the topographic relief of the entire source area is more prominent at Arou and Guantan sites compared with that at the Daman site, the z 0m influences from changes in wind direction may depend on the heterogeneity of the topography. For areas with complicated terrain, the density and distribution of rough elements presented to each wind direction differ; therefore, the drag forces near the earth's surface also differ, resulting in z 0m variations under different wind directions. However, in the area surrounding the Daman site, maize crops are planted in the homogeneous underlying surface; the densely distributed plants tend to flatten the small rough elements, leading to few z 0m effects from wind direction. In conclusion, wind direction is a main driving factor for z 0m in mountainous areas but can be ignored when modelling z 0m over flat surfaces underlying homogenous vegetation.

The Atmospheric Stability Factor
In previous studies, most researchers tend to filter z 0m under neutral conditions to simplify the z 0m estimation models, which are built on the assumption that z 0m under neutral conditions can represent the average level and that z 0m is not significantly different under different atmospheric conditions. However, based on the results of our factor analysis, this assumption might be incorrect. As listed in Table 6, the Obukhov lengths are binned by stability class to simplify analysis. Although the definition of the stability bins varies with local environmental conditions, these definitions have previously been used by authors [37,38]. The different atmospheric stability conditions in August are counted and presented by the frequency pie chart in Figure 7. The stable condition percentage is highest at the Arou and Daman sites, while unstable conditions form the highest proportion at the Guantan site. However, the neutral condition comprises only a small portion of all sites; therefore, the z 0m value under the neutral condition is not representative, especially for instantaneous flux calculations. Table 6. Classification of atmospheric stability.

Atmospheric Condition Definition
Stable 0 < L < 1000 Unstable −1000 < L < 0 Neutral L < −1000 or L > 1000  Figure 8 shows the distribution of half-hourly z 0m in August for stable, unstable and neutral conditions. Obviously, z 0m presents clear differences under different stability conditions at each site.
Because neutral conditions accounted for the smallest percentage, there are not enough samples to ensure that z 0m is normally distributed under these conditions. A comparison of these results shows that the averaged z 0m under stable conditions has the highest value, followed by that under neutral conditions, and the lowest values appear under unstable conditions. Under stable conditions, the z 0m values have a border value range and are distributed in a nearly normal fashion. Actually, atmospheric stability is closely related to the alternation of day and night. In the daytime, affected by the solar radiation, turbulent exchanges occur frequently between vegetation and the atmosphere and unstable atmospheric conditions are prominent; therefore, z 0m maintains low levels. At night, the atmospheric turbulence weakens, and z 0m increases significantly under the more stable conditions. In addition, because the relative frequencies of z 0m are similarly distributed among these three sites, we deduce that the influence on z 0m caused by atmospheric stability has nothing to do with the underlying surface.

Conclusions
This study statistically investigated the driving factors that influence z 0m at different time scales using three AWSs in the Heihe river basin that collect data for different land surface patterns. Aerodynamic factors, including wind speed and Obukhov length, are the main driving forces for the half-hourly changes in z 0m at the Arou Site and Guantan site, contributing 44.03% and 42.87%, respectively, at those sites, while topographic factors linked to the wind deflection angle contribute 38.33% and 38.66% to the half-hourly changes in z 0m at the two sites, respectively. Regarding the daily z 0m at the Arou and Guantan sites, only wind speed and wind deflection angle are correlating factors. The driving factors for half-hourly z 0m at the Daman site include both vegetation and aerodynamic factors, with contributions of 53.14% and 30.29%, respectively. Similar results were found for the daily z 0m at Daman site, except that atmospheric stability was not considered. It should be noted that vegetation factors such as NDVI and LAI have no significant relationship with either the half-hourly z 0m or the daily z 0m at the Arou and Guantan sites, but these vegetation factors are the primary driving factors for the Daman site, probably because of fewer changes in the botanical characteristics of grass and forest over growing seasons compared with the rapid changes due to crop growth in farmland.
Researchers should recognize that meteorological factors cannot be neglected in z 0m estimation models, especially for instantaneous flux simulations.
However, there are still some uncertainties in the results. First, this study focused only on the changes of z 0m during the vegetative growth period; consequently, which factors are most influential in other seasons is still unclear. Second, the factors tested in this paper are not comprehensive due to data source limitations; therefore, more factors that may influence z 0m should be taken into account. Third, the displacement height d, as the other important roughness parameter, is discussed relatively little in this paper; and the correlation between instantaneous d and z 0m may influence the accuracy of the observed z 0m . Thus, the driving factors of d and z 0m need to be examined at more sites with a wider variety of underlying surface and climate features to improve the robustness of the conclusions.