A simplified method on estimation of forest roughness by use of aerial LIDAR data

In the last decade, satellite‐based measurements combined with local land cover information have produced datasets with a very detailed land cover description. CORINE Land Cover (CLC) dataset is owned and maintained by the European Environmental Agency (EEA) and published at the agency's website. Another remote sensing tool, developed in the same period, is the terrain LIDAR scanners with very high resolution and porosity information. In the current study, LIDAR scans of mainland Denmark with 0.4 m resolution were used to estimate the aerodynamic roughness of large forests, the borders of which were defined with the help of the CLC dataset. The results are compared with available in situ measurement results from the scientific literature. There was a generally good agreement between calculated and measured displacement height values but less so for aerodynamic roughness values due to the employed spatial averaging process. The results reveal a promising application that can be used for forest parameterization within modeling tools.


| INTRODUCTION
Until recent years, wind park developers used to prefer wind farm sites away from forested areas due to high turbulence in and around forests, which causes reduction in power production and shortens the life span of wind turbines. However, as prime areas of low relief are gradually being depleted, there is a surge of interest in the use of forested areas for purposes of wind farm siting, 1 and with that, a need to better quantify the effect of forests on the wind flow in relation to wind turbines. Some studies suggest that if a turbine is located deep within a forest, the logarithmic wind profile gets affected only/primarily through displacement height 2 while additional turbulence is limited; on the other hand, if the turbine is located at the edge of a forest high turbulence is present. In reality, even if wind turbines are located some distance away from the forest the internal boundary layer created by the sudden roughness change from a forested area to nonvegetated area creates a tall internal boundary layer and affects the area of interest. 3 Therefore, at any point in or around a forest, one should be able to correctly identify the forest parameters at the flow modeling stage in order to achieve the best possible prediction. Two of the key forest parameters in wind flow modeling are the aerodynamic roughness and the displacement height of the forest, which are used to predict the wind speed at hub height through logarithmic wind profile. 4 Other properties of a forest, such as tree height, tree type, or crown height, can be identified by visual inspection, but the roughness length of a mixed vegetation area is not easy to determine without measurement. Several studies have tried to measure the aerodynamic roughness of forests with conventional | 3275 BINGÖL masts located deep inside forests or with a combination of new generation wind profiler LIDARs and masts at forest edges. 3,5 Although these campaigns succeeded in offering better informed assumptions for aerodynamic roughness, they did also identify some variables that have a large impact on the final calculation, such as the type and heterogeneity of the trees, 6 gaps within the forest, and the distance of trees from each other. 7 Therefore, arriving at a basic rule of thumb has proven to be difficult. A decade ago, studies and modeling tools assumed that the roughness of a forest is between 0.8 and 1.0 m, where the highest value was commonly used in wind engineering practice and the scientific literature. 8 Studies after 2010 suggest a higher range of values, between 1 and 1.5 m. An averaged roughness value of 1.5 m is adopted by the industry standard modeling tool WAsP, developed by the Technical University of Denmark, Department of Wind Energy. 9 Most recent study about the subject employs the CFD approach and allometric relationships to retrieve forest aerodynamic parameters from terrain scanning LIDAR directly. 10 In this study, a new approach for calculating forest roughness by use of aerial LIDAR is presented, where the forest borders are identified through a satellite land cover dataset. The suggested methodology is based on a spatial averaging method derived from the currently known forest roughness models explained in further sections. In the next chapter, a description of the point cloud dataset of Denmark measured by an airborne LIDAR and the land use classification system CORINE Land Cover (CLC) are introduced. This is followed by the underlying theory and the numerical method for calculating the aerodynamic roughness values and for spatial averaging within each forest zone. Finally, the results are compared with available in situ measured forest roughness values from the scientific literature.

| DATASETS
The Airborne LIDAR scanning technique was introduced in since 1996, 11 primarily/initially for military use. The major bottlenecks of the technology that prevented wide civilian use were the resolution of collected data points and the processing speed of the point cloud data. If the data were collected with a low resolution to begin with or if it was analyzed based on a limited number of points from the entire data cloud, the output was not satisfactory. With the generation of fast computers at the turn of the century, it became possible to analyze large point cloud datasets, and therefore, the airborne LIDAR scan datasets started to become a desired engineering tool. Currently, the technology is used to create high-resolution records of topographical information. 12 Figure  1). a One can easily query each dataset and obtain the terrain elevation above mean sea level and the vegetation characteristics above the terrain with high accuracy.
Although the airborne data give detailed information on the tree characteristics, it does not directly provide forest borders, which are necessary for the spatial averaging of roughness values. In order to generate this information, a separate step is required, utilizing information from another large-scale, open-source study. In the last decade, satellitebased measurements in combination with other global or local land cover information have produced datasets with a very detailed land cover description. The CLC dataset is owned and maintained by the European Environmental Agency (EEA) and is published at the agency's website. b The first version of the database was finalized in the early 1990s as part of "the European Commission programme to COoRdinate INformation on the Environment" (CORINE). The database is also referred to as CORINE Land Cover (CLC) dataset. The polygonal borders of each forest in the CLC dataset were superimposed on the LIDAR generated point cloud in order to perform spatial averaging calculations. Methodological details regarding image processing, how polygon borders were defined, and the description for each land cover class are available in a number of publications 13,14 and the official technical reports of CLC written by the developers. 15,16

| METHOD AND THEORY
In the lower part of the atmosphere, above the ground, where the wind is affected by the earth's surface, the vertical wind profile is defined as 17 : where z is the measurement height above ground level, u* is the friction velocity, mostly a product of the covariance of fluctuations of the u and w components of the wind, κ is the von Karman constant, which varies between 0.35 and 0.43 but it is accepted as 0.4 for boundary layer meteorology, 18 d is the displacement height if there is a tall obstacle similar to forest vegetation present, and z o is the aerodynamic roughness.
It is important to understand aerodynamic roughness, z o , and the displacement height, d. Aerodynamic roughness is the height above the ground or above the displacement height where the wind speed is zero. This is the height where the flow is stopped due to friction of the surface that it is close to. The displacement height is the effect of an obstacle, such as trees or buildings, that does not let the wind go through it directly but rather shift it up. Over a vegetation, this reaction will shift the wind profile upwards to a certain level but not at the exact height as the tree. There will be a region at the top of the canopy where the wind can pass through, and therefore, it will slow down to zero speed somewhere below the tree height. That shift of the profile is called the displacement height. A scaled version of the vertical profiles with and without displacement is presented in Figure 2. The vertical profiles are derived from Equation 1 with realistic values up to two tree heights. The thick continuous vertical profile on the left corresponds to a flat terrain and neutral conditions, where the profile goes down to z o height. The double line profile on the right, above the trees, is the shifted profile upwards at a height of d. Such a profile can be found not at the edge of the forest but rather inside the forest. There is a transition zone before and after the forest edge that should be formulated differently but this is beyond the scope of this study.
Several studies 3,4 have indicated an inverse relationship between d and z o . One of the most common parameterizations for this relationship is Raupach's formula, 19 which is based on even earlier theory. 20 Massman et al later suggested a more detailed approach. 21 As a result, a combination of these methods, briefly described below, is used in recent studies. 4 (1) First, three terms are needed for the forest parameterization: tree height (h), tree crown width (b), and distance between trees (D) which leads to well-known frontal area density functions derived by Raupach et.al 22 : The λ is used to derive the displacement height of a tree if the drag coefficient c dI of the tree is known with the Equation 3. For c dI, a constant value of 7.5 is suggested. 19 Following Raupach 22 : where γ is the ratio between the wind speed at the height h, U(h) and the friction velocity, u * . ψ is a function describing the deviation from the inertial sublayer. 19

| Simplifications
In order to achieve a uniform method for all Denmark with the same quality LIDAR data, a few simplifications had to be made, mostly regarding tree allometry. There are three forest classes in CLC 15 : (a) Class 311: Broad-leaved, (b) Class 312: Coniferous forest, and (c) Class 313: Mixed of Classes 311 and 312. Aligned with this grouping, three assumptions were made. The aim of the first three assumptions is to calculate the frontal density (λ) by means of the forest density (ξ). The last assumption refers to the distinction between short vegetation forms and trees.
1. Any nonhomogeneous forest can be defined as homogeneous by assuming equidistant spacing (D) of trees with averaged parameters. Every forest parameter per single tree can be averaged, that is, tree height, tree width, and distance between trees, so that a forest such as the one schematically depicted in Figure 3-Left can be averaged to a forest represented by Figure 3-Right. 2. Every single tree in a nonhomogeneous forest can be treated separately, and tree parameters can be computed. Later, the borders of the forest can be used to average these values spatially and generate a single set of parameters describing the whole forest. 3. Tree allometry studies give a relationship between the tree height h and crown width b (See Appendix). 4. Any vegetation below 1.3 m height is ignored in the calculation because vegetation starts to display tree characteristics above such height. 23 When the above four assumptions are made, the relationship between forest displacement height (d) and ground The frontal density (λ) can be derived by knowing only the tree height (h) and the forest density (ξ) while the displacement height (d) becomes a function of h and calculated λ. Finally, the aerodynamic roughness can be computed as the 1/3 of the difference between h and d.

| RESULTS
Denmark has around 3700 m 2 forested land in the CLC classifications 311, 312, and 313 which are 20%, 43%, and 37% of the total, respectively. To apply the assumptions for all the available 2927 polygons, a processing algorithm was written, including the following steps: • Surface data of each forest is calculated from the point cloud data of Denmark, using the borders of the CLC polygons. • Empty areas within the borders of each forest are assumed to represent vegetation forms below 1.3 m high. The empty area is compared to the total area and ξ is derived. • The rest of the point cloud is treated as tree canopy height.
For each tree height (h), a specific crown width (b) is calculated based on the forest class. This step produces the relationship between λ and ξ.  Figure 4). The algorithm and the methodology detect the sparse and low vegetation successfully. In Figure  5, the forests close to the village have several gaps within them and the roughness and displacement height values are detected believed to be reasonable by the author. The validation cases for known data from literature are given in next subsection. Even though the calculated results are within the limits of known and expected values, it is necessary to validate the results with in situ measurements. Two sites from Denmark, Corselitze and Sorø, were chosen to compare the results. The sites have available allometry data collected as h, b, and D and derived the λ from several previously published studies, 4,24,25 where a full description of the sites can also be found. Surface

| 3279
BINGÖL data from the aerial LIDAR data were obtained for these two forests following the steps described above. The histograms of the data from both forests are given in Figure 6. The peaks for the lowest tree heights correspond to the nonvegetated or low vegetation (≤1.3 m) areas within or at the borders of the forests. The Sorø forest is characterized by rather variable tree heights which is clearly seen both in the histogram and in the aerial photographs, as the Southern part has markedly lower vegetation. On the contrary, Corselitze appears to be a much more uniform forest, with only two small nonvegetated patches. When the methodology described in this study is applied for both of the forests (Figure 7), averaged d and z o are calculated (Table 1). In order to make the results comparable with the single point measurement values from the literature, maximum values of the tree height (h max ) and the derived displacement height (d max ) are also presented. Single point measurements make use of the highest group of trees around the measurement location; therefore, measured tree height is also marked as h max . The density of the two sites (ξ) is similar but the frontal density (λ) of Corselitze is a little higher. Calculated displacement heights (d) are almost half the values of the results from the point measurements but this is due to the spatial averaging for the whole forest area. If maximum tree height (h max ) values are used to calculate the displacement height (d max ), similar displacement heights with the in situ measurements (d') are observed.

| CONCLUSION
By means of areal LIDAR scan and satellite land cover datasets, a simplified methodology was developed to calculate the roughness (z o ) and displacement height (d) of the forests of mainland Denmark. The proposed algorithm uses point cloud datasets and borders from CLC land cover datasets to calculate spatially averaged results for the whole forest. The method gives an opportunity to wind engineers to make realistic assumptions on these important forest parameterization values. The new method can be used with modeling tools which can employ both the displacement height (d) and aerodynamic roughness (z o ). The method was applied to the entire Danish mainland and was compared with point measurement results from two forests available in the scientific literature. The first site, Sorø, displays nonhomogeneous forest characteristics and includes two different tree types and many cleared sections within the forest borders. The second site, Corselitze, on the other hand is mostly a homogeneous forest with only a few cleared areas. While the new method is averaged throughout the forest borders, results from research studies were obtained from in situ vertical wind profile measurements at point locations. Therefore, the calculated displacement height (d) differs from the point data results, while the maximum displacement height (d max ) is significantly close the measured value.
However, to the extent that the calculated value of displacement height (d) is used to compute the roughness length (z o ), it is believed to be a representative of the overall displacement height of a forest site. This is further supported by the generally good agreement between calculated and measured roughness values z' o being only 15%-25% lower than z o , presumably due to the spatial averaging process which takes nonvegetated areas into account. Importantly, the calculated values are closer to the average roughness value of 1.5 m for forests, which is the latest recommendation by the WAsP   26 Again, due to the spatial averaging one can use the methodology to parameterize any forest around for a wind turbine location away from forested area. Also, any wind turbine erected within the forest, again away from the edges, can use the methodology to a certain level by keeping in mind the displacement height is averaged. The use of the methodology at the forest edges (in or out) should give a large error due to the spatial averaging.

ACKNOWLEDGMENTS
The study is supported by the Technological Research Council of Turkey (TÜBİTAK)-ERANET project 216M505 and the H2020-Eurostars programme, Grant Number E!10037. The author also would like to thank the Danish Geodata Agency and the European Environmental Agency for their long-time efforts on publishing and maintaining the remote sensing datasets that are used in the study.

RELATIONSHIP BETWEEN CROWN WIDTH AND TREE HEIGHT
An allometry study is used to derive the averaged relationship between the crown width (b) and the tree height (h). 27 The data taken from the study separated for CLC311 and CLC312 forest types are shown in Table A1. Both datasets are combined for CLC313 forest type. The available data are fit into a second degree polynomial equations with two coefficients, b(h) = c 1 ln(h) + c 2 which is a simplified single-variable type commonly used in allometry studies. 28 The results of the linear fit are given in Figure A1. These fit functions are used to calculate the crown width based on the tree height.