Aerodynamic roughness length estimation

Introduction Conclusions References


Introduction
Roughness length of land surfaces is an essential variable for the parameterisation of momentum and heat exchanges.The growing interest in the estimation of the surface energy balance components from passive remote sensing leads to an increasing development of models e.g.(Bastiaanssen et al., 1998;Roerink et al., 2000;Su, 2002;Colin et al., 2006b), some of which propose detailed parameterisation of resistances to heat transfer using advanced algorithms to retrieve roughness length for heat (z 0h ) from kB −1 formulations (Massman, 1997;Blümel, 1999).However, as complex as the parameterisation can be, the actual benefit from such formulations depends on an adequate estimate of the roughness length for momentum (z 0m ).Numerous formulations to derive this parameter from NDVI can be found in many studies e.g.(Moran, 1990;Bastiaanssen, 1995), but are commonly used out of recommended bounds and on highly heterogeneous land surfaces, sometimes leading to a significant degradation of turbulent flux estimates (Colin et al., 2006a).These approaches would benefit from the combined use of passive remote sensing and land surface structure measurements from Light Detection And Ranging (LIDAR) techniques.Since the very early use of laser altimetry (Ketchum Jr., 1971), sensor performances have significantly improved, allowing airborne profiler to be used to derive the roughness of the surface (Menenti and Ritchie, 1994).More recently, satellite and airborne imaging LIDAR systems have paved the way to the mapping of vegetation properties over forest areas (Hofton et al., 2002), sometimes associated with complex topography (Dorren et al., 2007), but also on low vegetation like salt-marshes (Wang et al., 2009) or semi-arid steppes (Streutker and Glenn, 2006).
The objective of this paper is to explore the use of imaging LIDAR measurements for the estimation of the aerodynamic roughness length over a heterogeneous landscape of the Heihe river basin, a typical inland river basin in the northwest of China.This investigation is part of the Watershed Allied Telemetry Experimental Research (WATER) project (Li et al., 2009), which is a simultaneous airborne, satelliteborne, and ground-based remote sensing experiment aiming at improving the observability, understanding, and predictability of hydrological and related ecological processes at a catchment scale.LIDAR points were used to extract a Digital Surface Model (DSM) and a Digital Elevation Model (DEM) from a seamless mosaic of three flight passes over an irrigated area covered by field crops, small trees arrays and tree hedges, with a ground resolution of 1 m and a total surface of 7.2 km 2 .As a first step, the DSM is used to estimate the plan surface density and frontal surface density of obstacles to wind flow and compute a displacement height and roughness length following the work done by (Raupach, 1994) and (MacDonald et al., 1998).In a second step, both the DSM and DEM are introduced in a Computational Fluid Dynamics model (CFD) to calculate wind fields from the surface to the top of the Planetary Boundary Layer (PBL), and invert wind profiles for each calculation grid and compute a roughness length.Examples of the use of these three approaches are presented for various wind directions together with a cross-comparison of results on heterogeneous land cover and complex roughness element structures.

Theoretical background
The wind velocity profile over the land surface and in neutral atmospheric stratification conditions is commonly approximated by a simple logarithmic expression of the form: where u * is the friction velocity, k the von Karman constant, d 0 the displacement height and z 0m the aerodynamic roughness length.The later is usually expressed as a constant ratio of the canopy height for homogeneous surfaces like continuous low vegetation canopies, with a consensus for values of around z 0m h v ≈ 0.1 (Brutsaert, 1982).However, the homogeneity assumption is generally never met.Therefore, such kind of approximation is of limited interest for most of the environmental studies.It has long been demonstrated from field work and wind tunnel experiments that the drag affecting the airflow over a heterogeneous land surface is related to roughness elements density and size (Counehan, 1971;Wooding et al., 1973).This was expressed in the formulation proposed by (Lettau, 1969): where h is an effective averaged obstacle height, and λ f the frontal area index defined as: The frontal area index expresses the ratio of the frontal surface A f (perpendicular to the flow) over the total surface covered by roughness elements A T .A well-known formulation based on the combined use of h and λ f was proposed in (Raupach, 1994) to calculate the displacement height d 0 and the roughness length z 0m .Raupach's formulation of the displacement height is: and the formulation of the roughness length is: where ψ h expresses the influence of the roughness sublayer, C s is the drag coefficient for an obstacle free surface, C R the drag coefficient for an isolated obstacle, and C dl a free parameter (Raupach, 1994).In this study, we used values recommended by Raupach, i.e. the values of ψ h = 0.193, C s = 0.003, C R = 0.3, C dl = 7.5 and u * U max = 0.3 These analytic expressions of the displacement height and of the roughness length normalized by the averaged obstacle height were derived from numerous field and wind tunnel experiments, and proved to provide a good fit with experimental data (Raupach, 1994).(Theurer, 1993), quoted in (MacDonald et al., 1998), noted that z 0m and d 0 could be approached by combining the frontal area index with the plan area index defined as: where A p is the plan surface of the roughness elements within the same total surface A T .The plan area index λ p is related to the importance of intervening spaces between roughness elements.For an array of roughness elements of equivalent height, an increase of the plan area index will lead to an increase of the displacement height, and a decrease of the roughness of the obstacle array.When the plan area index tends to 1, the top surfaces of the obstacles, forming an obstacle canopy, tend to form a homogeneous surface with a very limited resistance to airflow.This explains the nonmonotonic behaviour of z 0m with λ f .For a given value of the frontal area index, an increase of λ p leads to a decrease of the drag effect of the roughness elements.Therefore the Lettau's formulation of z 0m is known to fail for plan area index values higher than 0.2-0.3 because of mutual effects of high frontal area index and limited intervening spaces (MacDonald et al., 1998).This was expressed by (MacDonald et al., 1998), who proposed two formulations for z 0m and d 0 based on Lettau's Hydrol. Earth Syst. Sci., 14, 2661-2669, 2010 www.hydrol-earth-syst-sci.net/14/2661/2010/ concept to account for a larger variety of geometrical configurations of roughness elements, and show an appropriate behaviour over the entire range of density indexes.The ratio of the displacement height over the roughness element height is expressed as: The convexity can be controlled by α.Experiments lead to recommend a value of α = 4.43 for staggered arrays of roughness elements and α = 3.59 for squared arrays (MacDonald et al., 1998).This ratio is then incorporated in the calculation of the ratio of the roughness length over the roughness element height following: The expression includes the obstacle drag coefficient C D = 1.2, and an extra β coefficient to best fit the relation with experiments.In the following study this coefficient is not used (β = 1).This formulation proved to reproduce the peak of z 0m h v for λ f = 0.15-0.30,which is consistent with wind tunnel experiments.
Beside the use of the plan area and frontal area indexes, the direct use of both the DEM and DSM in a Computational Fluid Dynamics (CFD) solver is explored.The CFD solver called Canyon, embedded in the WindStation software (Lopes, 2003), allows for numerical simulations of turbulent flows over complex topography, and can account for the geometry of surface roughness elements through the Digital Surface Model, as obtained from LIDAR data.The solver follows a control-volume approach, and solves for mass conservation, momentum conservation following Navier-Stokes equations, and also energy conservation for non-neutral situations.3-D wind fields obtained in output of the CFD express the combined effect of topography and roughness elements on the airflow, and result from the solving of the transport equation.Values of wind speed of a given profile not only characterise local effects of the vegetation structure, but the total surface stress resulting from the upstream roughness elements on a distance called the length scale (Menenti and Ritchie, 1994).This length scale is usually considered to be of 1-2 order of magnitude of the height of the measurement of the wind speed.Therefore an aerodynamic roughness length can be obtained from the wind profile of each computation grid by inverting Eq. ( 1) with values within the ground and a given elevation.As mentioned earlier, this is only true in neutral atmospheric stability conditions.

Study area
The HeiHe River Basin is a typical inland river basin in the northwest of China.Second largest inland river basin of the country, it is located between 97 • 24 -102 • 10 E and 37 • 41 -42 • 42 N, and covers an area of approximately 130 000 km 2 (Fig. 1).Experiments conducted in the scope of the WA-TER project consisted in simultaneous airborne, satelliteborne and ground-based remote sensing measurements aiming at improving the observability, understanding and predictability of hydrological and ecological processes at catchment scale (Li et al., 2009).Observations focused on six different areas with landscapes ranging from desert steppe and Gobi desert to grassland and irrigated farmland.Airborne data used in this study were acquired over the Yingke area.The Yingke Oasis, located to the south of Zhangye city ( 100• 25 E, 38 • 51 N, 1519 m a.s.l.), is a typical irrigated farmland.The primary crops are maize and wheat, with fields often separated by tree hedges.This site was selected for its interest in investigating crop evapotranspiration, bio-geophysical and structure parameters of crops, interaction between groundwater and surface water, and irrigation.

Airborne LIDAR
The WATER field campaign has been completed with an intensive observation period.Twenty-five missions were flown in 2008 with different sensors.This study is based on the use of a LiteMapper 5600 imaging LIDAR, whose major characteristics are a wavelength of 1550 nm, a pulse of 3.5 ns at 100 kHz and a scan angle range of ±22.5 • .The spatial density for a flight height of 800 m above the ground is 4 impacts per square meter.After correction of the raw data to account for the attitude of the plane and for the precise location of the sensor, each return signal can be translated into an accurate 3-D georeferenced point.The resulting point cloud is then processed to extract the minimum and maximum elevation values from the population of points contained in each square meter grid.The lowest elevation point is used to derive the elevation of the ground.After removing local aberrations, the resulting map is a Digital Elevation Model (DEM), expressed in terms of altitude above mean sea level (a.m.s.l), in meters.If the surface is a solid block, or a bare soil, the minimum and maximum elevation values are equivalent.However, for scarce vegetation structures, the difference between the lowest and highest elevations reflects the depth of the vegetation canopy.Therefore, it is possible to separate the land surface topography, represented by the DEM, from the elevation of the top of the vegetation canopy.The later is called the Digital Surface Model (DSM), and is also expressed in terms of altitude a.m.s.l.Both digital models have a spatial resolution of 1 m.The LIDAR flight was operated the 20 June 2008, and the scene covers an area of 7.2 km 2 .A 3-D view of a subset of the entire dataset is presented in Fig. 2.

Meteorological data
The Yingke Oasis experimental site is permanently instrumented with an Automatic Weather Station (AWS).The 19 456 457 station records air temperature, wind speed and direction at 2 and 10 m, and air pressure, relative humidity, precipitation, net radiation, soil heat flux, soil temperature and water content every ten minutes.Moreover, latent heat flux, sensible heat flux and water vapour concentration are obtained from eddy correlation systems with an integration step of 30 min.
Six atmospheric soundings were performed during June and July 2008 with GPS-tracking balloons.The instruments onboard have measured air temperature, relative humidity, air pressure, wind horizontal component and direction, mixing ratio and some information about localization and altitude.They were used to identify neutral atmospheric stability conditions, as well as the height of the top of the planetary boundary layer (PBL).The height of the top of the boundary layer is mainly useful in the CFD model in situations of non-neutral atmospheric stability.In such cases, the vertical boundary may be considered as an impervious slip wall.In our case studies, as we always have neutral stability situations, this parameter is not significant, but must still be provided.

Implementation of the approach from (MacDonald et al., 1998)
The canopy height obtained by difference between the DSM and DEM gives the distribution of roughness elements over the entire area.Considering a subset grid of this area with a surface A T , and the total surface of roughness elements A p within this subset, it is possible to compute a plan area index for the grid.The separation of pixels between roughness elements and intervening space was performed by defining a height threshold from the vertical distribution of pixels.It should be noted that a h ≈ 0 within a LIDAR grid can either mean that there is no vegetation within this grid, or that the canopy is homogeneous and dense enough to prevent impulses to reach the soil underneath.However, in both cases it was considered that such grids belong to intervening spaces in between bigger roughness elements.In the following calculations, the threshold was defined as 12 cm.In the same way, pixels of the same subset grid of surface A T can be projected on a plan surface orthogonal to the airflow, giving a total integrated surface A f used to compute the frontal area index.21 A tool was developed for these purposes.Considering a given grid size, a number of wind direction (2, 4, 8. . . ) and an input pixel size, the tool will sequentially compute the plan area index, the frontal area indexes for each wind directions, and the associated displacement heights and roughness lengths following the two formulation of (Raupach 1994) and(MacDonald et al., 1998).The tool can optionally generate views of frontal surfaces for various wind configurations, with associated roughness elements view within a grid, as illustrated in Fig. 3.In this figure, a view of the obstacles (top) is associated with a view of the frontal profile of all the obstacles inside a grid of 100 by 100 m (bottom).The integral of this profile gives the frontal surface for the wind orientation considered.In the Fig. 3, we present the results for 4 different wind orientations.For such an area, covered with field crops and tree hedges, the variation of the frontal surface with the wind orientation is significant.
Depending of wind direction, the shape of the frontal surface opposed to wind flow can change significantly, as illustrated on Fig. 3, with frontal area index ranging from 0.078 to 0.108.In this particular example, the effect of the orientation of the airflow as compared to the orientation of the tree hedges explains most of the variation, with high values when the wind is perpendicular to the wind flow (45 • and 90 • ), and lower ones when it becomes nearly parallel (0 • and 135 • ).

Configuration of the Computational Fluid Dynamics model
The Canyon CFD requires the input of a Digital Elevation Model and at least one local set of wind profile properties for initialization, i.e. wind speed and direction at two levels, and the height of the top of the Planetary Boundary Layer (PBL).
As mentioned earlier, the later parameter is mostly required in non-neutral atmospheric conditions (Lopes, 2003), which are never considered here.It can use a roughness element height map whenever available, or assumes this height constant on the entire scene.In this study, the Digital Surface Model is used to document the height of the elements on the entire scene.Therefore the model can account for both the topography and the surface stress from roughness elements.
It should be noted that the Yingke area is almost planar, with a very slight slope from West to East leading to an altitude difference of nearly 30 m over the 2400 m swath of the LI-DAR path.The AWS wind speed and direction measurements at 2 and 10 m are used to initialize the profile, together with the PBL height obtained from nearly simultaneous atmospheric soundings.In the following experiments, meteorological contexts are limited to neutrally stratified PBL conditions, leading to the selection of five simulation periods listed in Table 1.In each of these five simulation periods, the top of the PBL was identified at 700 m above ground.

Wind field computation
Wind fields were computed with a ground resolution of 25 m and 15 levels from 5 to 870 m above ground.A control of the computed wind speed at 2 and 10 m with original values from the AWS reveals that speeds can significantly vary.
Table 2 shows differences between measured and simulated wind speed values.The differences are more important near the surface, with a mean underestimation of 40% for CFD wind speeds at 2 meters, but reduce at 10 m, with a mean underestimation of 10% and mean overestimation of 20%.
As quoted in Sect.2, this is due to the fact that the wind speed provided in input is used as an initial guess for an iterative solving of the transport equation in the three dimensions (Lopes, 2003).The use of lower resolution wind fields from a numerical weather prediction model to define a 3-D forcing of the wind speed should give better results, but couldn't be experimented in this case due to a lack of such kind of dataset.
Output wind fields are affected by a border effect on the upstream boundaries of the scene (e.g. on the lower left image of Fig. 4).This imposes to discard results within the first 150 m north and east of the fields.It could however be overcome following a nested scale approach, with use of a lower resolution regional DEM to compute a first initialization field to be used in place of the AWS initialization measurements.This couldn't be performed at this stage of the study.

Roughness length processing
LIDAR data were processed to compute the plan area index of the scene, the frontal area indexes for each wind directions, and associated displacement height following both the approaches from Raupach and MacDonald.The average height of roughness elements within the scene leads to choose a grid size of 100 m, i.e. ten times the height of most of the obstacles to airflow, also tree hedges usually reach 30 m, and up to 38 m for some trees.However, the following calculations stick to the 100 m grid to preserve some granularity.Values of λ p were found to range between 0.08 to 0.64 for tree arrays and some building groups.This leads to d 0(Raupach) / h v values mainly between 0.35 and 0.45, while d 0(MacDonald) / h v values range from 0.2 for bare soil, and up to 0.7 for dense low tree arrays.The lowest λ f values are of 0.025, but can reach 0.2 for grids containing tree hedges.These values are globally rather low because the area doesn't contain regular arrays of high elements, but rather one-line obstacles like the alignments of trees.Values of z 0m are very similar from one orientation to the other, e.g. with a variation of the order of ±8.10 −3 m between values obtained with a 51 • and 270 • airflow.However, significant variations with wind direction in frontal area index, and as a consequence in roughness length, can be obtained for grids containing tree line structures, as illustrated on Fig. 3.The difference of values of z 0m between the formulations of Raupach and MacDonald are more significant, and related to the larger values of displacement height obtained over densely vegetated surface from the MacDonald's formulation.Indeed, z 0m(Raupach) ranges from 0.015 to nearly 0.51, with a maximum z 0m(Raupach) / h v of 0.142, while z 0m(MacDonald) from 0.015 to 0.195, maximum z 0m(MacDonald) / h v of 0.120.CFD based roughness obtained from the inversion of equation (1) using wind fields give a rather different view of the surface drag effect.Although computations are made at a 25 m ground resolution, a grid value expresses the effect of surface stress upstream on the entire footprint of the profile used in the inversion, while the geometrical approach can only account for the frontal density of obstacles within the calculation grid.Here it is assumed that the footprint for the selected neutral conditions is ten times the height of the profile.To obtain results at local scale, the wind field levels from the ground up to 30 m are used to compute the roughness length, for an assumed footprint size of 300 m. Results presented in Fig. 4 illustrate very well in particular the shelter effect of tree hedges, and simulations differ significantly from one wind direction to the other.Here z 0m values are of the order of 0.02-0.03for low vegetation areas, 0.12-0.2for  corn fields, but can reach values as high as 0.8 and even 1.1 nearby tree hedges areas, depending of the orientation of the airflow as compared to the orientation of the hedges.

Discussion
A strict grid per grid comparison between geometrical and CFD based results may not be relevant.Indeed geometrical approaches account for airflow orientation, but they cannot reproduce the footprint effect of upstream roughness elements.These approaches are designed for regular arrays of roughness elements.It can very well account for heterogeneity in terms of grassland with staggered arrays of trees, or any configuration where local heterogeneity tends to a mesoscale homogeneity.However, in such a complex land cover context, the CFD approach proves to give a much finer view of interactions between the airflow of the structure and orientations of roughness elements of significant height.That said, geometrical and CFD based z 0m tend to converge on large, open areas covered either by bare soil, grassland, low field crops, and even to some extent on some corn fields.For an example, Fig. 5b shows in green output grids where z 0m(Raupach) match z 0m(CFD) within an interval of ±0.05 m.Compared to the Digital Surface Model presented on Figure 5a, and considering that the wind direction is 51 • , it appears rather clearly that beside areas affected by a significant shelter effect, both approach tend to give comparable re- sults.This suggests that geometrical formulation could give more comparable results on natural heterogeneous land covers present in the region, like the sparse grassland and low trees land covers.
It should be mentioned that the values obtained from the CFD wind fields for grids containing tall trees might not be correct.In these calculations, it was decided to use the levels of wind fields within the first 30 m to stick to the 300 m footprint, while in some areas some trees can reach up to 38 m.Further investigations are needed to check the quality of the wind speed estimates in the lower part of the boundary, but it seems clear that in cases where roughness elements can reach such a height, the footprint size should be reconsidered, e.g. by using the first 60 m of wind profiles.
It should also be noted that results from the two approaches could only be compared because of the very low variation of the topography over the scene.The use of the difference between the Digital Surface Model and the Digital Elevation Model in the computation of the frontal area index cannot account for a more significant variation of the elevation, while the combined use of the DEM and DSM in the CFD could still give consistent results.
Finally, it must be emphasised that both the geometrical and CFD approaches assume roughness elements to be solid blocks.In both cases, the Digital Surface Model is used to derive an obstacle height that is used either to estimate a frontal surface and intervening spaces, or as a first assumption on local roughness.But none of these approaches will account for the porosity of vegetation structures.In particular, it is obvious that such a representation of obstacles like the tree hedges will lead to overestimate their effect on the flow, while these trees will mainly oppose a resistance to the airflow on levels with the highest foliage density.Therefore, the computation of the roughness length over complex land cover would require accounting for the vertical structure of the canopy.This could be achieved by use of the full waveform of the LIDAR measurements, instead of statistics on points used here.Although the full waveform was retrieved during LIDAR acquisition over the Yingke oasis, these dataset still require further pre-processing, and couldn't be exploited in this study.

Conclusion and perspectives
Hydrological and micro-meteorological studies based on the modelling of surface heat exchanges from radiometric observations can benefit from the contribution of very highresolution LIDAR digital elevation and digital surface models over a complex land cover.The geometrical characterisation of the surface topography but also the structure of the roughness elements paves the way for a more accurate modelling of aerodynamic processes, and in particular a detailed estimate of the surface roughness.The implementation of the geometrical approaches to compute the plan area index and the frontal area index, together with the formulations from Raupach and MacDonald in a single tool is of very general purpose and could be used either on vegetated or urban areas, provided that the local heterogeneity tends to some homogeneity at the fetch scale.On the other hand, the combined used of the DEM and the DSM in a CFD model proves to account for the complexity of the land cover, in particular for staggered structures of tall roughness elements.However, the spatial meaning of the values is different from the gridded geometrical approaches, as a 25 m resolution grid actually accounts for the upstream surface stress within its own footprint.It is also emphasized that both the geometrical and CFD based approach rely on a simple representation of the roughness elements, and do not account for the porosity of foliage structures to the airflow.The definition of the exact footprint of such computations still needs to be investigated.And a cross comparison of results from the CFD based approach with ground measurements at footprint scale could provide a first validation of the results.Moreover, a general analysis of the structure of the landscape along the airflow should allow for an adequate definition of the footprint size and related wind fields levels to be used in the inversion.Results would also benefit from a nested scale computation of the wind fields.The use of coarser DEM over a larger area for the initialization of the high-resolution computations should remove any border effects.Finally, the use of such approaches over other land cover types, but also more accentuated topographies within the Heihe basin could give an extended view of the adequacy of both approaches in various contexts.

Figure 1 :Fig. 2 .
Figure 1: (left) location of the HeiHe River Basin in the People's Republic of China and 458 location of the three experimental areas within the basin (Source: Li et al 2009); (right) 459 detailed view of the Yingke Oasis.The extent of this scene is the same as the coverage of the 460 LIDAR data, ie.7.2 km 2 .This image and the LIDAR data were acquired simultaneously, 461 using a CCD camera onboard the aircraft.462 463

470 471Figure 3 :
Figure 3: A 100x100 meters subset covered with field crops and tree hedges is presented to 472 illustrate the results of a geometrical processing.The four images at the top present the same 473 scene considered from four different view angles, corresponding to wind directions of 0°, 45°, 474 90° and 135° (from North).The colours express the vegetation height from 60cm (blue) to 25 475 meters (red), and the tick marks are pixel coordinates.The four bottom images illustrate the 476 frontal area of this subset from each of the four orientations of the wind.Abscissa marks are 477 pixel coordinates, while the ordinate marks represent the obstacle height in meters.In this 478 figure, Lp refers to the plan area index, and Lf to the frontal area index.479

Figure 4 :
Figure 4: Roughness length maps derived from the LIDAR data over the Yingke area (7.2 482 km 2 ) for wind flows from N-E (51°), W-NW (295°) and W (270°), presented from left to 483 right, and related results following the approaches from MacDonald, Raupach, and from the 484 CFD, from top to bottom.Arrows represent wind directions accounted in both geometrical 485 and CFD based calculations.486

Fig. 4 .
Fig. 4. Roughness length maps derived from the LIDAR data over the Yingke area (7.2 km 2 ) for wind flows from N-E (51 • ), W-NW (295 • ) and W (270 • ), presented from left to right, and related results following the approaches from MacDonald, Raupach, and from the CFD, from top to bottom.Arrows represent wind directions accounted in both geometrical and CFD based calculations.

23 489Figure 5 :
Figure 5: (a) Roughness element height from the DSM (in meters); (b) areas where both 490 z0m(Raupach) and z0m(CFD) match at ±0.05 meters for the calculation with a N-E wind are 491 represented in green.Both figures cover the 7.2 km 2 area of interest.492 493 494 Fig. 5. (a) Roughness element height from the DSM (in meters); (b) areas where both z 0m(Raupach) and z 0m(CFD) match at ±0.05 m for the calculation with a N-E wind are represented in green.Both figures cover the 7.2 km 2 area of interest.

Table 1 .
Wind speed and direction measured at Yingke AWS for the selected neutrally stratified Planetary Boundary Layer conditions.The height of the top of the PBL was identified at 700 m in each case.

Table 2 .
Comparison between wind speed measured at the AWS and simulated wind speed values obtained with the CFD model.