Articles | Volume 6, issue 6
Wind Energ. Sci., 6, 1379–1400, 2021
Wind Energ. Sci., 6, 1379–1400, 2021

Research article 04 Nov 2021

Research article | 04 Nov 2021

Satellite-based estimation of roughness lengths and displacement heights for wind resource modelling

Satellite-based estimation of roughness lengths and displacement heights for wind resource modelling
Rogier Floors1, Merete Badger1, Ib Troen1, Kenneth Grogan2, and Finn-Hendrik Permien3 Rogier Floors et al.
  • 1DTU Wind Energy, Risø Campus, Technical University of Denmark, Frederiksborgvej 399, 4000 Roskilde, Denmark
  • 2DHI GRAS A/S, Agern Alle 5, 2970 Hørsholm, Denmark
  • 3Siemens Gamesa Renewable Energy A/S, Borupvej 16, 7330 Brande, Denmark

Correspondence: Rogier Floors (


Wind turbines in northern Europe are frequently placed in forests, which sets new wind resource modelling requirements. Accurate mapping of the land surface can be challenging at forested sites due to sudden transitions between patches with very different aerodynamic properties, e.g. tall trees, clearings, and lakes. Tree growth and deforestation can lead to temporal changes of the forest. Global or pan-European land cover data sets fail to resolve these forest properties, aerial lidar campaigns are costly and infrequent, and manual digitization is labour-intensive and subjective. Here, we investigate the potential of using satellite observations to characterize the land surface in connection with wind energy flow modelling using the Wind Atlas Analysis and Application Program (WAsP). Collocated maps of the land cover, tree height, and leaf area index (LAI) are generated based on observations from the Sentinel-1 and Sentinel-2 missions combined with the Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2). Three different forest canopy models are applied to convert these maps to roughness lengths and displacement heights. We introduce new functionalities for WAsP, which can process detailed land cover maps containing both roughness lengths and displacement heights. Validation is carried out through cross-prediction analyses at eight well-instrumented sites in various landscapes where measurements at one mast are used to predict wind resources at another nearby mast. The use of novel satellite-based input maps in combination with a canopy model leads to lower cross-prediction errors of the wind power density (rms = 10.9 %–11.2 %) than using standard global or pan-European land cover data sets for land surface parameterization (rms = 14.2 %–19.7 %). Differences in the cross-predictions resulting from the three different canopy models are minor. The satellite-based maps show cross-prediction errors close to those obtained from aerial lidar scans and manually digitized maps. The results demonstrate the value of using detailed satellite-based land cover maps for micro-scale flow modelling.

1 Introduction

Wind turbines on land represent a cost-competitive source of renewable energy (Global Wind Energy Council2019). More than 95 % of the global installed wind energy capacity of ≈733 GW in 2020 is installed on land (, last access: 21 April 2021). In northern Europe's temperate climates, a vast amount of the land surface is covered by forest. The exploitation of wind power within such forests has become more widespread as the hub height of modern wind turbines has exceeded the forest height (Enevoldsen2016).

Wind resource assessment is typically performed with linearized modelling in a wind farm siting tool (e.g. WAsP or WindPRO), which contains several submodels to predict the flow based on an input wind histogram (Troen and Petersen1989). Based on Monin–Obukhov similarity theory (Eq. 1), the flow model can be used to predict the wind speed, U, for any height, z, above the ground:

(1) U = u * κ ln z - d z 0 - Ψ m ,

where u* is the friction velocity, κ is the von Kármán constant (=0.4), z0 is the aerodynamic roughness length, 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). In addition to z0, the zero-plane displacement height, d, is traditionally used at forested sites to account for the canopy, forcing the mean flow to be displaced upwards over forests (Thom1971). A dense forest appears more smooth (i.e. a lower z0) and has a larger d than a sparse forest with clearings (Shaw and Pereira1982).

Values of z0 and d can be assessed through visual inspection of the forest at a given site combined with digital maps. In practice, a background value for z0 and d is set, and adjustments are made for specific areas where the roughness and displacement height differ from the background values. Manual assessments are subjective and time-consuming, and they can lead to a high level of uncertainty of the estimated wind resource (Kelly and Jørgensen2017).

Fully automated assessments of z0 and d can be achieved based on global or regional land cover data sets derived from satellite observations. Each land cover class is assigned a value of z0 and d via a land cover table (Jancewicz and Szymanowski2017; Floors et al.2018). State-of-the-art flow modelling tools offer embedded access to such land cover maps and to the associated roughness translation tables, which the user may modify. Due to the coarse spatial resolution of global land cover data sets (grid spacing of 100–1000 m), the finer-scale variability within a forest, such as smaller clearings on the order of 10–100 m, is not resolved. Further, the available land cover-to-roughness translation tables may not be fully representative for the site in question due to the data sets' global nature. For example, commonly used tables (see Appendix A) show very low roughness lengths for forest classes, which might not be accurate for areas with heterogeneously forested land cover (Floors et al.2018).

Enevoldsen (2017) and Floors et al. (2018) have demonstrated that tree heights and forest densities retrieved from aerial lidar scans can be used to parameterize z0 and d over the forest. This approach is more physical than the ad hoc assignment using land cover data sets. It sets new requirements to the flow modelling tools used for wind energy siting because (i) a canopy model is needed to estimate z0 and d from the tree height and density, and (ii) the data processing routines in the flow modelling tool need to be more efficient due to the finer resolution of z0 and d that result from lidar-derived canopy measurements. Several countries in northern Europe have released national aerial lidar scans, and in addition, dedicated aerial lidar campaigns may be carried out to obtain location-specific data sets for wind farm planning. Lidar campaigns are rare due to the high cost of the instrumentation and deployment and the time-intensive data processing requirements. The temporal frequency of such observations is therefore low.

A wealth of new satellite observations with unprecedented spectral properties and spatial and temporal resolutions have become available, e.g. through Copernicus (, last access: 27 September 2021). Forest monitoring is a key objective of several new missions because information on deforestation and forest degradation is important in connection with climate change mitigation. However, key metrics for wind-resource assessment such as the forestry canopy height are still missing or only available at a low spatial resolution, but they can be derived through post-processing of the data from available sensors. Here, we hypothesize that such post-processed values of the forest canopy height and density retrieved from satellites at a high spatial resolution can also be used to estimate the wind resource for a site with the same accuracy as aerial lidar scans but at a lower cost.

This paper uses tree heights and densities retrieved from satellite observations to derive roughness lengths and displacement heights that are then used to estimate the wind energy resource available at forested sites. Based on wind observations from meteorological masts at eight sites worldwide, we evaluate the prediction errors on the wind power density and the mean wind speed, and we compare them with predictions based on global land cover maps, aerial lidar scans, and manual digitization for the parameterization of z0 and d.

2 Background

2.1 Forest parameters from satellites

Firstly, a goal in this paper is to derive tree heights from satellites. The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) carries the Advanced Topographic Laser Altimeter System (ATLAS), used to observe the land surface height with great precision. Amongst many other parameters, the mission delivers global forest canopy heights (Neuenschwander and Pitts2019), which are well correlated with canopy heights from aerial lidar scans (Li et al.2020) and field measurements (Huang et al.2019). A mission with similar capabilities to ICESat-2 is the Global Ecosystem Dynamics Investigation (GEDI). However, GEDI only collects data up to 52 north and south of the Equator and therefore does not offer global coverage. Although representing a major advancement in estimating 3D forest structure, these current spaceborne laser observations are restricted to very narrow footprints, which are insufficient to map canopy heights over larger areas. Recent studies have shown that the laser-derived canopy heights can be extrapolated through different combinations with other satellite observations and machine learning techniques (Csillik et al.2020; Fagua et al.2019). Typically, the canopy heights from laser measurements are extrapolated using textural information from active microwave sensors (e.g. Sentinel-1 or ALOS PALSAR) and multispectral information from passive sensors (e.g. Landsat or Sentinel-2) (Li et al.2020).

Secondly, in connection with flow modelling for wind energy, the leaf area index (LAI) can be used as a proxy for the forest density (Raupach1994). LAI is defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as one-half the total needle surface area per unit ground area in coniferous canopies (Chen and Black1992). Several space-borne sensors, operating in the visible and near-infrared range, monitor LAI and other vegetation properties daily. For example, the Terra and Aqua satellites each carry a Moderate Resolution Imaging Spectroradiometer (MODIS). From the two instruments in combination, 4 d composites of LAI are generated routinely. This product has a relatively coarse resolution and a pixel size of 500 m. Guzinski and Nieto (2019) have developed a method for downscaling of LAI estimates.

2.2 Forest roughness models

Different models can be used to estimate z0 and d from forest canopy heights and LAI. Here, we consider the objective roughness approach (ORA), the Raupach model, and the scalar distribution (SCADIS) model.

2.2.1 The objective roughness approach (ORA)

The relation between the canopy height, h, and z0 and d has been recognized and discussed by many authors (Thom1971; De Bruin and Moore1985; Enevoldsen2017). Because z0 is usually proportional to h, an easy way to obtain z0 and d is by relating them linearly to the canopy height,

(2) z 0 = c 1 h ,


(3) d = c 2 h .

Maps based on lidar scans using the constant values c1=0.1 and c2=2/3 were shown in Floors et al. (2018) to reduce the risk of making large errors (>25 %) in predicted power density by 40 %–50 % compared to the best land cover database at a forested site.

2.2.2 The Raupach model

Raupach (1992) developed a model to predict the bulk drag coefficient over a rough surface with a given canopy height. The main model parameter resulting from his analysis is the frontal area index λ, which for isotropically oriented elements is given by

(4) λ = 0.5 LAI .

Raupach (1994) discussed some simplifications to the original model and suggested

(5) d h = 1 - b ,


(6) b = 1 - exp ( - a ) a


(7) a = 2 c d 1 λ .

cd1 was experimentally found to be equal to 7.5. Finally, the roughness length is estimated as

(8) z 0 h = b exp - κ min C S + C R λ , c max - Ψ h .

Raupach (1994) suggested based on empirical evidence that CS=0.003, CR=0.3, cmax=0.3, and Ψh=0.193.

2.2.3 The SCADIS model

A criticism of the Raupach (1994) model is that some of the constants are highly dependent on the structure of the canopy and may not be universally applicable. To address this, one can use a one-dimensional version of a kω turbulence model, which depends on the leaf area density profile within the canopy (Sogachev et al.2002; Sogachev and Panferov2006). An additional advantage of this model is that it can estimate z0 and d in non-neutral atmospheric conditions.

The leaf area density (LAD) profile is frequently described (e.g. Meyers and Tha Paw U1986) by a beta probability density function,

(9) LAD = LAI ( z / h ) α - 1 ( 1 - z / h ) β - 1 ,

where α=9 and β=3. The chosen constants are representative of a temperate pine forest, and we use them throughout this work due to the absence of information on the canopy profile from satellite data. Sogachev et al. (2017) validated the SCADIS model for flow over a 3D area with forested terrain and found that simulations that explicitly resolve the drag of the canopy and those that use an effective z0 compare well.

2.3 Flow modelling and cross-prediction

Maps of z0 and d can be used in combination with terrain elevation maps as input to wind energy flow models. The WAsP methodology (Troen and Petersen1989) can be applied to analyse the observed wind climate (a histogram of wind speeds for each wind direction sector) at a particular mast location and height and predict the wind climate at a nearby location, assuming the same large-scale atmospheric forcing. The effects of land surface roughness and orography are first removed from the observed wind climate. Equation (1) and the geostrophic drag law (Blackadar and Tennekes1968) are then used to estimate the geostrophic wind speed distribution, which is assumed to be valid some kilometres away from the mast. The geostrophic wind speed distributions can be transformed to Weibull distributions over idealized flat terrain at a number of specified heights and z0 (see Sect. 8.7 in Troen and Petersen1989). The result of this process is a so-called generalized wind climate (GWC), which is characterized by the Weibull parameters A and k and the probability density per wind direction sector. The GWC can be used to predict the wind resource near the measurement site by adding local roughness and orographic speed-up factors. The process of predicting the wind climate from one position to another is called a cross-prediction. The accuracy of such predictions depends strongly on the quality of the roughness and elevation maps used as input (Floors et al.2018).

3 Sites and data

In this section, we present eight sites that are considered in this study and the different types of input data which are needed for cross-prediction analyses at the sites.

3.1 Measurement sites

Eight measurement sites are selected for cross-prediction analyses (Table 1). The sites represent vastly different land cover types and complexities. The sites Ryningsnäs, Finland, and Sweden are located in forests, where the displacement height becomes important, and the uncertainty in wind resource assessment is known to be high (Kelly and Jørgensen2017). The site near Mérida has low (h≈5 m) and irregular forest in all directions. The sites Cuauhtémoc and Humansdorp are surrounded by very open landscapes where the WAsP model is expected to perform well. The sites Risø (Giebel and Gryning2004) and Østerild (Peña2019) are characterized by a mixture of forest and open areas near the mast, which is more challenging in terms of flow modelling. The sites in Sweden and Finland are from confidential projects, and the exact locations of these can therefore not be disclosed.

Table 1Sites and masts used for cross-prediction analyses. The characteristic land cover is indicated. The EPSG code identifies the coordinate reference system, including the map projection, datum, and zone number (see text), and the x, y, and z columns identify the positions in projected map coordinates.

Download Print Version | Download XLSX

3.2 Surface roughness maps

3.2.1 Standard land cover data sets

Four coarse-resolution land cover databases (spatial discretization between 100 and 1000 m) that are regularly applied for wind energy modelling are used in connection with this analysis. These are the Global Land Cover Characterization (GLCC 1000), the MODIS MCD12Q1 V6 product (MODIS 500), the C3S Global Land Cover product (GLOB 300), and for the European sites the Corine Land Cover inventory (CORINE 100). Further details about the land cover products are given in Table 2.

USGS EROS Archive (1993)Friedl and Sulla-Menashe (2019)ESA CCI (2015)2021

Table 2Summary of the different coarse-resolution land cover products (spatial discretization between 100 and 1000 m) used for creating the roughness maps.

Download Print Version | Download XLSX

The coarse-resolution land cover data sets come with standard roughness conversion tables (Appendix A, Tables A1A5). The z0 for forest classes tends to be too low for temperate forests (Floors et al.2018). Recent studies of Dörenkämper et al. (2020) and Badger et al. (2015) have highlighted this shortcoming and have suggested higher z0 for forested areas. To reflect the uncertainty of z0 assignments in forested areas, we include the original vs. revised values in Tables A1A5. In connection with the C3S Global Land Cover product, the z0 classification is based on the data set of 2009, and some of the 23 land cover classes have been split into sub-classes since then. Therefore, z0 of each subclass is assumed to be identical to z0 of the class it inherits from (see Table A3). It is mostly classes with forests that have been split up, and one could get a better estimation of z0 by a more detailed analysis of the canopy structure in these subclasses. However, this approach is not attempted here.

3.2.2 Novel Sentinel data sets

High-resolution satellite-based data packages are developed for an area of 40×40 km around each of the eight sites described above. Each data package has a regular grid spacing of 20×20 m and includes (1) land cover classification, (2) LAI, and (3) forest canopy height. These data are derived from satellite imagery using machine learning methods. Specifically, random forests (Breiman1996) are used to classify land cover and for down-scaling LAI, while support vector machines (Boser et al.1992) are used for canopy height estimation. Both algorithms are from a family of supervised machine learning methods that are routinely applied to satellite imagery to extract features of interest. This is done by supplying the algorithms with a set of samples from the satellite imagery that have known labels (e.g. land cover class or canopy height). These samples are used to train the machine learning models whereby the algorithms learn to predict each feature of interest across the entire satellite image. The primary data source used for the production of the three data layers is Copernicus satellite imagery from the Sentinel-1 and Sentinel-2 missions. Sentinel-1 is a c-band (5.6 cm) synthetic aperture radar sensor, while Sentinel-2 is an optical sensor providing data in the visible, near-infrared, and shortwave infrared parts of the spectrum.

The land cover classification is based on five land cover classes most relevant for wind modelling (Appendix A, Table A5). For each site, training data are collected for each land cover class by manually labelling areas in satellite imagery of known land cover. These training data are used as a dependent variable in a random forest model to predict land cover for each 20×20 m grid cell. Independent variables used for the land cover classification model include observations from the Sentinel-1 and Sentinel-2 sensors. LAI has been estimated for all grid cells identified as forest in the land cover classification through down-scaling of coarse-resolution LAI from the MODIS sensor to the 20×20 m grid of the Sentinel-2 observations (Guzinski and Nieto2019).

Forest canopy heights are estimated for all grid cells identified as forest in the land cover classification. Canopy height estimates from the ICESat-2 sensor are used as a dependent variable to train a support vector regression model to predict forest canopy height for each 20×20 m grid cell. Observations from Sentinel-1 and Sentinel-2 are used as independent variables in the regression model. ICESat-2 provides a globally available data set from the US National Aeronautics and Space Administration (NASA) that uses space-borne lidar technology to estimate ground and vegetation heights. It does this with six laser beams, scanning a swath of terrain 9 km wide, with each beam having a footprint of 17 m diameter. Because of the Sentinel sensors' importance for obtaining the three layers described above, they are labelled “Sentinel” throughout the rest of this paper.

3.2.3 Reference data sets

To validate the novel Sentinel maps, we consider two types of roughness maps for reference: (1) from manual digitization and (2) from lidar scans for sites where such data are readily available from previous works. For an experienced wind resource engineer, the most accurate way to obtain a land cover data set is by visual inspection in the field or using photographs followed by digitization of the land cover areas with the biggest impact on the flow. A representative z0 is then assigned to each of these areas. Maps digitized manually are available for Cuauhtémoc, Humansdorp, Mérida, and Risø (see Table 3). Another approach to obtain an accurate roughness map is by using lidar scans. Tree heights can be retrieved from the cloud point data delivered by airborne lidars. The lidar data sets are described further in Sect. 3.3 as they are also used to map the terrain elevation. Once h is obtained from lidar scans, the procedure described in Sect. 3.2.1 in Floors et al. (2018) is used to estimate z0 and d for Ryningsnäs, Sweden, and Østerild. The reference maps are denoted as “manually digitized or lidar scans” throughout the rest of the paper (see Table 3).

Table 3Overview of the input data that are used for each site to estimate wind resources. The meaning of the different abbreviations are discussed in the text.

Download Print Version | Download XLSX

3.3 Elevation maps

In addition to the land surface roughness, a flow model needs inputs of the terrain elevation. Lidar scans can be used to retrieve the terrain elevation and also to estimate the canopy height at forested sites as described above (Popescu et al.2003; Floors et al.2018). Elevation data from aerial lidar scans are available for five of the eight sites analysed here (Table 3). The Swedish sites are covered by a national laser campaign conducted in 2013 (Lantmäteriet2021). The retrieved elevation data are available with a 20×20 m grid spacing. Lidar scans are available for the Danish site Østerild at the same spatial resolution, whereas the Finland elevation model is a digital terrain model produced by the National Land Survey of Finland (Maanmittauslaitos, MML) and available at a 10×10 m grid. For Risø, elevation data are obtained from the 2.5 m contour lines from “Danmarks Højdemodel” (Styrelsen for Dataforsyning og Effektivisering2016).

Elevation data for the two Mexican sites, Cuauhtémoc and Mérida, and the site in South Africa, Humansdorp, are obtained from the Shuttle Radar Topography Mission (NASA JPL2013) at 90 m resolution. We note that wind engineers should be careful when using databases such as SRTM in combination with a displacement height. Because the SRTM product shows the height of the surface, it includes the tree height. If this is not taken into account, the turbine or mast might be placed at too high elevation, and there is a risk of double-counting of the effects of the forest.

4 Methodology

The WAsP model is a microscale flow model that is frequently employed for wind resource assessments (Troen and Petersen1989). It contains submodels for orography, roughness changes, obstacles, and stability effects. In the following, we explain how the roughness and orographic submodels within WAsP are modified to better utilize the high spatial detail of z0 and d obtained from the Sentinel data set. We then describe the cross-prediction method that is applied to test the novel Sentinel data sets against more conventional input data for wind energy flow modelling.

4.1 Model setup

The WAsP model consists of a graphical user interface and a core model that is written in the programming language Fortran. In the following, we refer to the core model code that is directly accessed via a set of Python codes. The routines are not yet available in the graphical user interface, but a new WAsP version scheduled for release in late 2021 will include them. One of the main advantages of the WAsP core model is its speed: it typically takes seconds to calculate the annual energy production (AEP). The newly implemented routines are parallelized using OpenMP in the Fortran language. Because in the WAsP core each grid point is independently calculated, the problem is easily distributed across central processing units. The standard settings of the WAsP core model (corresponding to WAsP version 12.6) are used here unless otherwise specified. In the following, we describe the modifications to the roughness and orographic submodules. The different parts of the model chain described in this section are shown in Fig. 1.

Figure 1Process diagram of the submodels in WAsP. The rounded rectangles denote data structures, and the square boxes denote a method.


4.1.1 Forest submodel

The ORA and Raupach models are used to obtain z0 and d and were implemented as described in Sect. 2.2. The SCADIS model is run without Coriolis force, because we want to obtain z0 and d for the logarithmic wind profile in the surface layer. The height of the model domain is specified as 4h. The height of the first model level is specified to be z1=0.5z0b, where the background roughness length z0b=0.3 m. The SCADIS model is run from an initial logarithmic profile (z0=z0b) with a wind speed of 10 m s−1 at the domain top. The time step in the model is set to 300 s. The model integration is terminated when the wind speed at the canopy top changes less than 0.001 m s−1. Finally, z0 and d are then found by fitting dUdz=u*κ(z-D) to U and u* obtained in the range from 1.5h to 2.5h.

4.1.2 Roughness submodel

Previous versions of WAsP use an algorithm, which finds all z0 crossings between the roughness lines defined in a map and a number of rays extending from the point of interest. In this study we are not only interested in z0 of a land cover patch, but we also want to take into account the effect of a displacement height, d. This information has to be passed from the land cover contour lines to the roughness model. The roughness submodel in WAsP uses vector lines as input, and therefore, the satellite-based raster maps of land cover are first converted to a vector-map format.

The new routines expect the input of “identifier” lines; i.e. each line contains an identifier (ID), which represents a certain land cover left and right of the line. This ID is then represented in a land cover table, which prescribes the corresponding z0, d, and a description of that land cover class (Appendix A, Tables A1A5). Keeping the land cover table separate from the map has several advantages, like the possibility to perform sensitivity studies with respect to z0 and d, which has been difficult so far because one has to modify the contour map itself in an external program.

Using the approach above, we can process the standard land cover maps (see Sect. 3.2.1). The novel Sentinel maps contain additional layers with h and LAI. h is discretized into bins of 5 m and the LAI into bins of 1. The result is a Sentinel land cover raster map which typically contains 10–40 different classes of forest types plus the classes specified in Table A5.

For the forest class the centre of the h and LAI bins are used to estimate z0 and d according to the three roughness models described in Sect. 2.2. Because the forest submodel operates on a table and not on the contour lines in the map itself, the conversion of h and LAI to z0 and d is fast (i.e. the speed of the computation scales with the number of entries in the land cover table).

In order to introduce land cover, a new routine, here referred to as a “spider grid analysis”, is developed. It uses a polar zooming grid, similarly to the orographic flow model (Troen1990). The advantage of using a zooming grid is that it concentrates the resolution where it is most needed, and we can use arbitrarily distributed points. The latter is for example beneficial for calculating the wind climate at exact positions of wind turbines. The distance to the first radial segment r0 in the zooming grid is defined by the user (default r0=25 m), and each next segment has a grid spacing that is 5 % larger than the previous one. The number of azimuthal bins (i.e. wind direction sectors) can also be specified by the user and by default is set to 12, i.e. using a sector width of 30. The first sector is always centred at the north.

For each cell in the zooming grid, the fraction of the total area fi that each of a total of N land cover types in the land cover table occupies is determined, and the roughness length is calculated,

(10) ln z 0 = i = 1 N f i ln z 0 i .

The displacement height is taken into account similarly to z0 and is calculated for each cell in the zooming grid as

(11) d = i = 1 N f i d i .

Figure 2Example of a zooming grid analysis at the Ryningsnäs site up to 3 km distance showing the displacement height h (a), the roughness length z0 (b), and the corresponding reduced z0 rose (c) after finding the most significant z0 changes.


The zooming-grid analysis of d and z0 from the Sentinel ORA 20 map over Ryningsnäs is shown in Fig. 2a and b, respectively, as circles with a radius of 3 km. For the 0 sector for an area up to 2 km away from the mast, d≈10 and z0≈2 m, while for the 150 sector at around 2 km distance there is a lake that causes lower z0. Another detail visible from close inspection of Fig. 2a is the very low values of d and z0 in the 150 sector at distances less than 100 m away from the mast. This is because of a clearing in the forest in that direction (also seen in other sectors), which has important implications for the flow modelling at Ryningsnäs, as will be further discussed in Sect. 6. We note that Ryningsnäs was extensively investigated by Bergström et al. (2013), who fitted a logarithmic wind profile to the measurements to obtain z0 and d. For the 150 sector, d was close to zero due to the clearing, while z0≈2 m. Their value is presumably higher than the z0 of the clearing because it was determined at 30–40 m and therefore has a different footprint area.

Figure 2c shows a so-called reduced roughness rose, which shows the most significant z0 changes in all directions for the same area as the zooming grid analysis. This roughness rose captures the main features of the area around the site, but it also misses some features that perhaps would have been identified using the most accurate and detailed method of manual digitalization. The large number of relatively small clearings and forest patches in all directions clearly makes it challenging to find the most significant z0 changes.

To account for roughness changes, WAsP calculates sector-wise speedup factors for a certain point. Because the effect of a roughness change on the wind speed in a sector is distance dependent, with nearby areas having a higher impact, the distance to each z0 area is multiplied with an exponential weighting function as described in Floors et al. (2018). From these weighted values z0w, the ones that explain most of the variance of z0w are stored for further processing (up to a maximum of nmax). This is done for computational efficiency, so that equations that take into account the effect of internal boundary layers can be used. These equations are given in Sect. 8.3 of Troen and Petersen (1989). The output of these equations are sector-wise speed-up factors, which are used to remove the effect that microscale features may have on the wind observations. Apart from the speed-up factors, z0w is also used to compute a geostrophic (sometimes referred to as an effective or mesoscale) roughness length z0G.

Similarly to z0, d also has to be filtered to use in Eq. (1) and the geostrophic drag law. Therefore, a triangular weighting function is applied to the zooming d-grid. The average dG in each sector is found by taking the triangular weighted average up to a distance xd=10d0, with d0 defined as d at the first cell in the zooming grid at a distance x0 from the origin. The triangular weight w is 1 at x0 and w=0 when x>xd. Physically the reason for applying this filter is that it takes some time for a new logarithmic profile to develop. Thus, taller trees will need a fetch to lead to an effective dG that is applied in Eq. (1). The sector-wise z0G and dG can finally be used in connection with Eq. (1) and the geostrophic drag law to find a geostrophic wind climate from an input histogram. This procedure is unaltered from the description in chap. 8.7 of Troen and Petersen (1989) and is therefore not further described.

4.1.3 Orographic submodel

The submodule for orography computes speed-ups due to elevation, and WAsP uses the Bessel expansion on a zooming grid (BZ) model (Troen1990). The input to the orographic model is a map with elevation lines, which are processed into a polar zooming grid, with the highest radial resolution at the centre point of the grid. The resolution of the zooming grid depends on the radius R being large enough so that the entire map of height contours is contained inside the circle with radius R. Contour lines that are more than 20 km away from the site are ignored.

Here we are interested in the effects of z0 on flow modelling, and therefore the highest-quality terrain elevation map is chosen for each of the eight sites (see Sect. 3.3). We then study the impact of varying the land cover maps only, while keeping the elevation map constant.

The resulting zooming grid of d obtained by using Eq. (11) is added to the terrain elevation zooming grid. Using this grid, the methodology described in Troen (1990) is used to calculate sector-wise wind-speed-independent orographic speed-ups. Similarly to the roughness speed-ups, these are used as local perturbations to the flow, which can be used to obtain a wind climate that is representative for a large area.

4.2 Cross-prediction analyses

We first perform cross-prediction analyses using the standard land cover databases (GLCC 1000, MODIS 500, GLOB 300, and CORINE 100) in combination with their respective un-modified translation tables leading to z0 (Appendix A). Then, we consider the same input data sets with the modified tables from Dörenkämper et al. (2020) and Badger et al. (2015). Finally, we perform the cross-predictions using the novel Sentinel data sets, the manually digitized maps, and the lidar scans (see Table 3).

At some sites, there are multiple masts measuring at the same time (see Table 1), and the wind climate from one mast is used to predict the wind resource at the neighbouring mast. At other sites, only one mast is available and only vertical cross-predictions are possible, i.e. predictions of the wind distribution from one height to another at the same mast. We exclude observations within the roughness sublayer where the WAsP model does not apply, so only measurements between 20 and 200 m are used. In total there are 914 cross-predictions for the CORINE 100 (only Europe) database and 950 for all the other databases and Sentinel-based maps. This number is based on all possible cup anemometer pairs per site in Table 1, subtracting the combinations where a cup anemometer input is used to predict the wind climate at the exact same point (a so-called self-prediction).

Observed wind climates are generated for each height and mast from the 10 min time series of wind speed U and wind direction θ, which are discretized into histograms using a bin width of 1 m s−1 and 30, respectively.

The power per unit area swept by wind turbine blades, P (also referred to as the power density), is given by

(12) P = 0.5 ρ U 3 ,

where ρ is the air density. In WAsP, the wind distribution is described by sector-wise Weibull distributions, and therefore we can conveniently find P from the Weibull parameters A and k for Nd sectors as

(13) P = ρ i = 1 N d ϕ i A i 3 Γ 1 + 3 / k i ,

where ϕ is the sector frequency, Γ is the gamma function, and ρ is calculated according to the methods described in Floors and Nielsen (2019). Throughout this paper we choose to evaluate relative errors of power density,

(14) ε P = 100 P mod P obs - 1 ,

where P is obtained from the modelled (mod) or observed (obs) Weibull distributions. The statistics reported throughout the rest of the paper are the bias, εP, and the root-mean-square error (RMSE), εP2, where the overbar denotes a mean. We also include relative errors in wind speed U, which we obtain as in Eq. (14) but instead of P, U is computed from the Weibull distributions as

(15) U = i = 1 N d ϕ i A i Γ 1 + 1 / k i .
5 Results

In the following, we show results of the three different forest models that were implemented. We present the outcome of cross-prediction analyses per site, and subsequently, we aggregate results for all eight sites investigated.

5.1 Model response to LAI

The behaviour of the three forest roughness models for a canopy height h of 10 m as a function of LAI is shown in Fig. 3. The ORA model does not depend on LAI and is therefore constant, i.e. d/h=2/3 and z0/h=0.1. The Raupach and SCADIS models both show an increasing d as the LAI increases. z0 has the opposite behaviour and decreases to lower values after an initial maximum of around LAI  1.

Figure 3Response of the different roughness length and displacement height models to leaf area index (LAI) for a forest with h=10 m. The SCADIS model is set up with α=9 and β=3, which are typical values for forests with most of the canopy density in the upper part of the canopy layer (Sogachev et al.2017).


The main differences between the SCADIS and Raupach models occur for a relatively low LAI ∼1, using the SCADIS model z0≈2 m while using the Raupach model z0≈1 m. The same holds for d, where using the Raupach model yields a d that is nearly twice that of the SCADIS model. For the more commonly occurring LAI >3, the differences are minor. When a different canopy profile is specified, i.e. using α=5 and β=3, there are larger differences between the two models (not shown).

5.2 Roughness length maps

Figure 4 illustrates the different z0 representations for a 6×6 km area around the Ryningsnäs mast. The GLCC database has the coarsest resolution of 1000 m, and therefore, it fails to capture any detail near the mast (Fig. 4a). The whole area is represented by land cover class 14, i.e. evergreen needle-leaf forest, which corresponds to z0=1.5 m (Table A1). Likewise, the MODIS database at 500 m resolution does not capture any detail within the selected area (Fig. 4b), and everything is represented as evergreen needle-leaf forest with z0=1.0 m (Table A2). GlobCover, with a resolution of 300 m, represents most of the land cover around the site as evergreen needle-leaf forest with z0=1.5 m (Table A1) and can capture some of the lakes and open areas (Fig. 4c). The Corine land cover database captures more details due to its higher resolution of 100 m (Fig. 4d). Most of the area around the mast is classified as “mixed forest” with z0=1.1 m (Table A4). The lakes to the southeast of the masts and the open area to the west of the masts are captured. The Sentinel ORA 20 map (Fig. 4e) shows the result of combining the land cover map with the tree height and LAI layers according to the procedure presented in Sect. 3.2.2 and using the conversion of h to z0 and d from ORA (Sect. 2.2.1). The Sentinel SCADIS and Raupach maps look very similar to this and are therefore not shown. It can be seen that the map generally contains a wider range of z0, ranging from 0.0002 to over 3 m.

Figure 4Roughness lengths obtained from a 6×6 km square around the Ryningsnäs site (red point) from standard land cover data sets (a–d), novel Sentinel data sets (e), and aerial lidar scans (f). The time periods for which the different data sets are obtained are shown in Table 2.

Figure 5Roughness lengths obtained from a 6×6 km square around the meteorological mast (red point) from Cuauhtémoc (a), Humansdorp (b), Mérida (c), Risø (d), Ryningsnäs (e), and mast 1 at Østerild (f) from the Sentinel ORA 20 maps.

Finally, the z0 map obtained using the ORA approach with h obtained from lidar scans is shown in Fig. 4f. Spatially comparing the Sentinel-based tree heights (Fig. 4e) with tree heights derived from aerial lidar scans (Fig. 4f) reveals that the areas with lower and higher h correspond very well. However, z0 derived from Sentinel is generally lower than from the lidar scans. In Fig. 5 we show the roughness length around all the sites except the Finnish and Swedish sites, which are confidential.

Figure 6Maps showing the terrain elevation (a), the displacement height (b), the power density (c), and the relative difference in power density compared to a terrain map where d=0 (d) at Ryningsnäs for westerly winds.

As an example of the new routines presented in Sect. 4.1.2, we use the observed wind speed histogram at the Ryningsnäs mast at 80 m and predict the wind resource at 80 m for every point in the 6 × 6 km area shown in Fig. 4 using a raster with a resolution of 40 m. Figure 6 shows the resulting maps of the terrain elevation, displacement height, power density, and relative difference before and after d is taken into account.

The elevation changes are modest at the Ryningsnäs site (Fig. 6a), but there is a hill to the west. In Fig. 6b, the displacement height in each point is shown. For clarity, we only focus on d for the most frequently observed westerly sector. The clearing located southeast of the mast (which is in the centre of the map) is visible, resulting in lower displacement heights. In addition, displacement heights of more than 10 m are visible further from the mast. The result of combined elevation, displacement, and roughness description on the emergent power density in the area is shown in Fig. 6c. The higher power density due to the flow speed-ups over the hill in the west is clearly visible. Because the changes in power density due to the introduction of a displacement height are difficult to discern, Fig. 6d shows the difference in percent between P with and without using a displacement height. Using d causes a decrease in power density between 0 % and 12 %. The area with the highest reduction in power density in the centre indeed corresponds closely to the area with the largest concentration of high displacement heights in Fig. 6b.

5.3 Cross-predictions at wind energy sites

The rms of εP per site is shown in Fig. 7. We have grouped the maps from the lidar scans and those that were manually digitized into one class, because we consider all these as reference maps, i.e. the best possible we can achieve. We find rather large errors for Mérida, Ryningsnäs, and Finland. All are forested sites with mixed vegetation patterns. This indicates a need for better flow modelling at such inhomogeneous sites. Humansdorp and Cuauhtémoc, in contrast, are characterized by a relatively simple land cover, which results in a low rms of εP. The cross-predictions performed using the novel Sentinel maps as input mostly lead to lower rms of εP than those using the standard land cover databases. However, at Mérida and Østerild the rms of εP from the Sentinel-based maps is not lower than those from the standard land cover databases. The former is a very complex site that is characterized by forest in all directions. One possible explanation for the higher errors at this site can be a clearing in the forest that is shown to the southeast in the Sentinel maps, but which does not appear in reality. At Østerild the dominating wind direction is from the west, which is subject to many wind breaks. These are not generally detected in the Sentinel-based roughness maps but contribute significantly to a higher z0. Therefore, the background roughness for grassland that was used in the Sentinel maps (see Table A5) might be too low at Østerild. The three forest roughness models lead to very similar results for all sites.

Figure 7The RMSE in power density (εP, %) obtained when each available wind sensor is used to predict the power density at the location of other wind sensors located at the same site. Values are shown for all sites. Different colours refer to different roughness length maps used as input to the calculation.


If we consider only the cross-predictions for sites where aerial lidar scans are available (i.e. Ryningsnäs, Sweden, Finland, and Østerild), we can calculate that the combined rms of εP≈10.9 % for the aerial lidar scans (lidar scans ORA 20) and 10.7 % for the Sentinel ORA 20 maps, respectively. At Ryningsnäs, Finland, and Østerild the aerial lidar scans yield lower rms of εP than the Sentinel maps, whereas at the Swedish site results are more comparable.

When we consider the sites with manually digitized maps only (i.e. Humansdorp, Mérida, Cuauhtémoc, Risø), the rms of εP is 10.9 % for the manually digitized map, 13.1 % for the Sentinel ORA 20 map, and 17.3 % for the GLOB 300 map. Thus, averaged over these four sites, satellite-derived estimates of z0 do not yield better power predictions than those based on manually digitized maps.

5.4 Aggregated results of cross-predictions

We aggregate results from the cross-predictions for the eight sites to obtain the average performance of each set of input data in connection with flow modelling in WAsP.

Table 4RMSE in power density (εP, %) and wind speed (εU, %) for all sites combined. For GLOB 300 and CORINE 100 both the original and revised land cover databases are shown (see Sect. 3.2 and Appendix A, Tables A4 and A3).

Download Print Version | Download XLSX

The rms of εP for the cross-predictions using standard land cover databases with original and revised roughness translation tables is shown in Table 4. The GLCC 1000 generally leads to the highest errors, when both the original and revised land cover tables are used. The lowest rms of εP of all is achieved with the GLOB 300 maps in combination with the revised land cover table. The rms of εP is reduced for all land cover databases when the revised land cover tables are used instead of the original ones.

Figure 8RMSE in power density (εP, %) estimated with different input maps for all eight sites combined.


We can now compare results generated with the standard land cover databases with results generated with the novel Sentinel data layers as input. Figure 8 shows the prediction errors when the WAsP model is run with seven different inputs: the three global standard land cover databases with revised roughness translation tables plus three types of Sentinel maps based on the three different forest roughness models described in Sect. 2.2. The CORINE 100 database is not available at all sites and is therefore not shown. The reference map for each site is chosen based on availability of lidar scans or a manually digitized map (see Table 3).

Maps generated with an advanced forest model (Sentinel ORA 20, Sentinel Raupach 20, and Sentinel SCADIS 20) yield lower rms of εP (≈11 %) than maps based on the standard land cover databases. Table 4 shows that also the mean bias of εP is lower in the Sentinel maps than in the standard land cover databases GLCC 1000, MODIS 500, and GLOB 300. Similarly for wind speed, the rms values of εU are smaller for the model runs using Sentinel maps. The mean bias in εU becomes slightly negative but is generally close to zero for all WAsP model runs.

When running WAsP it is often more instructive to evaluate horizontal cross-predictions, because the roughness rose will then be different for two points. To investigate this, our data were filtered to include only horizontal cross-predictions (i.e. sites with more than one mast). However, this did not change the large differences in εP and εU between the model simulations based on standard land cover databases and Sentinel maps. In addition, the effect of the measurement height was investigated by only including cross-predictions between 50 and 200 m, but again, this did not change the general picture that emerges from Table 4.

6 Discussion

For the first time, forest parameters retrieved from the Sentinel and ICESat-2 satellites have been used for wind energy flow modelling. Our results demonstrate that the spatial variability of the land cover within forests (e.g. vegetation with different height and density, clearings, lakes) is resolved far better by these novel products than it is by standard land cover products with pan-European or global coverage. The high level of spatial detail in the novel Sentinel data layers is almost comparable to products derived from aerial lidar scans. This is promising in connection with wind energy flow modelling because these data layers can be produced for any site in the world and updated frequently. The cost of retrieving the satellite-based forest data layers is far lower than the price of dedicated airborne lidar campaigns thanks to the open access to satellite data archives by Copernicus and NASA. From Fig. 7 it is clear that using the Sentinel maps of tree height rather than standard land cover maps has the highest benefit at sites where masts or turbines are surrounded by forest (Ryningsnäs, Finland, and Sweden). Taking the displacement height into account for such sites leads to significantly lower εP. Because there are seven masts at the Swedish site, leading to a large number of cross-predictions, this site has a large impact on the aggregated results (Table 4).

Using roughness derived from Sentinel canopy height maps had the greatest impact at the Ryningsnäs site: reducing the rms of εP by more than 50 % compared to using roughness derived from land cover maps (Fig. 7). This improvement of εP is due to the mast location very close to the forest edge in the westerly sector (i.e. winds coming from 265–285). Because of the large impact of sector-wise displacement heights on the wind profile, the results are highly dependent on specifying the correct h. This indicates that for masts located very close to forest edges, sensitivity studies are recommended, and the Sentinel maps should be carefully calibrated.

The observed difference in h in the novel Sentinel product compared to the lidar scans at Ryningsnäs may partly be due to the uncertainty in sensing and retrieval methods that are required to convert satellite observations to canopy heights and partly due to temporal differences between the satellite and lidar products. The standard canopy height product from ICESat-2, which we used to train our model for forest height estimation, contains tree heights over a 17×100 m along-track transect and this is different from the spatial resolution of the lidar scans. Additionally, tree growth in temperate forests can be up to a few metres per year. Therefore, it is important to use the most updated h parameterization when modelling wind resources. Further research is needed to fully understand and improve the absolute accuracy of tree heights retrieved from satellite observations.

Likewise, the absolute accuracy of forest densities expressed through the LAI should be thoroughly quantified in future work. For example, the LAI used for our analyses was downscaled from a very coarse-resolution product, and the consequence of this downscaling is unknown. Differences between outputs of the three forest canopy models are small in this study. It is, however, possible that the SCADIS model is too advanced for the type of analysis performed here given the large uncertainty of the LAI input and the absence of canopy density profiles. More detailed studies at a higher number of sites, preferably with observed h and LAD profiles, are needed to validate the SCADIS model further.

Our analyses show only minor differences in rms of εP when cross-predictions are performed using Sentinel vs. aerial lidar scans to estimate d and z0. These findings are promising in light of the lower cost and the global coverage of the satellite-based data layers. Manual assessments of the land surface roughness do not necessarily lead to better results. For most sites (three out of five) automated assessments based on the novel Sentinel data sets led to lower error values than manual assessments. In addition, automated procedures speed up the wind resource assessment process and remove the subjective judgement of a siting engineer. Using high-resolution input data for the automated retrieval of z0 is essential. We have demonstrated how the small-scale roughness changes within forests are poorly resolved by standard land cover data sets, and as a consequence, z0 and predictions of the wind power density come with large uncertainties even after revision of the roughness translation tables for forest classes.

Our results indicate that all methods considered in this paper underestimate the roughness length over forest (i.e. higher z0 values lead to better predictions). This is possibly due to roughness length aggregation, i.e. finding z0 for a larger area that correctly prescribes the momentum flux in surface-layer similarity in connection with Eq. (1). In this paper, the larger-scale sector-wise roughness length z0G is obtained via simple logarithmic averaging of the cells in the spider-grid rose applied in Eq. (10), but this is known to be flawed in heterogeneous conditions. For example, Taylor (1987) suggested accounting for sub-grid variance of z0 when calculating an effective roughness length. Vihma and Savijärvi (1991) compared model approaches for several landscape configurations in Finland and found that z0eff was always higher than a simple logarithmic average. This was also confirmed by Hasager and Jensen (1999), who found that the difference between the logarithmic average and z0eff was higher for landscapes with small patches and with a half-to-half mix of rough and smooth patches. For small patches, one can argue that the form drag of the forest will always lead to a much higher z0eff than implied through logarithmic averaging. However, an analysis of the effect of z0 aggregation was not attempted in this study, but its impacts are likely large for heterogeneous sites like Østerild, Finland, and Ryningsnäs. Bottema et al. (1998) reviewed a range of z0 aggregation methods, which could be investigated in future work.

In this study, cross-prediction is used without any manual interference. For an actual wind energy project, much smaller errors would be achieved using the WAsP model because the site engineer usually has access to measurements at several heights, which can be used to fit the model nearly perfectly to the observed mean wind profile. However, such a fitting process can be deceiving, because a good fit is influenced by both z0 and atmospheric stability conditions. It may thus me subject to compensating errors. The effect of atmospheric stability is beyond the scope of this study, and the default settings in WAsP are used, but it is acknowledged that there can be deviations at some sites.

7 Conclusions

We have tested a novel satellite-based product for land surface parameterization in forested areas and quantified the effect of using this product for wind resource modelling. The novel satellite-based product is based on observations from Sentinel-1, Sentinel-2, and ICESat-2, and it contains collocated layers of land cover, tree heights, and LAI at a 20 m spatial resolution. These maps are converted to maps of roughness lengths and displacement heights using three different forest modules of varying complexity. The simplest way to convert the tree height to z0 and d is by multiplying with a constant (e.g. 1/10). Secondly, a physically more advanced approach, which also considers the effect of LAI, is implemented. The third module is a 1D version of the SCADIS model, and it takes the effects of varying LAI and canopy density profiles into account. All three forest modules are used in the WAsP model, which is frequently used for predicting the wind resources from wind measurements.

We show that the high land cover variability of forested landscapes is poorly resolved by global and pan-European land cover databases such as GLCC, MODIS, CORINE, and GlobCover. The novel satellite-based product leads to more detailed maps of z0 and d, which are spatially comparable to aerial lidar scans or manually digitized maps. Retrieving canopy heights from satellite data alone is a rapidly evolving branch of research – driven by the recent release of global calibration data sets such as ICESat-2 and GEDI. Further research is ongoing to refine the satellite-based products. This can be achieved through collection of larger and more accurate training data sets for the tree height retrieval, the integration of ground-based tree height observations, improved geolocation of satellite-derived products, and advancements in machine learning algorithms.

Cross-predictions are performed at eight sites with tall masts to evaluate the effect of using different input data sets in connection with flow modelling in WAsP. The novel maps from Sentinel lead to a reduction of the rms of relative errors in power density at most sites and on average by ≈3 % compared to the best performing roughness map obtained from a coarse-resolution land cover database. This is even after the roughness lengths for specific land cover categories in the coarse-resolution products are improved. Differences between the three forest modules are minor, showing that the sensitivity of the WAsP model to different approaches to obtain z0 and d is low.

The rms of relative errors in power density found for the Sentinel ORA maps (10.7 %) is comparable to those obtained from aerial lidar scans (10.9 %). This finding is very promising because the novel satellite-based maps of z0 and d can be generated at a lower cost and a higher temporal resolution than aerial lidar campaigns. Processing of the satellite-based maps is fully automated. For sites that show a potential for wind power projects, the new routines and products could replace current practices of land cover analysis, which is time-consuming and plagued by subjective assessments.

Appendix A

The land cover tables for the land cover databases that were investigated in this paper are supplied in Tables A1A5. The references to the land cover databases can be found in Sect. 3.2.1.

Table A1Land cover table used for the GLCC database (Thøgersen2021).

Download Print Version | Download XLSX

Table A2Land cover table used for the MODIS-based maps (Thøgersen2021).

Download Print Version | Download XLSX

Table A3Land cover table used for the GlobCover/ESA-CCI database (Thøgersen2021).

Download Print Version | Download XLSX

Table A4Land cover table used for the CORINE database (Thøgersen2021).

Download Print Version | Download XLSX

Table A5Land cover table used for the Sentinel-based maps.

Download Print Version | Download XLSX

Code availability

The numerical results are generated with proprietary software.

Data availability

Sample data packages containing the Sentinel forest data layers for selected sites are available at (EMD international2021). The wind measurements from the mast in South Africa (WM08) are available at (CSIR2021). Data from the two Mexican masts are available at (INEEL2021, registration required). GLCC data can be downloaded at (USGS EROS Archive1993)​​​​​​​. MODIS data can be downloaded at (Friedl and Sulla-Menashe2019). ESA CCI data can be downloaded at (ESA2021). CORINE land cover data can be downloaded at (Copernicus Land Monitoring Service2021).

Author contributions

RF generated the results, contributed to the WAsP model code, and drafted the manuscript. MB coordinated the work and revised the manuscript. IT developed most of the methods and code in the WAsP model and contributed to the preparation of the manuscript. KG created the Sentinel data layers and wrote a section of the manuscript. FHP helped with data processing and preparation of the manuscript.

Competing interests

The high-resolution forest data layers generated from Sentinel are sold by DHI GRAS A/S, where Kenneth Grogan is employed. The WAsP software is maintained and sold by DTU Wind Energy, where Rogier Floors, Ib Troen, and Merete Badger are employed​​​​​​​.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work has received funding from the H2020 e-shape project (grant agreement 820852), from Innovation Fund Denmark through the InnoWind project (6172-00004B), and from the Ministry of Foreign Affairs of Denmark administered by Danida Fellowship Centre through the “multi-scale and model-chain evaluation of wind atlases” (MEWA) project (17-M01-DTU). We acknowledge the providers of global and pan-European land cover data sets​​​​​​​: US Geological Survey for GLCC and MODIS, ESA Climate Change Initiative and in particular its Land Cover project as the source of the CCI-LC database, and ESA Land Monitoring Service for CORINE. ESA and the European Commission are acknowledged for Copernicus Sentinel data. NASA are acknowledged for ICESat-2 data, and in particular for regular consultation through the Early Adopter Program. Ebba Dellwik from DTU is thanked for preparation of aerial lidar data and Andreas Bechmann for his review of the paper. We acknowledge the Wind Atlas for South Africa project for the wind observations near Humansdorp (CSIR2021) and the Mexican Wind Atlas project for observations at the Mexican sites (INEEL2021). Technicians of DTU are acknowledged for the maintenance of the Risø and Østerild masts.

Financial support

This research has been supported by the Innovationsfonden (grant no. 6172-00004B), the Horizon 2020 Framework Programme, H2020 Societal Challenges (e-shape (grant no. 820852)), and the Danida Fellowship Centre (grant no. 17-M01-DTU).

Review statement

This paper was edited by Rebecca Barthelmie and reviewed by three anonymous referees.


Badger, J., Hahmann, A., Larsen, X. G., Badger, M., Kelly, M., Davis, N., Olsen, B. T., and Mortensen, N. G.: The Global Wind Atlas, Tech. rep., DTU Wind Energy, Roskilde, Denmark, available at:​​​​​​​ (last access: 26 October 2021), 2015. a, b

Bergström, H., Alfredsson, H., Arnqvist, J., Carlén, I., Dellwik, E., Fransson, J., Ganander, H., Mohr, M., Segalini, A., Söderberg, S., Bergström, H., Alfredsson, H., Carlén, J., Dellwik, I., Ganander, J., and Mohr, H.: Wind power in forests: Winds and effects on loads, Tech. rep., Uppsala University, Stockholm, Sweden, available at:​​​​​​​ (last access: 26 October 2021), 2013. a

Blackadar, A. K. and Tennekes, H.: Asymptotic Similarity in Neutral Barotropic Planetary Boundary Layers, J. Atmos. Sci., 25, 1015–1020,<1015:ASINBP>2.0.CO;2, 1968. a

Boser, B. E., Guyon, I. M., and Vapnik, V. N.: A Training Algorithm for Optimal Margin Classifiers, in: Proceedings of the 5th Annual ACM Workshop on Computational Learning Theory, 27–29 July​​​​​​​ 1992, Pittsburg, USA, ACM Press, 144–152, 1992. a

Bottema, M., Klaasen, W., and Hopwood, W.: Landscape Roughness Parameters for Sherwood Forest – Validation of Aggregation Models, Bound.-Lay. Meteorol., 89, 317–347,, 1998. a

Breiman, L.: Bagging predictors, Mach. Learn., 24, 123–140, 1996. a

Businger, J. A., Wyngaard, J. C., Izumi, Y., and Bradley, E. F.: Flux-Profile Relationships in the Atmospheric Surface Layer, J. Atmos. Sci., 28, 181–189,<0181:FPRITA>2.0.CO;2, 1971. a

Chen, J. and Black, T.: Defining Leaf-Area Index for non-flat leaves, Plant Cell Environ., 15, 421–429,, 1992. a

Copernicus Land Monitoring Service: CORINE Land Cover 2018, Copernicus Land Monitoring Service [data set], available at:, last access: 26 October 2021 a, b

Council of Scientific & Industrial Research (CSIR): WASA data, available at:, last access: 26 October 2021. a, b

Csillik, O., Kumar, P., and Asner, G. P.: Challenges in estimating tropical forest canopy height from planet dove imagery, Remote Sens., 12, 1160,, 2020. a

De Bruin, H. A. and Moore, C. J.: Zero-plane displacement and roughness length for tall vegetation, derived from a simple mass conservation hypothesis, Bound.-Lay. Meteorol., 31, 39–49,, 1985. a

Dörenkämper, M., Olsen, B. T., Witha, B., Hahmann, A. N., Davis, N. N., Barcons, J., Ezber, Y., García-Bustamante, E., González-Rouco, J. F., Navarro, J., Sastre-Marugán, M., Sīle, T., Trei, W., Žagar, M., Badger, J., Gottschall, J., Sanz Rodrigo, J., and Mann, J.: The Making of the New European Wind Atlas – Part 2: Production and evaluation, Geosci. Model Dev., 13, 5079–5102,, 2020. a, b

EMD international: INNOWIND data layers, available at:, last access: 26 October 2021. a

Enevoldsen, P.: Onshore wind energy in Northern European forests: Reviewing the risks, Renewable and Sustainable Energy Reviews, 60, 1251–1262,, 2016. a

Enevoldsen, P.: Managing the Risks of Wind Farms in Forested Areas: Design Principles for Northern Europe, PhD thesis, Aarhus university, Aarhus, Denmark, 2017. a, b

ESA: ESA-CCI Land Cover, ESA [data set], available at:, last access: last access: 26 October 2021. a

European Space Agency (ESA) Climate Change Initiative (CCI): Land cover classification gridded maps from 1992 to present derived from satellite observations, v2.0.7, available at:!/dataset/satellite-land-cover?tab=overview​​​​​​​ (last access: 26 October 2021), 2015. a

Fagua, J. C., Jantz, P., Rodriguez-Buritica, S., Duncanson, L., and Goetz, S. J.: Integrating LiDAR, multispectral and SAR data to estimate and map canopy height in tropical forests, Remote Sens., 11, 2697,, 2019. a

Floors, R. and Nielsen, M.: Estimating Air Density Using Observations and Re-Analysis Outputs for Wind Energy Purposes, Energies, 12, 2038,, 2019. a

Floors, R., Enevoldsen, P., Davis, N., Arnqvist, J., and Dellwik, E.: From lidar scans to roughness maps for wind resource modelling in forested areas, Wind Energ. Sci., 3, 353–370,, 2018. a, b, c, d, e, f, g, h, i

Friedl, M. and Sulla-Menashe, D.: MCD12Q1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V006, distributed by NASA EOSDIS Land Processes DAAC, NASA [data set],, 2019. a, b

Giebel, G. and Gryning, S.-E.: Shear and stability in high met masts, and how WAsP treats it, available at: and-stability-in-high-met-masts-,-and-how-it-Giebel-Gryning/4624ab387a8135a4437aae4dd1df1276b3b4f302#citing-papers​​​​​​​ (last access: 26 October 2021), 2004. a

Global Wind Energy Council: Global Wind Energy Report: Annual Market Update 2019, available at:​​​​​​​ (last access: 26 October 2021), 2019. a

Guzinski, R. and Nieto, H.: Evaluating the feasibility of using Sentinel-2 and Sentinel-3 satellites for high-resolution evapotranspiration estimations, Remote Sens. Environ., 221, 157–172,, 2019. a, b

Hasager, C. B. and Jensen, N. O.: Surface-flux aggregation in heterogeneous terrain, Q. J. Roy. Meteor. Soc., 125, 2075–2102,, 1999. a

Huang, H., Liu, C., and Wang, X.: Constructing a finer-resolution Forest Height in China Using ICESat/GLAS, Landsat and ALOS PALSAR data and height patterns of natural forests and plantations, Remote Sens., 11, 1740,, 2019. a

Jancewicz, K. and Szymanowski, M.: The Relevance of Surface Roughness Data Qualities in Diagnostic Modeling of Wind Velocity in Complex Terrain: A Case Study from the Śnieżnik Massif (SW Poland), Pure Appl. Geophys., 174, 569–594,, 2017. a

Kelly, M. and Jørgensen, H. E.: Statistical characterization of roughness uncertainty and impact on wind resource estimation, Wind Energ. Sci., 2, 189–209,, 2017. a, b

Lantmäteriet:, last access: 18 March 2021. a

Li, W., Niu, Z., Shang, R., Qin, Y., Wang, L., and Chen, H.: High-resolution mapping of forest canopy height using machine learning by coupling ICESat-2 LiDAR with Sentinel-1, Sentinel-2 and Landsat-8 data, Int. J. Appl. Earth Obs., 92, 102163,, 2020. a, b

Meyers, T. and Tha Paw U, K.: Testing of a higher-order closure model for modeling airflow within and above plant canopies, Bound.-Lay. Meteorol., 37, 297–311,, 1986. a

NASA JPL: NASA Shuttle Radar Topography Mission Global 3 arc second number, distributed by NASA EOSDIS Land Processes DAAC, NASA [data set],, 2013. a

National Institute of Electricity and Clean Energies (INEEL): MEWA data, available at:, last access: 26 October 2021. a, b

Neuenschwander, A. and Pitts, K.: The ATL08 land and vegetation product for the ICESat-2 Mission, Remote Sens. Environ., 221, 247–259,, 2019. a

Peña, A.: Østerild: A natural laboratory for atmospheric turbulence, J. Renew. Sustain. Energy, 11, 063302,, 2019. a

Popescu, S. C., Wynne, R. H., and Nelson, R. F.: Estimating plot-level tree heights with lidar: Local filtering with a canopy-height based variable window size, Comput. Electron. Agr., 37, 71–95,, 2003. a

Raupach, M. R.: Drag and drag partition on rough surfaces, Bound.-Lay. Meteorol., 60, 375–395,, 1992. a

Raupach, M. R.: Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index, Bound.-Lay. Meteorol., 71, 211–216,, 1994. a, b, c, d

Shaw, R. H. and Pereira, A. R.: Aerodynamic roughness of a plant canopy: A numerical experiment, Agric. Meteorol., 26, 51–65,, 1982. a

Sogachev, A. and Panferov, O.: Modification of two-equation models to account for plant drag, Bound.-Lay. Meteorol., 121, 229–266,, 2006. a

Sogachev, A., Menzhulin, G. V., Heimann, M., and Lloyd, J.: A simple three-dimensional canopy – planetary boundary layer simulation model for scalar concentrations and fluxes, Tellus B, 54, 784–819,, 2002.  a

Sogachev, A., Cavar, D., Kelly, M. C., and Bechmann, A.: Effective roughness and displacement height over forested areas, via reduced-dimension CFD, Tech. rep., DTU, Roskilde, Denmark, 2017. a, b

Styrelsen for Dataforsyning og Effektivisering: Danmarks Højdemodel, DHM/Terræn, Tech. Rep. August, Styrelsen for Dataforsyning og Effektivisering, available at: (last access: 26 October 2021), 2016. a

Taylor, P. A.: Comments and further analysis on effective roughness lengths for use in numerical three-dimensional models, Bound.-Lay. Meteorol., 39, 403–418,, 1987. a

Thøgersen, M.: EMD Wiki, available at: (last access: 26 October 2021), 2021. a, b, c, d

Thom, A. S.: Momentum absorption by vegetation, Q. J. Roy. Meteorol. Soc., 97, 414–428,, 1971. a, b

Troen, I.: A high resolution spectral model for flow in complex terrain, American Meteorological Society, Boston, MA, USA, 417–420, 1990. a, b, c

Troen, I. and Petersen, E. L.: European Wind Atlas, Risø National Laboratory, Roskilde, Denmark, 1989. a, b, c, d, e, f

USGS EROS Archive: Land Cover Products – Global Land Cover Characterization (GLCC), USGS [data set],, 1993. a, b

Vihma, T. and Savijärvi, H.: On the effective roughness length for heterogeneous terrain, Q. J. Roy. Meteor. Soc., 117, 399–407,, 1991. a

Short summary
Wind turbines are frequently placed in forests. We investigate the potential of using satellites to characterize the land surface for wind flow modelling. Maps of forest properties are generated from satellite data and converted to flow modelling maps. Validation is carried out at 10 sites. Using the novel satellite-based maps leads to lower errors of the power density than land cover databases, which demonstrates the value of using satellite-based land cover maps for flow modelling.