A Method for Estimating the Aerodynamic Roughness Length with NDVI and BRDF Signatures Using Multi-Temporal Proba-V Data

Aerodynamic roughness length is an important parameter for surface fluxes estimates. This paper developed an innovative method for estimation of aerodynamic roughness length (z0m) over farmland with a new vegetation index, the Hot-darkspot Vegetation Index (HDVI). To obtain this new index, the normalized-difference hot-darkspot index (NDHD) is introduced using a semi-empirical, kernel-driven bidirectional reflectance model with multi-temporal Proba-V 300-m top-of-canopy (TOC) reflectance products. A linear relationship between HDVI and z0m was found during the crop growth period. Wind profiles data from two field automatic weather station (AWS) were used to calibrate the model: one site is in Guantao County in Hai Basin, in which double-cropping systems and crop rotations with summer maize and winter wheat are implemented; the other is in the middle reach of the Heihe River Basin from the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) project, with the main crop of spring maize. The iterative algorithm based on Monin–Obukhov similarity theory is employed to calculate the field z0m from time series. Results show that the relationship between HDVI and z0m is more pronounced than that between NDVI and z0m for spring maize at Yingke site, with an R2 value that improved from 0.636 to 0.772. At Guantao site, HDVI also exhibits better performance than NDVI, with R2 increasing from 0.630 to 0.793 for summer maize and from 0.764 to 0.790 for winter wheat. HDVI can capture the impacts of crop residue on z0m, whereas NDVI cannot.


Introduction
The aerodynamic roughness length (z 0m ) is defined as the height at which the wind speed becomes zero under neutral conditions [1].Most currently used land surface models require estimates of aerodynamic roughness length to characterize the momentum transfer between surface and atmosphere [2][3][4].
The traditional method used to calculate z 0m is based on measurements of wind profiles at different levels over the ground under neutral atmospheric conditions, by applying the profile equation derived on the basis of Monin-Obukhov similarity theory [3,5,6].Typically, according to many long-term observations, modelers have assumed that the momentum roughness is identical at all locations that have similar type of surface coverage [7] and have constructed look-up tables of z 0m values over different surfaces [8].However, these look-up approaches ignore the inherent temporal and spatial variability of a certain land use type and its concomitant effects on momentum transfer [9].Aerodynamic roughness mainly depends on the geometric features and distributions of the roughness elements [10].For vegetated surfaces in particular, the mean canopy height, the canopy structure and the plant density are key variables [11].In the past two decades, remote sensing has emerged as an effective way to retrieve surface information and parameterize aerodynamic roughness on the global or regional scale [12][13][14], and several models have been developed for z 0m as a function of vegetation physical structural parameters, such as leaf area index (LAI) [15,16], canopy area index (CAI) and frontal area index (FAI) [17,18].Besides, optical parameters such as normalized difference vegetation index (NDVI) estimated by remote sensing have been widely used for z 0m estimation [18][19][20][21].For cropland, NDVI is closely correlated with z 0m during crop growing period [18,19], Gupta et al. and Moran et al. described the relationship between z 0m and NDVI as z 0m = exp(a + b NDVI) [20,22].
Nevertheless, the vegetation indexes calculated from single observed angle have limited capacity to retrieve three-dimensional vegetation structures that are closely related to z 0m .One potential solution is to bring in multi-angular optical remote sensing [23][24][25][26].Sunlight hitting the vegetation canopy is scattered unevenly because of surface roughness, which is related to the canopy's shape and height.Multi-angular observations can capture the uneven scattering of sunlight by vegetation, which can be described by the bidirectional reflectance distribution function (BRDF) related to the biophysical structural information [24,[26][27][28].The hotspot-darkspot index (HDS) calculated from multi-angular optical remote-sensing data is capable of representing the geometric structures of vegetation [29].A high HDS index results from a strong contrast between the hotspot and darkspot reflectances, which is a manifestation of the roughness of the canopy surface [26].
Moreover, the annual variability of crop phenology and daily variability of crop development are major sources of uncertainty for z 0m assessments over vegetation.The newly launched Proba-V satellite collects daily reflectance data in the visible and near-infrared (NIR) ranges.The optical design of Proba-V consists of three cameras, providing wide-view angles of ±51 • from nadir that measure abundant directional effects in the surface reflectance.Thus, the BRDF of each pixel can be derived within a few days of observation [30].Accordingly, this paper aims to develop a new method for estimation of aerodynamic roughness length based on NDVI and BRDF parameters over farmland.The calibration of the new model is done using field-based measurements of wind profiles data from different climate systems.

Site Description
In this paper, two sites with typical climates and crop types are selected: the middle reach of the Heihe River Basin with arid climate system in northwestern China and Guantao County with semi-moisture climate system on the North China Plain, Hai basin.
(1) The middle reach of the Heihe River Basin (98 • 57 -100 • 52 E, 38 • 39 -39 • 59 N) is located in the western part of Gansu Province near Zhangye City, which is the main irrigation agriculture economic zone and main water consumption area in the Heihe River Basin. Figure 1a displays the land-cover classification results in this area determined by Zhong et al. in 2014 using multiple classifiers and multi-source remotely sensed data [31,32].In recent years, many research projects have addressed land-surface processes, hydrology and water resources, such as HEIFE [33], WATER project [34], and Hi-WATER [35].Large quantities of ground observation data from meteorological, hydrological and energy flux stations and associated study results have been accumulated, including data from the Yingke automatic weather station (AWS), which were used in this paper (Figure 1b).The Yingke AWS is located in a typical oasis in the middle reach of the Heihe River Basin on very flat terrain and is surrounded by the Gobi Desert.The observational equipment includes a 40-m high tower with the instruments facing north to measure the wind speed, wind direction, atmosphere temperature and temperature humidity at different heights [36].The main crop near the site is irrigated spring maize, which has a maximum height of approximately 1.8 m [37].47 N) and has a warm temperate, semi-humid, continental monsoon climate.Guantao is a major agricultural county, and the main cultivated crops are winter wheat, maize and cotton.Arable land accounts for 63% of the whole county's area.We created the land-cover map of Guantao County using multi-temporal and high-resolution remote sensing images (GF-1, HJ-1A/1B CCD and TM images), Support vector machine (SVM) was selected for classification with guidance regarding the local crop phonology [38,39] (Figure 2a).The AWS in Guantao was installed with support from the Global Environmental Facility (GEF) Hai Basin project in 2014 (Figure 2b), similar with Yingke site, the equipment on AWS tower includes meteorological sensors installed on different layers.This site is dominated by double crops with the rotation of winter wheat and summer maize (Figure 3).
Remote Sens. 2016, 8, 6 3 of 15 (2) Guantao County is located in the south of Hebei Province in the Hai Basin (115°06′-115°40′E, 36°27′-36°47′N) and has a warm temperate, semi-humid, continental monsoon climate.Guantao is a major agricultural county, and the main cultivated crops are winter wheat, maize and cotton.Arable land accounts for 63% of the whole county's area.We created the land-cover map of Guantao County using multi-temporal and high-resolution remote sensing images (GF-1, HJ-1A/1B CCD and TM images), Support vector machine (SVM) was selected for classification with guidance regarding the local crop phonology [38,39] (Figure 2a).The AWS in Guantao was installed with support from the Global Environmental Facility (GEF) Hai Basin project in 2014 (Figure 2b), similar with Yingke site, the equipment on AWS tower includes meteorological sensors installed on different layers.This site is dominated by double crops with the rotation of winter wheat and summer maize (Figure 3).(2) Guantao County is located in the south of Hebei Province in the Hai Basin (115°06′-115°40′E, 36°27′-36°47′N) and has a warm temperate, semi-humid, continental monsoon climate.Guantao is a major agricultural county, and the main cultivated crops are winter wheat, maize and cotton.Arable land accounts for 63% of the whole county's area.We created the land-cover map of Guantao County using multi-temporal and high-resolution remote sensing images (GF-1, HJ-1A/1B CCD and TM images), Support vector machine (SVM) was selected for classification with guidance regarding the local crop phonology [38,39] (Figure 2a).The AWS in Guantao was installed with support from the Global Environmental Facility (GEF) Hai Basin project in 2014 (Figure 2b), similar with Yingke site, the equipment on AWS tower includes meteorological sensors installed on different layers.This site is dominated by double crops with the rotation of winter wheat and summer maize (Figure 3).

In Situ Data
The AWS systems at these two sites collected data on wind speed, wind direction, air temperature and air humidity at 10-min intervals (Table 1).To guarantee the reliability of the wind profile results, the raw datasets were selected to obtain highly accurate aerodynamic roughness values based on the following criteria: (1) the wind speed is greater than 1 m/s; (2) the wind friction velocity u * exceeds 0.2 m/s; and (3) data obtained on rainy days are discarded.For further information about the operating principles of the instruments and the data processing, refer to Liu and coworkers' study [35,36].

Satellite Data
The Proba-V satellite was launched on 6 May 2013 and was designed to bridge the gap in space-borne vegetation measurements between the SPOT-VGT and the upcoming Sentinel-3 satellites.The technical specifications of the on-board sensors on the Proba-V satellite are listed in Tables 2 and 3 [40].In this study, 300 m daily top-of-canopy (TOC) reflectance products were downloaded from the web page of VITO's Product Distribution Portal (PDP) in HDF5 file format (http://www.vito-eodata.be/).

In Situ Data
The AWS systems at these two sites collected data on wind speed, wind direction, air temperature and air humidity at 10-min intervals (Table 1).To guarantee the reliability of the wind profile results, the raw datasets were selected to obtain highly accurate aerodynamic roughness values based on the following criteria: (1) the wind speed is greater than 1 m/s; (2) the wind friction velocity u * exceeds 0.2 m/s; and (3) data obtained on rainy days are discarded.For further information about the operating principles of the instruments and the data processing, refer to Liu and coworkers' study [35,36].

Satellite Data
The Proba-V satellite was launched on 6 May 2013 and was designed to bridge the gap in space-borne vegetation measurements between the SPOT-VGT and the upcoming Sentinel-3 satellites.The technical specifications of the on-board sensors on the Proba-V satellite are listed in Tables 2 and 3 [40].In this study, 300 m daily top-of-canopy (TOC) reflectance products were downloaded from the web page of VITO's Product Distribution Portal (PDP) in HDF5 file format (http://www.vito-eodata.be/).The Proba-V 300 m TOC products were acquired for 2014 in the middle reach of the Heihe River Basin and for 2015 in Guantao based on the crop phenology (Figure 3).Two images (Titles: X27Y03 and X28Y03) were needed to cover the middle reach of the Heihe River Basin, whereas one was sufficient for Guantao (Titles: X29Y03).The required data groups including the solar zenith angles (SZA), solar azimuth angles (SAA), viewing zenith angles of the visible and NIR (VNIR) detector (VZA), viewing azimuth angles of the VNIR detector (VAA), quality control, and the TOC reflectances of RED and NIR were extracted and converted to image files from the original HDF5 files by Spirits software [41].The images were mosaicked and clipped to the study area, projected to Albers conical equal area projections, and resampled using the bilinear interpolation method.

Ground Aerodynamic Roughness Length
The aerodynamic surface roughness length can be determined iteratively based on wind profile data.Using Monin-Obukhov similarity, roughness length and zero-plane displacement can be related via the logarithmic wind profile equation [42]: where u and θ are the wind speed and potential air temperature, respectively, at 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; θ 0 is the potential temperature near the surface; L is a function of the friction velocity, friction temperature and temperature (called the Monin-Obukhov length); and Ψ m and Ψ h are the stability functions.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 [43].
For Z/L < 0 (unstable conditions) [44], For Z/L > 0 (stable conditions) [45], In this study, we use least-square method to determine z 0m [46, 47].Equation ( 1) can be written as the following simplified form: where a = u * /k, x = ln(z − d) − Ψ m and b = −lnz 0m •a.Commonly, d has a relationship with the vegetation height (h) that d equals to 0.67 h over surfaces covered by dense and homogeneous vegetation [48], therefore, d is determined through iterative algorithm which increases from 0.1 m to 3 m at intervals of 0.1 m.For each given d, there is a correlation coefficient for the fitting Equation ( 7), and d is chosen as the optimum value when the correlation coefficient reaches maximum, then, z 0m can be calculated according to d value.

BRDF Parameters with Ross-Li Model
The semi-empirical, kernel-driven BRDF model is generally used to correct BRDF effects and retrieve surface albedo from multi-angle datasets [30,49,50].This model relies on the weighted sum of an isotropic parameter and two kernels of viewing and the illumination geometry to determine the reflectance of certain observed angles [27].The Ross-Li-Maignan model is formulated as follows.
In Equation ( 8), the surface reflectance (R) is expressed as a function of three components.The relative azimuth is ϕ = ϕ i − ϕ r .f iso is the reflectance acquired by nadir observation when the solar zenith angle is zero.K vol describes the volume scattering kernel, and K geo describes the surface scattering kernel.Volume scattering is caused by a horizontal layer of randomly distributed leaves, and surface scattering is caused by the shadows of natural objects.These two kernels are functions of only the solar and sensor geometries, including the solar zenith angle (θ i ), view zenith angle (θ r ) and relative azimuth angle (ϕ).f vol (λ) and f geo (λ) are the spectrally dependent BRDF kernel coefficients.In this study, K vol and K geo were calculated using the Li-Sparse model and the Ross Thick model, as described by Equations ( 9) and (11), respectively [28], and ξ is the phase angle calculated by Equation ( 12).
To calculate the kernel coefficients, a multiple linear regression fit was applied with 21-day series of reflectance and angle data.Because the model is relatively insensitive to noisy data [49], here, we use Proba-V's quality control file to determine each pixel's quality flag and, thereby, discard the cloudy and invalid pixels.For a specific pixel, if at least five cloud-free observations of the surface are available during a 21-day period, the kernel coefficients are calculated; otherwise, it is assigned an invalid value.Then, reflectance at any angle can be acquired in the hemisphere.

NDHD, NDVI and HDVI
In the principal plane, the reflectance exhibits large variations at different scattering angles.This information is very useful for measuring the canopy structure.The largest reflectance in the principal plane is chosen as the hot spot, and the smallest defined as the dark spot where shadows created by the vegetation canopy are maximally observed [24].In practice, usually we cannot obtain both hotspot and darkspot directly on the simulated curve as shown in Figure 4, so the darkspot (hotspot) is defined at the same view zenith angle with the hotspot (darkspot) but in the forward-scattering (back-scattering) direction, in order to extract the hotspot reflectance (ρ HS ) and the darkspot reflectance (ρ DS ).As Leblanc et al. proposed in 2001 [51], the normalized-difference hotspot-darkspot index (NDHD) is the normalized difference between ρ HS and the ρ DS : The normalized index selected for this research is to reduce the influence of canopy optical properties and enhance the importance of canopy structure.Analogously, the HDS index ((ρ HS − ρ DS )/ρ DS ), which was proposed by Lacaze et al. [29], has been widely used in recent research.
For HDS, the problem of mixed pixels becomes prominent at Proba-V's 300 m resolution: Some pixels containing mixed cropland and bare land demonstrate low ρ DS values, resulting in abnormal HDS values.Therefore, NDHD was applied here.
Daily NDVI images are generated from Proba-V's RED and NIR bands.Based on the quality control file, time series of 5-day composite NDVI at 300 m resolution are acquired to reduce the disturbance of cloudy pixels.Unlike the reflectance, the use of NDVI partially avoids the impact of BRDF signatures because the directional effects are similar in the visible and NIR bands, and therefore, the effect on the individual reflectances can be reduced by taking the ratio of the two bands [52].However, BRDF signatures, such as the NDHD, can express z 0m according to the point on the canopy structure, whereas NDVI cannot.Because both NDVI and NDHD are positively correlated with z 0m but their correlations rely on different aspects, a new vegetation index, the Hot-darkspot Vegetation Index (HDVI), which combines NDVI and NDHD is proposed and defined by Equation (14).
A linear relationship between HDVI and z 0m was attempted to map z 0m as follows.

Simulated Reflectance in Red and NIR Band
As shown in Figure 4, the pixel at which the Yingke station is located is taken as an example to demonstrate the simulated changes in the reflectances of the red and NIR bands in Proba-V's field of view.The solar zenith angle is fixed at 35 degrees in the principal plane.Based on the calendar of Yingke's maize, the first day of each crop-growing month is selected for this example.ρ )/ρ ), which was proposed by Lacaze et al. [29], has been widely used in recent research.For HDS, the problem of mixed pixels becomes prominent at Proba-V's 300 m resolution: Some pixels containing mixed cropland and bare land demonstrate low ρ values, resulting in abnormal HDS values.Therefore, NDHD was applied here.Daily NDVI images are generated from Proba-V's RED and NIR bands.Based on the quality control file, time series of 5-day composite NDVI at 300 m resolution are acquired to reduce the disturbance of cloudy pixels.Unlike the reflectance, the use of NDVI partially avoids the impact of BRDF signatures because the directional effects are similar in the visible and NIR bands, and therefore, the effect on the individual reflectances can be reduced by taking the ratio of the two bands [52].However, BRDF signatures, such as the NDHD, can express z0m according to the point on the canopy structure, whereas NDVI cannot.Because both NDVI and NDHD are positively correlated with z0m but their correlations rely on different aspects, a new vegetation index, the Hot-darkspot Vegetation Index (HDVI), which combines NDVI and NDHD is proposed and defined by Equation (14).
A linear relationship between HDVI and z0m was attempted to map z0m as follows.

Simulated Reflectance in Red and NIR Band
As shown in Figure 4, the pixel at which the Yingke station is located is taken as an example to demonstrate the simulated changes in the reflectances of the red and NIR bands in Proba-V's field of view.The solar zenith angle is fixed at 35 degrees in the principal plane.Based on the calendar of Yingke's maize, the first day of each crop-growing month is selected for this example.This figure shows that the reflectances of the two channels exhibit a decreasing trend from the back-scattering region to the forward-scattering region.The lowest point, which is known as the darkspot, can be identified on most curves as the point at which the view zenith angle is the same as the solar zenith angle.The hotspot did not appear in the simulated range of the view zenith angle, and the darkspot did not appear on the red-reflectance curves on 1 August and 1 September.
The NDHD of the NIR reflectance varies regularly with the crop growth period, with high values during the initial and terminal stages of crop growth and low values during the vigorous growth stage, in response to the variation of the crop-canopy structures.In contrast, the NDHD of This figure shows that the reflectances of the two channels exhibit a decreasing trend from the back-scattering region to the forward-scattering region.The lowest point, which is known as the darkspot, can be identified on most curves as the point at which the view zenith angle is the same as the solar zenith angle.The hotspot did not appear in the simulated range of the view zenith angle, and the darkspot did not appear on the red-reflectance curves on 1 August and 1 September.
The NDHD of the NIR reflectance varies regularly with the crop growth period, with high values during the initial and terminal stages of crop growth and low values during the vigorous growth stage, in response to the variation of the crop-canopy structures.In contrast, the NDHD of RED reflectance exhibits less variation (the nearly identical inclination of the five curves in Figure 4b) during the crop growth period, and thus, we calculate the NDHD by applying the NIR band.

Relationship between NDVI/HDVI and z 0m
As shown in Figure 3, the growing season of spring maize in the Yingke Site is from May to October.In this paper, daily averaged aerodynamic roughness values are calculated from AWS raw data with 10-min intervals.To match the NDVI data, the mean value of z 0m is obtained over five days with the assumption that z 0m changed little during that period.
For the Yingke site, the temporal variations of the aerodynamic roughness, as deduced from the values for each five-day period, reflect the process of crop growth, which exhibits a characteristic rise and fall according to the crop growth cycle (Figure 5).The NDVI values are low at the beginning of crop emergence and increase rapidly over the course of the subsequent month.The peak growth stage lasts 184-244 days.Then, the NDVI values decrease as the spring maize approaches maturity.Clearly, the variations of NDVI are in good agreement with those of z 0m , suggesting that NDVI is a reasonable predictor of z 0m .Compared with the NDVI profile during the maize growth period, the changing trend of NDHD generally exhibited almost opposite directions, initially increasing rapidly and maintaining a high level from Days 264 to 299.RED reflectance exhibits less variation (the nearly identical inclination of the five curves in Figure 4b) during the crop growth period, and thus, we calculate the NDHD by applying the NIR band.

Relationship between NDVI/HDVI and z0m
As shown in Figure 3, the growing season of spring maize in the Yingke Site is from May to October.In this paper, daily averaged aerodynamic roughness values are calculated from AWS raw data with 10-min intervals.To match the NDVI data, the mean value of z0m is obtained over five days with the assumption that z0m changed little during that period.
For the Yingke site, the temporal variations of the aerodynamic roughness, as deduced from the values for each five-day period, reflect the process of crop growth, which exhibits a characteristic rise and fall according to the crop growth cycle (Figure 5).The NDVI values are low at the beginning of crop emergence and increase rapidly over the course of the subsequent month.The peak growth stage lasts 184-244 days.Then, the NDVI values decrease as the spring maize approaches maturity.Clearly, the variations of NDVI are in good agreement with those of z0m, suggesting that NDVI is a reasonable predictor of z0m.Compared with the NDVI profile during the maize growth period, the changing trend of NDHD generally exhibited almost opposite directions, initially increasing rapidly and maintaining a high level from Days 264 to 299.15) and ( 16), respectively.
Figure 6 reports the field-observed aerodynamic roughness (z0m), plotted as linear functions of NDVI and HDVI in the Yingke Site.This figure clearly demonstrates that the NDVI of arable land is correlated with z0m (R 2 = 0.636, n = 33).In contrast, the data are less dispersed and the correlation is higher (R 2 = 0.793, n = 33) when z0m is plotted as a function of HDVI (Figure 6b).In addition, the NDVI and HDVI values in October (points from A~F in Figure 6) exhibit obviously different correlations with z0m: For the same value of z0m, NDVI showed much greater deviation from the fitted curve, while HDVI was distributed around the fitted curve for the six marked points.This finding indicates that NDVI decreased quickly with z0m during the crop harvesting stage when the leaf cells of spring maize senesce, and the leaves become yellow and withered, resulting in z0m less evaluated if only the changes of NDVI were taken into concern.In contrast, HDVI, which contains information on the crop-planting structure, has advantages for retrieving z0m in the late growth period.z 0m is presented as a five-day average, and NDVI is shown as the maximum value over five days.NDHD and HDVI were calculated with Equations ( 15) and ( 16), respectively.
Figure 6 reports the field-observed aerodynamic roughness (z 0m ), plotted as linear functions of NDVI and HDVI in the Yingke Site.This figure clearly demonstrates that the NDVI of arable land is correlated with z 0m (R 2 = 0.636, n = 33).In contrast, the data are less dispersed and the correlation is higher (R 2 = 0.793, n = 33) when z 0m is plotted as a function of HDVI (Figure 6b).In addition, the NDVI and HDVI values in October (points from A~F in Figure 6) exhibit obviously different correlations with z 0m : For the same value of z 0m , NDVI showed much greater deviation from the fitted curve, while HDVI was distributed around the fitted curve for the six marked points.This finding indicates that NDVI decreased quickly with z 0m during the crop harvesting stage when the leaf cells of spring maize senesce, and the leaves become yellow and withered, resulting in z 0m less evaluated if only the changes of NDVI were taken into concern.In contrast, HDVI, which contains information on the crop-planting structure, has advantages for retrieving z 0m in the late growth period.
The same method was used to determine the differences between NDVI and HDVI for winter wheat and summer maize at the Guantao site.The linear fitting equations based on the crop calendar of Guantao are shown in Figure 7.The R 2 values of the four fitted equations also indicate that HDVI is more significantly correlated with z 0m than NDVI.Comparable results are also obtained for the relationships of z 0m of Guantao's summer maize and NDVI/HDVI, with the correlation coefficient improving from 0.670 to 0.793.The HDVI-based model produces a slightly better fitting result than the NDVI-based model (the R 2 value ranges from 0.764 to 0.790) for winter wheat in Guantao.The same method was used to determine the differences between NDVI and HDVI for winter wheat and summer maize at the Guantao site.The linear fitting equations based on the crop calendar of Guantao are shown in Figure 7.The R 2 values of the four fitted equations also indicate that HDVI is more significantly correlated with z0m than NDVI.Comparable results are also obtained for the relationships of z0m of Guantao's summer maize and NDVI/HDVI, with the correlation coefficient improving from 0.670 to 0.793.The HDVI-based model produces a slightly better fitting result than the NDVI-based model (the R 2 value ranges from 0.764 to 0.790) for winter wheat in Guantao.Coefficients a and b of Equation ( 14) are acquired for different crop types from the linear fitting results shown in Figure 6 and Figure 7 at the Yingke site and the Guantao site, respectively.The results are shown in Table 4, along with the coefficients of determination (R 2 ).The robustness of z0m The same method was used to determine the differences between NDVI and HDVI for winter wheat and summer maize at the Guantao site.The linear fitting equations based on the crop calendar of Guantao are shown in Figure 7.The R 2 values of the four fitted equations also indicate that HDVI is more significantly correlated with z0m than NDVI.Comparable results are also obtained for the relationships of z0m of Guantao's summer maize and NDVI/HDVI, with the correlation coefficient improving from 0.670 to 0.793.The HDVI-based model produces a slightly better fitting result than the NDVI-based model (the R 2 value ranges from 0.764 to 0.790) for winter wheat in Guantao.Coefficients a and b of Equation ( 14) are acquired for different crop types from the linear fitting results shown in Figure 6 and Figure 7 at the Yingke site and the Guantao site, respectively.The results are shown in Table 4, along with the coefficients of determination (R 2 ).The robustness of z0m model performance is evaluated using the root mean square error (RMSE) and the mean absolute Coefficients a and b of Equation ( 14) are acquired for different crop types from the linear fitting results shown in Figures 6 and 7 at the Yingke site and the Guantao site, respectively.The results are shown in Table 4, along with the coefficients of determination (R 2 ).The robustness of z 0m model performance is evaluated using the root mean square error (RMSE) and the mean absolute error (MAE), which can measure the accuracy of the established model.The RMSE and MAE values between the observed z 0m and estimated z 0m for different crops range from 0.024 to 0.035 and from 0.020 to 0.028, respectively, which are lower than 0.1, indicating consistency between the observed z 0m and HDVI-based estimated z 0m .In particular, Durbin-Watson statistics are generally lower than 2, indicating that the error terms are positively auto-correlated.Results of the F-test suggest that the fitted models with p-values lower than 0.01, indicating a significant relationship between HDVI/NDVI data and field observed z 0m at the 99% confidence level for all crops, confirming the strong capability of the model.We use the relationship between HDVI and the aerodynamic roughness length z 0m to derive maps of the roughness length.To explore the seasonal and spatial variability of the aerodynamic roughness, we focus on the middle reach of the Heihe River Basin and estimate z 0m over the spring maize mask [31,32] of certain days from May to October, as shown in Figure 8.The value of z 0m ranges between 0 and 0.25 throughout the spring maize growing season in this area, reaches a peak in July and then declines beginning in August when the height of the maize plants reaches a plateau.The largest standard deviation value of all maize pixels occurred in June, indicating that z 0m exhibited more significant spatial differentiation during this period.This is most likely attributable to spatial differences in the sowing time, planting density and fertilization method of spring maize.The differences are also amplified by the rapid growth that occurs during this stage.

Discussion
With the development of remote-sensing technology, satellite-based algorithms are now routinely applied to retrieve terrestrial parameters, such as aerodynamic roughness.Because vegetation height is essential in most aerodynamic roughness models, usually field investigation or indirect methods are needed to acquire the spatial distribution of the vegetation height.Here, we try to wean the vegetation height observation off by developing a purely remote-sensing-based model for aerodynamic roughness estimation over arable land, improving the timeliness put forward for aerodynamic roughness.The introduction of BRDF signatures improves the accuracy and persuasiveness of the model because the multi-angle remote-sensing data have incomparable advantages for describing the crop canopy structure.
Because crops growing on farmland are flat and uniformly distributed, the topographic changes can be ignored, and the aerodynamic roughness can be separated into two parts: vertical roughness associated with crop growth conditions and horizontal roughness associated with crop planting structures.NDVI still plays the leading role in reflecting crop growth conditions, and NDHD expresses the variations in the vegetation spatial structure.Satisfactory results for aerodynamic roughness length estimation were obtained.According to our results, HDVI clearly exhibits better performance for aerodynamic roughness estimation than NDVI.Compared with NDVI, HDVI's range is wider, and it can exceed 1.The exponential algorithm describing z 0m [21,23] can be improved with linear model after the introduction of NDHD.
According to the results shown in Figures 6 and 7, this method is more efficient for maize than wheat.This finding can be explained by two reasons: (1) Wheat plants grow to approximately the same height and are usually closely spaced.Thus, the canopy of wheat is more flat and homogeneous than that of maize, and the spatial heterogeneity of wheat land, which affects the NDHD, changes less during the crop growth period.(2) The differences in the crop ripening stage and harvest period for the two crops may affect the results.Typically, maize stalk residues remain quite tall and are sparsely distributed over farmland after harvest (Figure 9a), as reflected by higher NDHD values.However, NDVI is difficult to interpret at this stage.For wheat, the approximate height of the residue is only 30 cm (Figure 9b).Therefore, the canopy structure is similar to that of bare land, and few BRDF signatures are captured by the NDHD.Thus, the plant structures and growth characteristics of different crops should be considered when modeling z 0m .
The method discussed here only applies to two crops (wheat and maize), because the aerodynamic roughness of different land use types differ significantly (by as much as orders of magnitude) [53,54].Results are based on limited datasets under particular meteorological conditions.Further research is needed to determine the model's applicability to other crop and vegetation types.This study was limited by inadequate AWS data for different underlying surface types.Although eddy covariance (EC) and large-aperture scintillometer (LAS) measurements have been used to acquire local roughness length values in previous studies [55,56], the estimation of z 0m using different instruments may introduce deviations in the same region because of the differences among AWS, EC and LAS footprints.The consistency and precision of ground data constitute the foundation of modeling based on remote sensing data.In our study, the two AWS sites were located in homogeneous and flat arable land; thus, the observed z 0m values are representative of the corresponding remote-sensing pixels, and the impact of topographic fluctuation is avoided.Actually, the surface condition of the vegetation cover is not the only factor that must be considered.The near-ground wind speed, wind direction and stability conditions also strongly influence the aerodynamic roughness [55,57].In our study, meteorological conditions were not accounted for, which may have affected the fitting results for the field-observed z 0m and HDVI.The effects of wind speed and wind direction on z 0m should be studied further, moreover, boundary layer conditions also play a role in surface-atmosphere exchanges parameterized through z 0m .It should be noted, however, that high-resolution surface meteorological data acquisition is also a key problem in determining the spatial distribution of z 0m .
changes can be ignored, and the aerodynamic roughness can be separated into two parts: vertical roughness associated with crop growth conditions and horizontal roughness associated with crop planting structures.NDVI still plays the leading role in reflecting crop growth conditions, and NDHD expresses the variations in the vegetation spatial structure.Satisfactory results for aerodynamic roughness length estimation were obtained.According to our results, HDVI clearly exhibits better performance for aerodynamic roughness estimation than NDVI.Compared with NDVI, HDVI's range is wider, and it can exceed 1.The exponential algorithm describing z0m [21,23] can be improved with linear model after the introduction of NDHD.
According to the results shown in Figures 6 and 7, this method is more efficient for maize than wheat.This finding can be explained by two reasons: (1) Wheat plants grow to approximately the same height and are usually closely spaced.Thus, the canopy of wheat is more flat and homogeneous than that of maize, and the spatial heterogeneity of wheat land, which affects the NDHD, changes less during the crop growth period.(2) The differences in the crop ripening stage and harvest period for the two crops may affect the results.Typically, maize stalk residues remain quite tall and are sparsely distributed over farmland after harvest (Figure 9a), as reflected by higher NDHD values.However, NDVI is difficult to interpret at this stage.For wheat, the approximate height of the residue is only 30 cm (Figure 9b).Therefore, the canopy structure is similar to that of bare land, and few BRDF signatures are captured by the NDHD.Thus, the plant structures and growth characteristics of different crops should be considered when modeling z0m.The method discussed here only applies to two crops (wheat and maize), because the aerodynamic roughness of different land use types differ significantly (by as much as orders of magnitude) [53,54].Results are based on limited datasets under particular meteorological Proba-V was designed to offer global coverage at spatial resolutions of 100 m, 300 m and 1 km [40].The 300 m TOC reflectance products used here can facilitate the observation of a certain terrestrial target on a daily basis in most locations on Earth.Similar to the angular observations collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument [49], multi-temporal data are extracted from a series of Proba-V passes.The BRDF model parameters were derived from the observed time series.Regression of the semi-empirical kernel-driven BRDF model should be performed over sub-periods that are short enough that the surface can be considered time invariant and long enough to contain sufficient sensor data for the regression [27].Obtaining sufficient cloud-free surface reflectance data is challenging over certain regions.The main limitation of this method is the assumption that the target is stable during the synthesis period.This assumption is obviously invalid for vegetation-covered land.Thus, we must assume that the variations in the BRDF shape are limited, so the period of 21 days was set for regression with multi-temporal Proba-V images to ensure the accuracy of the BRDF signatures.Consequently, multi-angular optical remote-sensing platforms with high spatial resolution and short revisit periods are still needed to further advance this research.Additionally, in the future, attempts to acquire multi-angular data from multi-source remote sensing data should be made.

Conclusions
This paper proposed a new model for estimating z 0m over arable land with a new vegetation index, the Hot-darkspot Vegetation Index (HDVI) with multi-temporal Proba-V 300-m TOC reflectance products.The linear relationship between HDVI and z 0m was found for the crop growth period.The results show that the relationship between HDVI and z 0m is more pronounced than that between NDVI and z 0m for spring maize at the Yingke site, with an R 2 value that improved from 0.636 to 0.772.At the Guantao site, in which double-cropping systems and crop rotations with summer maize and winter wheat are implemented, HDVI also exhibits better performance than NDVI, with R 2 increasing from 0.630 to 0.793 for summer maize and from 0.764 to 0.790 for winter wheat.The differences probably originate from the crop ripening stage and harvest period, indicating crop residue's impacts on z 0m captured by NDHD, when NDVI makes little sense to crop at that stage.NDHD from near-infrared reflectance contains more meaningful BRDF information associated with canopy structure against red reflectance.However, this study was based on limited available data, HDVI will be adopted and qualified over different vegetation types in long periods of time with the accumulation of robust validation exercise.

Figure 1 .
Figure 1.(a) Land cover map of the middle reach of the Heihe River Basin in 2014; and (b) photograph of the automatic weather station (AWS) tower in Yingke.

Figure 2 .
Figure 2. (a) Land cover map of Guantao County in 2015; and (b) photograph of the AWS tower in Guantao.

Figure 1 .
Figure 1.(a) Land cover map of the middle reach of the Heihe River Basin in 2014; and (b) photograph of the automatic weather station (AWS) tower in Yingke.

Figure 1 .
Figure 1.(a) Land cover map of the middle reach of the Heihe River Basin in 2014; and (b) photograph of the automatic weather station (AWS) tower in Yingke.

Figure 2 .
Figure 2. (a) Land cover map of Guantao County in 2015; and (b) photograph of the AWS tower in Guantao.

Figure 2 .
Figure 2. (a) Land cover map of Guantao County in 2015; and (b) photograph of the AWS tower in Guantao.

Figure 3 .
Figure 3. Crop calendars for the Yingke and Guantao Experimental Areas.

Figure 3 .
Figure 3. Crop calendars for the Yingke and Guantao Experimental Areas.

Figure 4 .
Figure 4. Simulated: (a) near-infrared (NIR); and (b) RED reflectances at Yingke station according to changes in the view zenith angle on the first day of the month, for five months.

Figure 4 .
Figure 4. Simulated: (a) near-infrared (NIR); and (b) RED reflectances at Yingke station according to changes in the view zenith angle on the first day of the month, for five months.

Figure 5 .
Figure 5.Time series profiles of NDVI, HDVI, NDHD and field observed z0m during the spring maize growing season at Yingke site in 2014.z0m is presented as a five-day average, and NDVI is shown as the maximum value over five days.NDHD and HDVI were calculated with Equations (15) and (16), respectively.

Figure 5 .
Figure 5.Time series profiles of NDVI, HDVI, NDHD and field observed z 0m during the spring maize growing season at Yingke site in 2014.z0m is presented as a five-day average, and NDVI is shown as the maximum value over five days.NDHD and HDVI were calculated with Equations (15) and (16), respectively.

Figure 6 .
Figure 6.Observed aerodynamic roughness length as functions of: (a) NDVI; and (b) HDVI.The data were collected during the growth period of spring maize (April to October) in 2014 at the Yingke site.Points A-F represent the values for October 1, 6, 11, 16, 21 and 26, respectively.

Figure 7 .
Figure 7. Observed aerodynamic roughness length as a function of: (a) the NDVI of summer maize; (b) the HDVI of summer maize; (c) the NDVI of winter wheat; and (d) the HDVI of winter wheat.The data were selected according to the growth period of the specific crop at the Guantao site.

Figure 6 .
Figure 6.Observed aerodynamic roughness length as functions of: (a) NDVI; and (b) HDVI.The data were collected during the growth period of spring maize (April to October) in 2014 at the Yingke site.Points A-F represent the values for October 1, 6, 11, 16, 21 and 26, respectively.

Figure 6 .
Figure 6.Observed aerodynamic roughness length as functions of: (a) NDVI; and (b) HDVI.The data were collected during the growth period of spring maize (April to October) in 2014 at the Yingke site.Points A-F represent the values for October 1, 6, 11, 16, 21 and 26, respectively.

Figure 7 .
Figure 7. Observed aerodynamic roughness length as a function of: (a) the NDVI of summer maize; (b) the HDVI of summer maize; (c) the NDVI of winter wheat; and (d) the HDVI of winter wheat.The data were selected according to the growth period of the specific crop at the Guantao site.

Figure 7 .
Figure 7. Observed aerodynamic roughness length as a function of: (a) the NDVI of summer maize; (b) the HDVI of summer maize; (c) the NDVI of winter wheat; and (d) the HDVI of winter wheat.The data were selected according to the growth period of the specific crop at the Guantao site.

Table 1 .
Details about the observation sites.

Table 3 .
Characteristics of Proba-V spectral bands.

Table 1 .
Details about the observation sites.

Table 3 .
Characteristics of Proba-V spectral bands.

Table 4 .
The coefficients used to calculate z 0m for different crop types and results.