From lidar scans to roughness maps for wind resource modelling in forested areas

Applying erroneous roughness lengths can have a large impact on the estimated performance of wind turbines, particularly in forested areas. In this study, a new method called the objective roughness approach (ORA), which converts tree height maps created using airborne lidar scans to roughness maps suitable for wind modelling, is evaluated via cross predictions among different anemometers at a complex forested site with seven tall meteorological masts using the Wind Atlas Analysis and Application Program (WAsP). The cross predictions were made using ORA maps created at four spatial resolutions and from four freely available roughness maps based on land use classifications. The validation showed that the use of ORA maps resulted in a closer agreement with observational data for all investigated resolutions compared to the land use maps. Further, when using the ORA maps, the risk of making large errors (> 25 %) in predicted power density was reduced by 40–50 % compared to satellite-based products with the same resolution. The results could be further improved for highresolution ORA maps by adding the displacement height. The improvements when using the ORA maps were both due to a higher roughness length and due to the higher resolution.


Introduction
In the past decade, there has been a significant increase in installed wind power capacity in forested areas across northern Europe (Enevoldsen and Valentine, 2016). This has been made possible due to the increased hub heights and improved technologies of modern turbines . However, these forested sites typically represent a challenging environment for wind turbines because of high turbulence and wind shear levels (Arnqvist et al., 2015). The forested landscape is rarely homogeneous, containing both forest edges and roughness changes that affect shear and turbulence and can produce sharp gradients in the flow (e.g. Poëtte et al., 2017). This, in turn, makes the accurate siting of wind turbines in forested areas critical. For all of these reasons, it is important to parametrize forest effects correctly in wind models aimed at the application of wind turbine siting. The aim of the current study is twofold: (1) to introduce a new and automated way of using highly detailed forest information with the WAsP method (Troen and Petersen, 1989) and (2) to investigate how prediction errors change when using different maps to represent the forest.
Starting with the second aim, a short background for the research in forest wind meteorology is given. In surface layer theory, the effect of a homogeneous forested area on the wind speed U is parametrized through a suitable value for the roughness length z 0 as well as a displacement height d in the well-known logarithmic wind profile: where u * is the friction velocity, z is the height above ground, and m represents the stability correction to the profile and depends on z/L, where L is the Monin-Obukhov length (Businger et al., 1971). Based on this equation, experimental studies on the wind profile over forested areas typically present estimates of z 0 and d as a pair (e.g. Thom, 1971;Hicks et al., 1975;Molder et al., 1999;Dellwik and Jensen, 2005;Arnqvist et al., 2015). A dependence among measurable forest characteristics (height and density) and z 0 and d was suggested by Jackson (1981) and Raupach (1994). According to Raupach (1994), a dense forest should be represented by a relatively lower z 0 and a higher d than a sparser forest with the same tree heights. For sites without forest density information, Garrat (1992) recommended setting z 0 = 0.1h and d = 2/3h, where h is the tree height. The Wind Atlas Analysis and Application program (WAsP, www.wasp.dk, last access: 23 November 2017) method (Troen and Petersen, 1989) uses a measured wind climate for the prediction of a wind climate at another location. It is also available as a flow model in the Wind-PRO software, which is developed by EMD International A/S (www.emd.dk/windpro, last access: 12 May 2018). WAsP includes z 0 in two different ways: in the absence of nearby obstacles, the difference between the predicted and observed wind climate is calculated through a geostrophic roughness model that uses Eq. (1) in combination with the geostrophic drag law for vertical and horizontal extrapolation. Heterogeneity in roughness is modelled through a roughness change model, which consists of an internal boundary layer model that is driven by lines of roughness changes on a map (Sempreviva et al., 1990). The speed-up effects due to changes in terrain height are taken into account by a flow-over-terrain model (Troen, 1990).
Using data from small beech forests, Dellwik and Jensen (2000) showed that the internal boundary layer growth closely resembled the prediction of the WAsP roughness change model. In a later study, (Dellwik et al., 2006), it was demonstrated that WAsP could correctly predict the wind climate at the forest sites from observations at a nearby agricultural area, when the z 0 value was set to the values found in the surface layer studies mentioned above, which were significantly higher than those that are typically used for forest studies. However, these studies investigated sites with nearsurface wind observations in an open landscape with small distinct forests, which differ greatly from the northern European sites that have recently started to be developed for wind energy. The northern European sites have a landscape that is dominated by forests with smaller areas of lower-roughness clearings and lakes, leading to study areas where both the observation and predicted sites are located in a high-roughness forested area. This study looks to determine if the roughness map has as much impact at such sites, as was demonstrated for the sites in Dellwik et al. (2006).
There is currently no consensus on how to take d in Eq.
(1) into account in WAsP (Enevoldsen, 2017). Many siting engineers use terrain height products from the Space Radar Topography Mission (SRTM) to describe the orography effects, which includes the vegetation height in the terrain height, meaning that for forested areas, the displacement height is already included to some extent (Kellndorfer et al., 2004;Crockford and Hui, 2007). In Dellwik et al. (2006), a terrain height map that did not include the forest height was used and the displacement height was included by adjusting the height of the wind observations in a post-processing step. However, this post-processing approach does not work close to forest edges since they cause the flow to speed up similar to a terrain height escarpment (Dellwik et al., 2014). In forested sites, measurement masts are often placed in clearings near forest edges, suggesting that the displacement height is particularly important. This study includes an investigation into the impact of adding d to terrain data on model predictions.
The performance of the flow modelling depends greatly on the availability of accurate maps of terrain height and z 0 . It is common for siting engineers to use z 0 maps derived from land use classifications, and many siting tools, such as Wind-PRO, provide these maps for download directly in the software. These maps are quite useful in that they can be obtained anywhere in the world and do not require additional time to locate and implement in the siting tool. However, because the land use categories were not developed for wind energy applications, a specific land use class may not correspond to a unique effect on the wind field. Additionally, neither the forest height nor forest density are typically part of land use classifications, making it difficult to apply the relationship suggested by Raupach (1994), for example, to forested areas. This information may be found either by a site visit or from the forest owners, but it is difficult to find consistent high-quality information over the large areas required for wind resource assessment. This difficulty was one of the starting points for this study and is connected with aim (1) stated above.
An attractive solution to obtaining high-quality information about forest characteristics can be found in the raw data from remote-sensing airborne laser scans (ALSs) for land surface mapping (Clark et al., 2004). During the last decade, ALS technology has seen a rapid expansion and ALS mapping campaigns have been performed for entire countries (e.g. Nilsson et al., 2017). Simultaneously, the price of ALS studies has decreased, making mapping campaigns for in-depth studies possible for both research and commercial wind projects (e.g. Boudreault et al., 2017). The raw data from the scans are stored in standard format, allowing siting engineers to use the same data source for land surface description for different sites (American Society for Photogrammetry & Remote Sensing, 2010). Boudreault et al. (2015) introduced a method that used the raw data from ALS campaigns in canopy-resolving flow models. Here, we translate the derived forest height from ALS data to a roughness map for use in a WAsP analysis. The ALS-based model results are compared with model results based on standard land use maps. The four land-use-based maps tested in this study had different resolutions, and included a varying level of detail. To investigate the importance of resolution, the roughness maps made from the ALS scans were created at resolutions matching the land use models. In addition to evaluating This study performs validation and analysis of the model results for the cross prediction of seven masts located in a predominantly forested area in central Sweden. At such high latitudes, icing on cup anemometers and wind vanes is a common problem for wind measurement (Bredesen et al., 2017). There were periods with frozen cup anemometers and wind vanes: the process for removal of these periods is described in Sect. 2. The processing of the ALS data and their conversion to roughness length and displacement height is presented in Sect. 3. Section 3 also contains descriptions of the land use classification based datasets. The WAsP flow model and the cross-prediction method used for the model assessment are then described in Sect. 4. Section 5 presents the results of the comparison among model performance among the different roughness products. Finally, we discuss uncertainties, opportunities, and possible improvements in Sect. 6.

Site and instrumentation
The site used in this study is situated on two forested ridges in central Sweden (Fig. 1, left). It is located approximately 140 km from the Baltic sea coast, approximately 3 • below the polar circle. Within an area of 10 × 10 km, the site had two 59 m and five 100 m tall meteorological masts, each of which had cup anemometers and wind vanes at several heights (Table 1). The measurement periods were different from mast to mast, but all masts were operational between 23 February 2009 and 18 February 2010. Three different types of cup anemometers were used; NRG #40 anemometers, manufactured by NRG Systems, Inc., were used for the profile measurements on all masts. The two 59 m masts each had two top-mounted WindSensor (Risø) P2546A anemometers (Pedersen, 2004), and the top-mounted anemometers on the 100 m masts were Thies First Class anemometers.

Treatment of observational wind data
An initial screening of the observed wind data (10 min averages) showed significant inconsistencies in both wind speed and wind direction. The erroneous data were prevalent during winter and most likely caused by ice growth on the instruments. During most of these periods, the wind speed would have a constant near-zero value, signifying that the cup was completely frozen. However, there were also times when the wind speed of a particular anemometer was much lower than expected given other wind measurements, suggesting that the anemometer was able to turn, but at a lower rate than if ice free. The following data screening steps were applied to clean the data of ice-contaminated measurements.
1. Periods when the cup anemometer was clearly malfunctioning were removed, requiring that 0 ≤ U ≤ 50 m s −1 and I u < 1, where U is the wind speed and I u = σ u /U is the turbulence intensity, and σ u is the standard deviation of the wind speed.
2. Periods with constant wind speed were removed, requiring U i = U i±1 , where i denotes a 10 min block average.
3. Ice-affected data were removed by comparing pairs of cup anemometers in each mast, requiring the relative difference from at least three other anemometers on the same mast to be within 3σ of the ice-free mean relative difference, where the relative difference is defined as RD = (U i − U j )/U i , σ is the standard deviation of RD, and the ice-free period is defined as May to September.
Data passing these criteria were associated with a quality control (QC) value of 1. These steps removed between 8 and 26 % of the data. In addition to the ice-related issues, it was found that the NRG cup anemometer wind speeds were systematically lower than the top anemometers. We attribute this difference to both mast shadowing and a systematic small instrumental error, possibly related to the temperature response of the instrument. Finally, there was a slight offset between the two top anemometers related to wind direction. To get an accurate wind profile, both the NRG offset error and the top wind speed measurement were corrected as follows.
1. Which top cup anemometer data to use was determined based on wind direction.
2. The expected top wind speed was calculated from the profile cups by linear extrapolation from the two highest profile cups.
3. A correction factor was defined as the ratio between the actual top cup measured wind speed and the profile estimated value for each 10 min period.
4. The correction factor was applied to all wind speeds in the profile.
The final data cleaning step required that all cup anemometers had simultaneous values of QC = 1 to ensure that the wind distributions used in the cross predictions were generated from the same period across all anemometers. After applying the filtering steps, 8764 10 min mean wind speeds were available at each height and mast, which corresponds to approximately 1.5 months of data. The effect of the data filtering is shown in Sect. 5.
After cleaning the data, an "observed wind climate" was created for use in WAsP. The observed wind climate is a histogram of wind speeds for different wind direction sectors, i.e. the frequency of the wind for each wind speed and direction bin. For this study, a data discretization of 1 m s −1 for wind speed and 30 • for wind direction were used.

Roughness maps based on land use classes
A common approach of obtaining roughness information for wind flow modelling is through the conversion of land use classifications, typically derived from satellite data to a roughness length through the use of look-up tables. In this study, four different sources of land use classifications are investigated (Table 2). These four datasets are provided for download by WindPRO. All four of the land use datasets used in this study are derived from satellite measurements; however, there are some significant differences. The CORINE dataset has the highest spatial resolution at 100 m, but only covers Europe. This dataset has frequently been used in mesoscale meteorological modelling (Pindea et al., 2002) due to the high spatial resolution and large number of classification classes. Five of the 44 CORINE classes are associated with forests. The GLOBCOVER dataset has been used in the Global Wind Atlas (Mortensen et al., 2017), and while it has a lower resolution than the CORINE data, it is a more recent dataset, provides global coverage, and, important for this study, has eight forest classes out of the 23 different land cover classifications. The MODIS Vegetation Continuous Fields is not a true land use dataset. Instead it provides information about three components of the ground cover: percentage of tree cover, percentage of non-tree vegetation, and percentage of bare ground (Carroll et al., 2010). Finally, the Global Land Cover Characterization (GLCC) dataset provides a coarse 1000 m resolution land cover dataset, and is also the oldest of the datasets.
When downloading the data using WindPRO, standard roughness conversion tables are used. The values for classes in our model domain are found in Table 3. While some of the datasets have more forest classes than others, when translated to roughness there are at most two different values of roughness length for the different forest types. Additionally, the values for the different roughness classes are relatively similar across the different products, with the exception of CORINE, which has a category of peat bogs with very low roughness that is not in the other datasets. It is also important to note that the roughness lengths in these tables are not related to the local forest height, just the generic class description.

Airborne laser scans
Although using land use maps is attractive due to its simplicity, the subjective conversion of land use classes to roughness length can significantly affect the uncertainty of wind resource estimations, particularly for forested sites (Kelly and Jørgensen, 2017). Therefore, a more objective approach is presented in this section.

Estimating elevation and forest height
The area at and around the site was scanned by an airborne laser in 2013 from June to October, with a small part of the northern regions being collected in 2012. This timing is fortunate, as the scan period is quite close to the observational period of this study. The scans were made as part of a national survey (Lantmäteriet, 2016). Around 50 GB of data in the .las format, covering a 40 × 40 km area, was downloaded from the Swedish national database for this study. The reflection points in the point cloud were classified into ground, water, and vegetation points by the data provider. The data were first processed from the ALS point cloud to a terrain height map, a forest height map, and a map of all the water areas (lakes and rivers) at a 20 × 20 m 2 grid resolution, using the approach of Boudreault et al. (2015). The terrain height is classified as the median elevation z m of all ground or points, inside the 20 × 20 m 2 grid points, depending on if the grid cell was determined to be land or water. The forest height h was estimated as h = max(z i − z m ), where h is the forest height and z i indicates the vertical coordinate of all points i inside the grid point. Figure 1 shows the elevation (panel a) and tree heights (panel b) derived from the ALS processing, and an overview of the location of the seven meteorological masts used in this study. From the tree height map, it can be seen that the area is predominately forest, with h ranging from 5 to 35 m, but it also contains areas of clearings. This makes the site ideal for testing the different maps because the flow at meteorological masts is impacted by both changes in roughness and the geostrophic roughness (see Sect. 1). The focus of this study Wind Energ. Sci., 3, 353-370, 2018 www.wind-energ-sci.net/3/353/2018/  is on the forest parameters, but the terrain height will also change the forest's impact on the wind conditions, affecting the results. Therefore, the significant terrain variability at this site was desirable. The baseline 20 × 20 m tree height map was based on the maximum tree height. Since the lidar beams do not necessarily hit the tree top, this tree height estimate underestimates the maximum tree height. A relationship between scanning density and maximum tree height was shown in Boudreault et al. (2015). For the site considered here, the scanning density was between 0.25 and 1 points per square metre, leading to an underestimation of the tree height of up to 2 m. Nilsson et al. (2017) created maps of the mean tree height from ALSs using the the 95 percentile of points located above 1.5 m above ground, leading to a somewhat lower tree height. Compared to the uncertainties involved in the determination of roughness for the case when no tree height information is available, we consider this uncertainty to be of little importance.

An objective approach to estimate the roughness length and displacement height
The relatively simple conversion recommended by Garrat (1992) of z 0 = 0.1h was used for converting directly between forest height and roughness length. To simplify the maps and reduce the number of roughness changes, for compatibility with WAsP, the forest height map was rounded to the nearest 5 m height H . Non-forested areas were identified as regions with vegetation heights lower than 2.5 m, and had their roughness set to 0.1 m, while water areas were prescribed a roughness value of 0.0001 m. In summary, the roughness was parametrized as for water areas. (2) www.wind-energ-sci.net/3/353/2018/ Wind Energ. Sci., 3, 353-370, 2018 A sensitivity experiment was carried out using conversions of z 0 = 0.2H and z 0 = 0.05H to investigate how the roughness-to-tree height ratio impacts the modelling error.
In addition to the calculation of roughness length, the forest height was used to calculate a zero-plane displacement height d, which was added to the terrain height for a subset of the WAsP modelling runs. The relationship among the parameters used in this study again came from Garrat (1992): The approach described above is called the objective roughness approach (ORA).

Changing map resolution
To investigate the impact of map resolution and provide a fair comparison with the satellite-based maps, the 20 × 20 m resolution tree height maps were downgraded to 100, 300, 500, and 1000 m resolution roughness and displacement height maps using the following steps.
1. The forest height was recalculated as the arithmetic mean of the forest height pixels in the higher-resolution map.
2. Water areas were kept if the coarser resolution pixel consisted of more than 50 % high-resolution water pixels.
3. The new forest heights were converted to roughness length and displacement height as in Eqs. (2) and (3).
4. The resulting displacement height was linearly interpolated to a 60 × 60 m resolution terrain height map. However, in order to simulate a 20 × 20 m resolution case, we linearly interpolated the 60 × 60 m resolution terrain height map to a 20 × 20 m resolution before adding the displacement height map.
All resulting maps from these steps were exported as a Golden Software ASCII grid (.GRD) file for compatibility with the WAsP Map Editor. Throughout the rest of the paper, the model simulations using these maps will be labelled as ORA followed by the resolution in metres (e.g. ORA20).
For the simulations in which the displacement heights were added to the terrain height, the anemometer heights (Table 1) need to be adjusted. This is carried out by subtracting the mast displacement height from the anemometer height. The mast displacement height was calculated the same way as it is in WindPRO, the arithmetic mean of all grid points within a 50 m radius of each mast. Table 4 shows the mast displacement heights for each ORA resolution. ORA maps that use a displacement height are labelled with an additional "D". All masts are located in small clearings, leading to generally low mast displacement heights, but at coarser resolutions the clearings are less resolved.

Conversion of raster maps to vector maps for WAsP
The WAsP model currently only accepts vector-based roughness and elevation maps as input. Therefore, the rasterbased maps produced by the process described above needed to be vectorized. This process was carried out using version 11.20.5, release D of the WAsP Map Editor for all ORA resolutions coarser than 20 m. The Map Editor allows for the import of raster-based maps, which it then either contours, for elevation data, or converts to roughness change lines. For elevation data, raster data with a 60 m resolution were converted to 10 m equidistant contours. For the roughness change lines, the option to keep the same roughness classes as in the raster file was used. To create the roughness change map, the Map Editor outlines all pixels with roughness change lines that represent the inside and outside roughness value, skipping those where both sides of the line would be the same. After processing the data, they were examined via a geographical information system (GIS) to ensure that the vector data provided a reasonable representation of the elevation and roughness raster data. For converting the 20 m resolution maps, the Terrain workshop, version 1.1 release D was used because the Map Editor was not able to efficiently process the higher-resolution data.

WAsP set-up
WAsP version 11.6 was used in this study to simulate the wind resource at the different masts. In WAsP, the effect of changes in roughness length on the wind speed around the site are modelled using the internal boundary layer zoominggrid model (Troen and Petersen, 1989;Sempreviva et al., 1990), while the speed-up effects from terrain elevation are simulated using the spectral model described in Troen and Petersen (1989) and Troen (1990). As input, WAsP requires an observed wind climate and vector maps of elevation and roughness. WAsP uses the following modelling chain to simulate the wind climate at one position, horizontal and vertical, from an observed wind climate at another position. This modelling chain is extensively described in Troen and Petersen (1989) and Kelly and Troen (2016, Eqs. 13-15). Here we briefly summarize the steps.
1. Fit a Weibull distribution to the observed wind climate for each sector, preserving the third moment of wind speed and the probability to exceed the mean wind speed.
2. Calculate the background wind profile using the geostrophic roughness and the observed wind speed, from which the effects of the roughness changes, terrain height, and atmospheric stability at the observational site have been removed.
3. Calculate the wind speed (generalized wind climate) for predefined heights and roughness lengths using the geostrophic drag law in combination with the background wind profile.
4. Apply the effects of geostrophic roughness, terrain height, roughness changes, and stability using this generalized wind climate for the predicted site.
The predicted wind climate can be quite sensitive to the choice of the standard heights and roughness lengths used in step 3. For this study, the predefined heights were set to 3, 10, 30, 60, and 120 m above the surface here, and the roughness lengths were set to 0.0, 0.1, 0.4, 1.0, 3.0 m, which covers all possible roughness length values that occur in the maps. In the WAsP model, large water bodies are required to have a roughness length of 0.0 m, which is then internally converted to a value of z 0 = 0.0002 m.
The impact of varying topography, roughness changes, and orography is modelled as correction factors that are applied to the scale parameter (A) of the Weibull distribution of the wind climate in steps 1 and 4 of the model chain. For computational efficiency, WAsP filters the roughness changes to www.wind-energ-sci.net/3/353/2018/ Wind Energ. Sci., 3, 353-370, 2018 include only those that have the most significant impact on the wind speed. This is carried out by creating a distanceweighted roughness length (ln z 0w ), which is found by multiplying log-transformed roughness values with the exponentially weighted distance from the mast to the roughness changes x k , where x d is a decay length, which is currently set to 10 km. Then the n items in the array with all transformed roughnesses and distances are found by fitting a step function that stops when either the maximum number of allowed steps n max is reached or the residual variance RMS max is below a specified threshold. More details about the algorithm can be found in Sect. 8.3 in Troen and Petersen (1989). The default values for n max and RMS max are 10 and 0.3, respectively. From this filtered set of roughness change lines, the internal boundary layer equations described in Sempreviva et al. (1990) and Floors et al. (2011) are applied to the reduced arrays to compute a speed-up factor in each sector.
To calculate the geostrophic wind, WAsP uses a corrected logarithmic wind profile, for which the correction factors for internal boundary layer and orographic effects have been removed, in combination with the geostrophic drag law. The roughness length that is used in these relations is the socalled geostrophic roughness length z 0G , which is computed sector-wise from the mean of ln(z 0 ) and with a distance weight similar to Eq. (4): where z 01 is the roughness length at the mast, k = 2, . . . , N − 1 denotes the kth roughness change from the mast location, and N is the last roughness change on the map. In addition to the standard inputs, WAsP uses the longterm distribution of heat fluxes at the site to model the effect of atmospheric stability. Since no heat flux observations were taken at the site, the default heat flux distribution, with mean and the standard deviation of −40 and 100 W m −2 , respectively, was used over land. Over water, these values were set to −8 and 30 W m −2 .
The performance of the different topographic maps was evaluated using cross predictions. A cross prediction is defined as the prediction of the flow from one observed wind climate, a specific mast and height, to another observed position, either another height on the same mast or an observed height on another mast. Cross predictions were made from all four heights at six of the seven masts, but at mast 6 the 80.7 m height was excluded due to the limited number of data with QC = 1. After excluding self-predictions, i.e. predictions and inputs at the same mast and height, there were a total of 702 different combinations. The relative errors for each cross prediction were computed as a percentage from the observed (obs) and modelled data (mod) as δ = 100(mod − obs)/obs for both wind speed (δU ) and power density (δP ). It is important to include power density in the evaluation since the production of wind turbines is determined by the available power. The total power density is calculated by summation of the frequency-weighted third moment of the Weibull distribution from each sector of the total number of sectors D, where ρ is a reference air density (here 1.225 kg m −3 ), f is the frequency of occurrence, and k is the shape parameter of the Weibull distribution.

Results
The results are presented in three subsections. First, the ORA roughness lengths will be compared with the roughness lengths from the land-use-based datasets. Second, the results of the wind data filtering algorithms will be shown. Finally, the results from the WAsP cross predictions are shown. forest edges and clearings are positioned differently across the different datasets.

Roughness maps
In the ORA maps the average roughness of the map varies between 1.42 and 1.46 m. The logarithmic average roughness value over these maps changed from 1.02 m (20 × 20 m) to 1.39 m (1000 × 1000 m) resolution. The impact of resolution is discussed further in Sect. 6.
It is clear from Fig. 2 that the roughness lengths from the ORA maps are larger than from the land-use-based maps. Additionally, the ORA data, in part due to the higher roughness values, have significantly more roughness changes. For example, CORINE only has two forest roughness lengths, 0.4 and 0.5 m, while the ORA data represent forest roughness lengths in six different bins from 0.5 to 3.0 m.
These features can be seen more clearly in the histograms of the roughness lengths for each map (Fig. 3). For example, the roughness lengths of the grid cells in the ORA maps are higher on average and are spread over more bins than the satellite-based maps. The GLCC dataset is dominated by a single land use class with z 0 = 0.5 m. The MODIS data have slightly more detail and grid cells with z 0 = 0.3 m are most frequent. The CORINE data have the highest resolution of www.wind-energ-sci.net/3/353/2018/ Wind Energ. Sci., 3, 353-370, 2018 the land-use-based products and have most grid cells with z 0 = 0.5 m. Note that the CORINE and ORA20 maps include more variation in the roughness, which causes the geometric mean to be smaller than the arithmetic mean.

Wind distributions and profiles
To demonstrate the impact of the quality control that was described in Sect. 2, three distributions of wind speeds at different stages of the quality control process are shown for the cup anemometer located at mast 1, 59 m (Fig. 4). Due to the confidentiality of the data, a normalized wind speedÛ was computed by dividing each 10 min measurement by the mean wind speed of the filtered data. The distribution before any QC filtering (Fig. 4a) shows a high frequency of measurements with very low wind speeds (< 0.2) as well as an increase in wind speed occurrences between 0.2 and 1.0, when compared to the distribution in Fig. 4b that includes only data with QC = 1. These differences likely reflect the icing of the cup anemometers. Because of the large number of nearzero wind speeds, the Weibull distribution from the WAsP method in Fig. 4a gives a rather poor fit. The fit was improved for the distributions including only data with QC = 1 (Fig. 4b), which include 59 % of the original dataset, and requiring that QC = 1 for all cup anemometers simultaneously (Fig. 4c), which retains only 16 % of the original dataset. The mean wind speed for the three distributions were 0.85, 0.97, and 1, highlighting that although there was a large change in dataset size when requiring QC = 1 for all cups simultaneously, there was a small change in mean wind speed compared to QC = 1 for a single location. This relationship was found for all masts and heights. All 27 anemometers showed an increase between 1 and 6 % when demanding that QC = 1 for all cup anemometers simultaneously, which shows that the wind climate in the final dataset used for the validation is representative for the site. Figure 5 shows the vertical profiles of mean wind speed for all anemometer positions using the observed wind climate from mast 1 at 59 m. The data have been normalized with the mean wind speed from the observed position. In addition to the land-use-based simulations, the ORA runs include the highest-resolution maps both with displacement heights included (ORA20D) and without (ORA20). The impact of introducing a displacement height can be observed near the canopy, where the wind speed is much lower than when a displacement height is not included. At masts 1, 3, and 7 the ORA20D simulations are closer to the observations than without displacement, while they are comparable at the other masts. At higher heights, the differences in mean wind speed among different tests are smaller. The worst predictions occur at mast 6; however, this site is at a higher hill than the other masts, so this difference could be due to orographic speed-ups.

WAsP sensitivity to different roughness maps
The first test performed was to investigate the sensitivity of the ORA method to different tree height to roughness length conversions in Eq. (2). In this sensitivity test, simulations were performed using maps where z 0 = H /5, z 0 = H /10, and z 0 = H /20 with and without applying a displacement length. For each pair of observed histograms, the histogram at the source location was used to predict the wind distribution at the target location and compared to the observed histogram at the target location (see Sect. 4). The resulting mean absolute error from all 702 cross predictions, |δU |, is shown in Fig. 6, for these tests as well as a sensitivity test of the CORINE land-use-based roughness.
For the roughness maps without a displacement height, it can be seen that increased roughness lengths result in lower errors, i.e. H /20 has the highest error and H /5, which corresponds to a mean roughness length of more than 2 m, has the lowest. When a displacement height is applied, better representing the location, the H /5 test has the largest errors, while the tests using H /10 and H /20 maps have very similar |δU |. When looking at power density, it was found that the ORA100D simulations with H /10 and H /20 had a |δP | of 9.48 and 10.74 %, respectively. This study shows that z 0 = H /10 is a reasonable modelling choice.
After finding that higher roughness lengths lead to better results, and noting that the roughness maps based on land use classes (GLCC1000, MODIS500, CORINE100, and GLOB300) are characterized by much lower roughness lengths than the ORA maps (Fig. 3), one could hypothesize that increasing z 0 in the satellite-based products to a more realistic value would improve the results. The CORINE100 data were chosen for this since it could be seen that the CORINE100 map had a realistic representation of the land use around the site (Fig. 2).
To test this, the roughness length of the forest classes (i.e. areas with z 0 = 0.4 m and z 0 = 0.5 m) in the CORINE100 map were increased by a factor of 3 to 1.2 and 1.5 m, respectively, which is approximately the roughness length observed in maps using the ORA approach (see Fig. 3). The simulation using this roughness map is labelled with "high z 0 " and had an error of |δU | = 3.5 %, which is closer to the ORA100 test's error of |δU | = 3.4 % than the standard CORINE100 test with an error of |δU | = 3.8 %. Therefore, the performance of the land-use-based maps could potentially be improved for this site by increasing the roughness length. However, the factor of 3 was chosen based on knowledge of the roughness length from the ORA maps; thus in the absence of tree height measurements it could be difficult to set the correct value.
The next sensitivity test was designed to investigate the effect of the data filtering on the roughness change model. This was achieved by changing the limit of RMS max to ≈ 0 and varying n max from 0 to 10, which causes the model to include only n max roughness changes after filtering. Note that    Figure 4. Histogram of the observed all-sector normalized wind speedÛ (see text) at mast 1 at 59 m a.g.l. without any filtering (a), after selecting data with QC = 1 with 59 % of the data (b), and during the times when all masts simultaneously had QC = 1 with 16 % of the data (c). The red line denotes the Weibull distribution that WAsP fits for this histogram.
using 0 roughness changes disables the internal boundary layer model, which means that all observed differences are due to a different z 0G (see Sect. 4). To isolate the impact of taking more roughness changes into account, all simulations were normalized by the |δU | from simulations with 0 roughness changes. Figure 7 shows the normalized |δU | for different numbers of roughness changes. It can be seen that for most roughness maps, the model errors generally decrease from 0 to 2 n max , and for larger n max the ratio is more constant. This indicates that in this case only two roughness changes contribute significantly to WAsP modelling at each mast. This is consistent with the best practice for hand digitizing roughness maps, which states that only the most important roughness changes need to be digitized (Mortensen, 2016). It should be noted that n max is likely site specific and that using n max > 2 could still give significant improvements at other sites.
By comparing the different ORA resolutions, the impact of resolution on the roughness change model can be examined. It was found that the simulation using the ORA20 map had the largest reduction and |δU |/|δU 0 | ≈ 0.86, while the coarser-resolution ORA1000 map did not have as large of a reduction in error ratio, |δU |/|δU 0 | ≈ 0.94 using > 2 roughness changes. This shows that resolution is important for roughness change modelling, which was expected since the higher-resolution maps can better resolve the precise location of a roughness change. Interestingly, the GLOB300 simulations have |δU |/|δU 0 | > 1, indicating that the map with n max = 0 performs best. This means that the model results can lead to larger errors when including roughness changes that do not represent the modelling area well.
To further illustrate the effect of using different maps, we can plot the |δU | as a function of z 0G (Fig. 8), where z 0G is the all-sector geostrophic roughness length, which was calcu- lated as the geometric mean of z 0G from each sector. Since there are no displacement heights available for the land-usebased maps, the displacement height was not used in this comparison. As expected, the satellite-based products have a much lower z 0G than the ORA maps. There is also a relationship between the model error and z 0G , in that the MODIS500 map has the highest |δU | of ≈ 4.4 %, with the lowest z 0G , whereas the ORA1000 map has the highest z 0G and the lowest |δU | of ≈ 3.4 %.
The effect of including a displacement height for the ORA maps at different resolutions is shown in Fig. 9. The roughness change model is used with the default settings (see Sect. 4); thus these results are impacted by both the roughness changes and the different values of z 0G . The ORA20D and ORA100D have a |δU | that is ≈ 0.2 % lower than the ORA20 and ORA100, respectively; they therefore benefit from including a displacement height. At coarser resolutions, adding the displacement height increases |δU | by ≈ 0.5 %. The displacement height at all masts using the ORA1000 map is significantly higher than for the ORA20 map (Table 4). This suggests that the displacement height for the coarser resolutions is likely too high, because the forest clearing around the mast is not resolved. Therefore, it is recommended that the displacement height correction should only be performed when the resolution of the map is sufficient to represent the effect of forest clearings near the mast.

Summary of using different roughness maps
In this section, the overall performance of WAsP simulations using satellite-based maps and the proposed ORA approach is summarized. Thus far, only |δU | was used for comparing the model results since the comparisons were largely investigations of the model behaviour. However, for an annual energy production estimation, it is more relevant to quantify map-related errors on |δP | (Eq. 6) to assess the potential performance of the wind turbine. Therefore, to evaluate the overall performance of the modelling with different maps, both |δU | and |δP | will be used. Figure 10 shows the distribution of |δP | for each of the cross predictions, binned in three categories for all four landuse-based maps and the highest-resolution ORA map with a displacement height. For this analysis, the cross predictions are broken into two types. Horizontal predictions (Fig. 10a) are those for which the observed wind climate and the predicted wind climate are from different masts, while the vertical cross predictions are those for which they come from the same mast (Fig. 10b).
Ideally, we want as many cross predictions with |δP | < 12.5 % as possible. For the horizontal cross predictions, there are ≈ 20 % more predictions in this category for the ORA20D run than for the land-use-based maps. Conversely, the number of very high errors, |δP | > 25 %, is greatly decreased for the horizontal cross predictions when using the ORA20D map, from ≈ 100 when using the CORINE100 and GLCC1000 map to ≈ 20 when using the ORA20D map.
For the vertical cross predictions (Fig. 10b), the errors are much smaller than for horizontal extrapolations because there are the same topographic errors for the input and the predicted anemometers. However, as in the horizontal cross predictions, the ORA20D run has fewer high-error predictions than the land-use-based roughness map; the ORA20D simulation had only 1 error > 12.5 %, whereas for all other land use products they occurred at least 11 times.
Even when comparing maps with the same resolution (not shown), i.e. the CORINE100 with the ORA100 simulation, the number of cross predictions with errors larger than 25 % decreases from 75 to 44, a 41 % decrease. Using the ORA300 compared to the GLOB300 simulation decreases these errors www.wind-energ-sci.net/3/353/2018/ Wind Energ. Sci., 3, 353-370, 2018  by 51 % and using the ORA500 compared to the MODIS500 decreased them by 40 %. Table 5 provides summary error statistics of the results from the 702 cross predictions for all land use, ORA, and ORA-with-displacement simulations. In addition to the mean absolute error of δU and δP , which has been shown previously, the mean bias and root-mean-square error (RMSE) have been included. The RMSE penalizes high errors more than the mean absolute error, and therefore identifies simulations that have a narrow error distribution.
All ORA-based simulations have lower |δP | than any of the simulations using land-use-based maps. The same is true for the RMSE of δP , where the improvement of the simulations with the ORA roughness maps is even larger due to the reduced number of large errors as seen in Fig. 10.
Furthermore, it was found that the ORA100D simulation has the lowest |δP |, RMS of δU , and RMS of δP , whereas the ORA20D has the lowest errors in |δU |. The ORA20D map results in lower errors due to a more accurate roughness change description (Fig. 7), and due to the application of a displacement height (Fig. 9). However, the z 0G was lower than that used for better-performing maps in the noroughness change tests (Fig. 8). From these results it is clear that the ORA maps provide better results, due to the higher z 0G and because of the increased detail of the forest structure. For this site, a best-performing map could likely be created by using the 20 m resolution map with a slightly higher z 0 .
The land-use-based land cover simulations perform similarly despite different resolutions and different descriptions of the land use. For example, the GLOB300 and CORINE100 Wind Energ. Sci., 3, 353-370, 2018 www.wind-energ-sci.net/3/353/2018/ simulations both have a rather high RMSE for both U and P , despite having a high resolution. This highlights that resolution is not a panacea but only improves performance when the surface is accurately represented. For example, the GLCC1000 simulation has almost no roughness changes (see Fig. 2), but has the lowest RMSE in U and P out of the satellite-based products. Also, we noted that the CORINE100 with a high z 0 (see Sect. 5.3) reduced |δP | to 12.63 % and the RMSE of δP to 15.76 % (not shown); however, both are still significantly higher than all of the ORA runs. This shows that tweaking z 0 of certain land use classes in existing satellite-based maps will only partially reduce the model errors.
In Table 5, the mean bias of U is very close to zero for all model runs. The ORA1000 simulation has the lowest mean bias of U . One would expect a difference between runs with a very low and very high roughness, for example the CORINE100 and ORA100 runs. However, it can be seen that for both of these runs, the bias is close to zero. This is likely because the cross prediction includes both upwards and downwards predictions. This results in errors that will cancel each other out; e.g. a too low prescribed roughness results in an over-prediction for an upward extrapolation but an under-prediction for a downward extrapolation.

Discussion
The presented results for the mean prediction error in Table 5 show a complex interdependency of roughness magnitude (Figs. 6 and 8), imprecision of roughness lines due to either limited map resolution or a coarse land use classification (Fig. 7), and inclusion (or not) of displacement height (Figs. 5 and 9). Nevertheless, the ORA approach showed consistently improved results relative to the roughness maps based on land use classification for the studied site. This demonstrated robustness of the ORA approach should make it attractive to use for turbine siting. Given the relatively small processing time needed to make the maps, and the fact that all maps could be loaded into WAsP and predictions finalized within a few minutes, the approach is promising. At this site, it turned out to be better to use a land use map with hardly any roughness change lines than to use products with more detailed, but incorrect, information (e.g. GLCC1000 versus GLOB300) when comparing only the commercially available maps (Table 5).

Ideas for further improvements
As argued by Jackson (1981) and Raupach (1994), forest density and forest height are likely to influence the optimal ratios of z 0 /H and d/H . Since forest density in terms of plant area index can also be estimated from airborne lidar data (Boudreault et al., 2015), it would be interesting to investigate a more refined roughness classification based on this parameter. It was also demonstrated that raising the roughness value in the CORINE data reduced the mean prediction error compared to the original roughness conversion (Fig. 6), indicating that this potential "quick fix" could result in significantly reduced errors. Without access to forest height information, it will however be hard to justify the exact level to which the roughness should be adjusted. Since imprecise roughness change lines can lead to an increase in prediction error, as seen when using GLOB300 data (Fig. 7), a simpler way of reducing the error could be to smooth the map and increase the geostrophic roughness in order to achieve lower error, as demonstrated for the ORA maps in Fig. 8. The results in Fig. 8 reinforce the recommendation to use high values of roughness length for forested areas, even for landscapes where all mast observations are affected by the forest.
When downgrading the resolution of the ORA maps, the arithmetic average of the roughness remained nearly identical because the starting point for all the ORA maps was the tree height, whereas the presence of clearings in the high-resolution ORA maps reduced the geostrophic roughness (Fig. 8) since the WAsP method uses the geometric average (see, e.g. Taylor, 1987, for a motivation and more background). In forested areas, the use of geometric averaging could be particularly problematic since -in realitythe presence of clearings, which have low roughness lengths, tend to increase the overall turbulence levels and increase the aggregated roughness of the landscape (Lopes et al., 2015). It should also be noted that the log-space averaging approach has limitations for low vegetation since it omits the significant non-equilibrium effects of roughness aggregation (Hasager and Jensen, 1999).
Because of the log-averaged aggregated geostrophic roughness, it would be of interest to systematically investigate the roughness value chosen for the clearings. The value of z 0 = 0.1 m, used in this study, was assessed to be low for a newly cleared forest area but reasonable to high for lowvegetation wetland. If a larger roughness length was used for the clearings, it would be expected that z 0G would also be increased, which would reduce the error associated with the geostrophic roughness closer to the level of the ORA1000 map. There is, however, a risk that increasing the value could affect the filtering of significant roughness changes, which could thereby increase the related error.
6.2 Does a better map make a significant difference?
The reduction in mean prediction error for both the mean wind speed and power density is only a few percent. However, as demonstrated in Fig. 10, this reduction reflects changes over a wide distribution of cross-prediction errors. We see the demonstrated reduction of the large errors as one of the main advantages of the ORA method since such large errors in the predicted wind climate could result in significant costs when estimating the performance of a wind turbine throughout its lifetime. Whereas large errors could re-sult in higher lifetime costs of turbine operation at any site, Enevoldsen (2017) revealed that wind turbines located in forested areas are more prone to experience fatigue loads, vibration errors, and lower annual energy productions than other sites. Therefore, an improved wind climate prediction could be of high consequence for the forested landscape.
With this in mind, it is recommended that effort is put into finding accurate and updated information on forest height and structure when evaluating an optimal location for wind energy exploitation in forested areas. For the reasons stated above, we also recommend investigating the sensitivity of the observed wind climate to the relationship between z 0 and H since the ratio used here may not be optimal for other forests.

Availability of tree height information data
To the authors' knowledge, databases of freely accessible ALS data already exist for Denmark, Finland, and England. In other countries, such as Sweden, data are accessible at reduced cost (compared to a new campaign) for international researchers and companies, whereas they are freely accessible for national researchers.
The processing of lidar data to tree heights is relatively simple, and can be performed on a normal desktop computer for small areas. In addition, future drone-mounted lidars and satellite-derived products (Stone et al., 2016) have the potential to provide high-accuracy tree height or land cover information for relatively low cost. This will enable better access to such data and provide for the ability to assess changes over time.

Conclusions
This study has demonstrated a new way of introducing highquality roughness maps derived from ALS scans of forest areas into the WAsP and WindPRO software. The steps involved can be reproduced by siting engineers with access to a reasonably powerful computer. By examining wind conditions at seven meteorological masts in a coniferous forested area in Sweden, it can be concluded that simulations using the tree-height-based roughness maps resulted in a closer agreement with observational data when compared to those that used standard land-use-based online sources. The lowest error was achieved when using the highest detail maps, 20 and 100 m resolutions, and including the displacement height via the terrain height (Fig. 9). A sensitivity test found that using z 0 = H /10 with a displacement height provided the best results for converting tree height to roughness length when the displacement was used, but the conversion of z 0 = H /5 performed better without a displacement height. In addition, while roughness length maps with a high resolution benefited significantly from applying a displacement height, it did not improve the results when using maps with a coarser resolution than 100 m. The highest-resolution ORA data provided the best results, but an improvement over the land-use-based map simulations was also achieved when lowering the spatial resolution to 1000 m. It can therefore be suggested that siting engineers and practitioners of wind resource assessment in forested areas should collect appropriate site data that describe forest height, type, and density.
Code and data availability. The wind data used in this study are confidential and unfortunately not made public.
Author contributions. ED and JA processed the tree height and cup anemometer data. ED, PE, and ND made the maps. RF performed the modelling and made the plots. RF, ED, ND, PE, and JA wrote the paper.
Competing interests. Rogier Floors, Neil Davis, and Ebba Dellwik work at DTU Wind Energy, where WAsP is maintained, developed, and sold.