ESTIMATION OF AERODYNAMIC ROUGHNESS AND ZERO PLANE DISPLACEMENT USING MEDIUM DENSITY OF AIRBORNE LIDAR DATA

This paper presents a framework to estimate aerodynamic roughness over specific height (zo/H) and zero plane displacement (d/H) over various landscapes in Kelantan State using airborne LiDAR data. The study begins with the filtering of airborne LiDAR, which produced ground and non-ground points. The ground points were used to generate digital terrain model (DTM) while the nonground points were used for digital surface model (DSM) generation. Canopy height model (CHM) was generated by subtracting DTM from DSM. Individual trees in the study area were delineated by applying the Inverse Watershed segmentation method on the CHM. Forest structural parameters including tree height, height to crown base (HCB) and diameter at breast height (DBH) were estimated using existing allometric equations. The airborne LiDAR data was divided into smaller areas, which correspond to the size of the zo/H and d/H maps i.e. 50 m and 100 m. For each area individual tree were reconstructed based on the tree properties, which accounts overlapping between crowns and trunks. The individual tree models were used to estimate individual tree frontal area and the total frontal area over a specific ground surface. Finally, three roughness models were used to estimate zo/H and d/H for different wind directions, which were assumed from North/South and East/West directions. The results were shows good agreements with previous studies that based on the wind tunnel experiments. * Corresponding author


INTRODUCTION
The information related to energy and mass transfer between land surface and the atmosphere is really important for both numerical weather prediction and climate studies (Su, et al., 2001).Usually, Monin-Obukhov similarity (MOS) theory was used for those applications with two important parameters (i.e.aerodynamic roughness length (z o /H) and zero-plane displacement (d/H)) (Tian et al., 2011).MOS theory serves an effective tool for turbulence statistics prediction in the horizontally-uniform surface layer for applying MOS formulation to relate surface fluxes of momentum, heat and water vapour to meteorological variables (Wood et al., 2010).z o /H can be defined as the height at which the wind speed become zero (Lefsky et al., 2002) while d/H can be explained such that: "When some roughness elements are closely arranged, there will be a higher surface formed by them, equivalent to which is the distance by which the height moves up, and this distance is called zero-plane displacement (d/H)" (Hu et al., 2015) Estimation of z o /H and d/H values can be done through two general methods i.e. experimental based method and remote sensing based method (Tian et al., 2011).Experimental-based method employed vertical wind profiles and micrometeorological theory provides locally estimated values (Schaudt and Dickinson, 2000).Nevertheless, this type of measurements consumes high cost of implementation, which limits its implementation over greater spatial scale (Paul-Limoges et al., 2013).This method provided z o /H values for the local surface features, which lead to the limited z o /H data in energy and water-heat exchange (Hu et al., 2015).The limitations of the experimental-based method can be overcome by remote sensing method.Advances in remote sensing technologies have introduced new concept in deriving three dimensional structure of land surface (Paul-Limoges et al., 2013).Furthermore, remote sensing method also can be utilized to portray the vegetation structures with their spatial variation (Hu et al., 2015).
Previous studies have shown that Airborne Light Detection and Ranging (LiDAR) has become one of the remote sensing technologies that able to estimate the height surface objects directly for z o /H estimation (Brown and Hugenholtz, 2012).The acquiring process of vegetation three-dimensional (3D) structure become more effective by means of LiDAR technology (Hu et al., 2015).It has been proved that this technology is a key method in retrieving various forest structural parameters (Kljun et al., 2015;Tian et al., 2011) as well as detailed estimation of vertical structural information of land surface (Rahman et al., 2015).For z o /H and d/H estimation, the ability of LiDAR technology to produce Digital Surface Model (DSM) of the area plays a significant role in calculating the plan surface density and frontal area of the trees (Raupach, 1994;Macdonald et al., 1998).
In general, the estimation of z o /H and d/H required detailed estimation of individual trees structure information (Rahman et al., 2015).The ability of airborne LiDAR technology in getting such valuable information of tree structure enhances the estimation process.However, some gaps have been identified where the performance of LiDAR system in tropical region still require thorough investigation that give an impact to the quality of the extracted parameters.This is due to the fact that higher vegetation density in tropical region reduces the laser penetration couple with error in LiDAR system configuration that significantly affects the basic processing of Digital Terrain Model (DTM) and Canopy Height Model (CHM) generations.The objective of this study is to explore the use of airborne LiDAR data for the estimation of the z o /H and d/H over various landscapes in Kelantan State.

MATERIALS AND METHOD
In general, the methodology of this study is divided into five main stages namely 1) data acquisition, 2) data pre-processing, 3) estimation of tree attributes, 4) estimation of aerodynamic roughness length (z o /H) and zero plane displacement (d/H), and 5) validation of results as shown in Figure 1

Description of Study Area
This study investigates the estimation of z o /H and d/H value using airborne LiDAR data in certain parts of Kelantan.Kelantan is located in the northern-eastern of Peninsular Malaysia with the total area of 5,830 square miles.This study focused on a specific region of Kelantan with the total area about 1160 square miles, which are covered with forest area, plantation area, water bodies, and residential area.Figure 2 shows the boundary and locations of study area.

Pre-processing of Airborne LiDAR Product
In the data pre-processing stage, the airborne LiDAR product has went through several steps of data pre-processing, which aims at generating CHM.CHM gives height information for tree canopies and this information is very crucial to estimate other related properties with support from the existing allometric equations.The CHM has been produced at 3m spatial resolution.

Estimation of Individual Tree Attributes
Accurate forest attributes such as number of trees, height, and diameter at breast height (DBH) should be acquired effectively for efficient and economical forest management (Chang et al., 2008).The estimation of tree attributes begins with individual tree crown delineation.The accuracy of tree crown delineation depends on the quality of CHM (Jakubowski et al., 2013) and method of crown segmentation such as inverse watershed (IW) segmentation method (Rahman and Gorte, 2009).This method was applied by inversing the CHM model and applying the morphological watershed segmentation, which assumes that the highest point of a tree to the lowest point of a tress crown that marks the edge of crown.Figure 3 shows the CHM surface, location of trees and the boundaries of each tree in the selected forest area.
The tree crown segments are then used in deriving individual tree attributes i.e. tree height ( ), crown width (CW), crown depth (CD) and diameter at breast height (DBH).Specific allometric equations were used in estimating those tree attributes.and tree locations (yellow points) The allometric equation (Equation 1) that relates DBH and maximum tree height was used to calculate individual tree DBH (Okuda et al., 2004).Individual tree crown width is obtained by simplifying tree crown diameter from the tree segments (Equation 2) (Morton et al., 2014).Equation 3is used to estimate individual tree crown depth from the tree height value (Morton et al., 2014). (1) (2) (3) where is the maximum height of tree with a forest stand.
All the identified tree attributes are further used in estimating the frontal area index, which is explained in the next section.

Frontal Area Index Calculation
This section discussed in details about the estimations of z o /H and d/H.The process started with calculating the tree frontal area index by using semi-automatically method via Model Builder with arcpy module in ArcGIS software.Figure 4 presents the general flowchart model for the model builder.The frontal area is calculated based on the assumption that tree crowns and tree trunk are characterized by round and vertical rectangle shapes respectively.The tree trunks are considered vertical perpendicular to the ground surface (Figure 7).
By referring Figure 4, four (4) main process involved in calculating frontal area index as follows: 1. Grid generation and tree location intersection.Figure 5. Modeler flowchart for grid generation and tree intersection

Tree trunk and crown generation based on different wind directions
Next, the process continues with tree trunk and crown generation based on different wind directions.The process was generalized as shown in Figure 6.All the grids were applied with this routine in getting the tree structure based on the identified tree attributes parameters.The 'x-coordinates' values were chosen as longitude value in order to signify the trees viewed from 0 0 to the North (North-South direction).Meanwhile 'y-coordinates' values were assigned to longitude value as if the trees were viewed from 90 0 to the North (East-West direction).Similar setting used for both crown and trunk generation.The tree trunk is reconstructed by generating a line that joins the crown base location and point that touch the ground surface.The length of the trunk is estimated by subtracting crown depth from the tree height.To illustrate the width of the trunk, the lines were buffered with the value of diameter at breast height (DBH) for each individual tree.The reconstruction of tree crowns was based on the assumption of round shape.In this case, the crowns were generated by buffering the crown centre points with reference to the crown radius values.The processes were repeated for all trees in each grid.Figure 7 summarizes the overall concept of an individual tree reconstruction.
Figure 7. Illustration of tree trunk and crown reconstruction

Dissolving merging trunk and crown for frontal area index calculation
Figure 8 depicts the processing steps for dissolving and merging tree trunks and crowns.The merged structure of individual tree crown and trunk is required in determining the frontal area index over specific grid.Thus, these individual trees then need to be dissolved in order to obtain only frontal area of the trees facet perpendicular to the flow direction (see Figure 9).Finally, the frontal area can be calculated in metre square unit of a grid.Figure 10.Flowchart for frontal area index gridding

Estimation of Aerodynamic Roughness Length (zo/H) and Zero Plane Displacement (d/H)
Aerodynamic roughness values were estimated using different models developed by previous researchers i.e.Lettau (Equation 4), Counihan (Equation 5), and MacDonald (Equation 6).The zero plane displacement is calculated using Equation 7. (4) (5) where is the frontal area index, is plane area of obstacles, is coefficient to best fit the relation with experiments ( for square-array data), is the drag coefficient and  is the Von Karman constant approximately set as 0.4.The frontal area index for a tree can be calculated using equation 8.
(8) where is the frontal area and is the total area covered by roughness element which determines the spatial resolution of the map.The value is calculated for two different possible wind directions (0 0 and 90 0 from North direction).The aerodynamic roughness length value and zero plane displacement were estimated with 500 m and 1000 m spatial resolutions.

RESULTS AND DISCUSSION
As mention earlier, the procedure for estimation of z o /H and d/H were done semi-automatically using available tools developed in ArcGIS software.Table 1 shows the common aerodynamic roughness values assigned conventionally for different landscape settings (Wieringa, 1992).The values are comparable with the estimated values from airborne LiDAR data.

Surface
. The first stage focuses on the data acquisition where airborne LiDAR product (i.e.DTM and DSM) acquired in raster data format.The data is then used in the pre-processing stage where generation of canopy height model (CHM) is performed with 3 m spatial resolution which similar to DTM and DSM.The third stage deals with estimation of tree attributes based on the delineated trees and allometric equations.The fourth stage is devoted for the z o /H and d/H estimations.In this study the estimation of z o /H and d/H were based on the derived individual tree geometrical and existing models i.e.Lettau, Counihan, and MacDonald models.Finally, the results are compared with the values of z o /H and d/H obtained from previous studies.

Figure 1 .
Figure 1.Flowchart of the Methodology

Figure 2 .
Figure 2. Map of study area, Kelantan

Figure 3 .
Figure 3. Sample of delineated tree crown segments (red lines)and tree locations (yellow points) 2. Tree trunk and crown generation based on different wind direction.3. Dissolving merging trunk and crown for frontal area index calculation.4. Gridding of frontal area index.

Figure 4 .
Figure 4. Flowchart of ArcGIS model builder for frontal area index calculation 2.4.1 Grid Generation and Tree Location Intersection Generally, specific grid sizes of ground surface are specially designed in the estimation of z o /H and d/H. Figure 5 shows the process which start by generating the grid with two sizes i.e. 500 m x 500 m and 1000 m x 1000 m.The tree locations (points) are intersected with the grid, which allows separate calculation of z o /H and d/H values for each location in the study area.

Figure 6 .
Figure 6.Flowchart of tree trunk and crown generation based on different wind directions

Figure 8 .
Figure 8. Flowchart of dissolving and merging of trunk and crown for frontal area calculation very rough surface (numerous obstacles) as guided based in Figure 11.The results show some effect of the grid size on the estimated z o /H Lettau's model.The greater grid size gives much lower z o /H values compare with smaller grid size especially for 90 o from the North direction.0 o from the North (1000 m) 90 o from the North (1000 m) (a) 0 o from the North (500 m) 0 o from the North (500 m) (b) Figure 11.Spatially distributed z o /H value with assumed wind directions 0º and 90º from north calculated using Lettau's model (the total area, are (a) 1000 m x 1000 m and (b) 500 m x 500 m) In general, the z o /H values obtained from Counihan's model (Figure 12) much higher than values estimated using Lettau's model.The maximum z o /H value for 1000 m x 1000 m grid falls under the very rough surface with numerous obstacles category.Similar estimation results obtained for 500 m x 500 m grid size where some of the areas fall under the category of closed surface (regular large obstacles coverage) and chaotic surface (500 m x 500 m grid size for 90 o from the North).The results marked some effect of the grid size on the estimated z o /H Counihan's model.

Figure 13
Figure 13 shows z o /H value acquired from the MacDonald's model where the z o /H value lower than Lettau and Counihan models except for 500 m x 500 m grid size of 90 o from the North which higher than Lettau's model.In general, MacDonald's model shows most of the area considered as roughly open whereas occasional large obstacles exists.Direction with 90 0 from the North direction (500 m grid size)

Table 1 .
(Wieringa, 1992)hness length with its corresponding earth landscape(Wieringa, 1992)Lettau's z o /H values range between 0 m to 0.647 m (Figure11).Based on Table 1, the estimated z o /H values can be grouped under the open sea which might come from waterbodies area to