Numerical modelling of the wind over forests: roughness versus canopy drag
Dalibor Cavar
Mark Kelly
Ebba Dellwik
Tobias Klaas
Paul Kühn
Parameterizing the effect of verticallydistributed vegetation through an effective roughness (z_{0,eff}) – whereby momentum loss through a threedimensional foliage volume is represented as momentum loss over an area at one vertical level – can facilitate the use of forest data in flow models, to any level of detail, and simultaneously reduce computational cost. Results of numerical experiments and comparison with observations show that a modelling approach based on z_{0,eff} can estimate wind speed and turbulence levels over forested areas, at heights of interest for wind energy applications (∼60 m and higher), but only above flat terrain. Caution must be exercised in the application of such a model to zones of forest edges. Advanced flow models capable of incorporating local (distributed) drag forces are recommended for complex terrain covered by forest.
Comprehensive investigation is needed to estimate the wind energy potential at any given site. Along with field measurement campaigns, different kinds of models are widely used nowadays. For sites of moderate or high complexity, RANS (ReynoldsAveraged NavierStokes) solvers have become a commonlyused tool in the wind energy industry (e.g. Sørensen et al., 2017). Increasing detail in the description of the site in the RANS model, along with improving the quality of the model, tends to increase the computing cost. Trees present at the site add a level of complexity to wind flow modelling. Due to rapid technological improvement in airborne lidar scanning (ALS), data on land surface and vegetation characteristics have become available, with increasing detail and quality. Methods have been developed to allow calculation of highly heterogeneous and realistic Plant Area Density (PAD) using raw ALS data (Lefsky et al., 1999; Richardson et al., 2009; Boudreault et al., 2015). To facilitate usage of ALS forest data in flow models with different types or levels of canopy dragforce prescription (e.g. linearized, RANS, and mesoscale models), Sogachev et al. (2017) recently suggested an objective method to translate PAD (momentum absorbed mostly by foliage distributed in the vertical) into spatially varying values of effective roughness, z_{0,eff} (momentum absorbed at one vertical level). The effective roughness can be used to connect the surface friction velocity and geostrophic wind speed in a way satisfying the geostrophic drag law (Blackadar and Tennekes, 1968); it also allows prediction of wind flow over vegetation using typical numerical models, without needing to incorporate local drag forces in each grid volume of a threedimensional model domain (see Fig. 1). In the present study, we explore the effect of different canopy representations on flow modelling over two forested sites: one located on a flat coastal plain (Østerild, Denmark) and the other in a hilly region (Rödeser Berg, Germany).
2.1 Models
Forests significantly reduce wind speeds and increase turbulence intensity. To address this, DTU Wind Energy has been developing RANS parameterizations for draginduced effects of vegetation for many years. The atmospheric boundary layer (ABL) model SCADIS (Sogachev et al., 2002, 2012), based on a finitedifference numerical approach, is a basic tool for investigation of vegetation effects on atmospheric flow. Due to its flexibility, SCADIS allows testing of any assumption in a straightforward way. However, the model was designed for a single processor and cannot yet be used on the massive computing scale, needed for the wind industry. However, DTU Wind Energy's other inhouse RANS (finitevolume) solver, EllipSys3D (e.g. Michelsen, 1992; Sørensen, 2003), was programmed for massively parallel computing, though it was not designed for the ABL. Based on fundamental investigations and research using SCADIS, EllipSys3D has been improved in order to simulate forested ABL flow problems. Both RANS models apply the twoequation closure approach based on coupled transport equations for the turbulent kinetic energy (TKE) k, and a supplementary lengthscale determining variable φ (which can be e.g. ε or ω, where ε is the dissipation rate of TKE and $\mathit{\omega}\equiv \mathit{\epsilon}/k$) (e.g. Pope, 2000).
2.2 Vegetation description in models
In the models with canopy dragforce prescriptions, detailed information about PAD can be directly used after modification of the RANS equations. Modification usually includes additional terms in the momentum and turbulence closure equations:

in the momentum (or velocity) component equations, for each spatial coordinate (direction) x_{i}, the gridvolume average momentum sink of canopy elements is usually parameterized by Raupach and Shaw (1982)
$$\begin{array}{}\text{(1)}& {F}_{i}={c}_{\mathrm{d}}\mathrm{PAD}\langle {\stackrel{\mathrm{\u203e}}{U}}_{i}\rangle S.\end{array}$$Here $S=(\langle {\stackrel{\mathrm{\u203e}}{U}}_{i}\rangle \langle {\stackrel{\mathrm{\u203e}}{U}}_{i}\rangle {)}^{\mathrm{1}/\mathrm{2}}\phantom{\rule{0.125em}{0ex}}$ is the magnitude of the spatial (angle brackets) and temporal (overbar) average of the wind speed, and c_{d} is the drag coefficient including the shelter effect.

In the turbulence closure equations (i.e. for k and φ), a wake production term and enhanced dissipation term due to the surface drag of canopy elements (e. g. Raupach and Shaw, 1982; Finnigan, 2000) is added. In contrast with conventional parameterizations for vegetationinduced momentum sinks, modification of the closure equations is still debated somewhat. In the present study, we use the canopy parameterizations described in Sogachev et al. (2012, 2015).
With an effective roughness approach, no modification of the RANS equations is needed; use of z_{0,eff} allows us to predict airflow over a (displaced) surface covered by the vegetation using a flow model, without explicitly resolving the canopy drag. A method for deriving z_{0,eff} from frontal area density consists of three basic steps (Sogachev et al., 2017):

using a statistical approach to rank frontal density data into landuse classes;

applying empirical methods to describe the mean properties (per each class):

utilizing the RANS solver to estimate the aerodynamic roughness, z_{0,eff} and as such to effectively describe each landuse class with a single parameter.
The SCADIS model using kε closure is employed in the last step here. SCADIS provides steadystate profiles of velocity components, turbulent kinetic energy, dissipation rate and the mean momentum flux for a specified geostrophic wind G, Coriolis parameter (i.e. latitude), and predefined canopy drag. Thus, a “geostrophic drag law” (e.g. Blackadar and Tennekes, 1968) is implicitly applied, which, combined with Rossby number similarity theory, gives the following geostrophic wind in the neutral, barotropic boundary layer:
Thus
In Eqs. (2) and (3) κ is the von Kármán constant (usually accepted to be 0.4), f is the Coriolis parameter, ${u}_{\ast ,\mathrm{0}}\phantom{\rule{0.125em}{0ex}}$is the friction velocity at the “virtual” surface (generally the mean treetop height). The geostrophic component coefficients A, B have historically been established empirically via ABL measurements (Melgarejo and Deardorff, 1974; Clarke and Hess, 1974). Our calculations used A=1.8 as in Troen and Petersen (1989) and wind energy practice (roughly equal to the mean of values found in literature). The value B=6.4 is fitted to match solutions of SCADIS with drag forces and z_{0,eff} and is slightly higher than the value B=4.5 given by Troen and Peterson (1989). In all calculations, G is taken as 15 m s^{−1} and f is specified in accordance to the location of the test site (see below).
2.3 Numerical experiments
The RANS solver EllipSys3D, which uses a kε turbulence model, was applied for numerical simulation of the wind field. For the present simulations, a polar (cylindrical) computational domain is chosen. The two main reasons for selecting a cylindrical domain are that the finest grid resolution is provided at the central part of the domain where it is most needed, and that the solution mesh is equally suited for ?ow from all directions. In the central part of the domain, a square grid zone is used around the target area. For both sites, the computational domain has a radius of about 25 km and a height of about 5 km. The target area for Østerild is about 10 km by 10 km with a horizontal resolution of 16 m (see below). The target area for Rödeser Berg is about 8 km by 8 km with a horizontal resolution of 10 m (see below). The vertical grid consists of nodes with variable size (i.e. about 1 cm near the ground and with maximal vertical resolution of 1 m in first 50 m above ground) allowing the canopy drag to be properly resolved. Details about modelling of complex terrain flow using EllipSys3D can be found in Koblitz (2013) and Sørensen et al. (2017). Calculations were done for the neutral flow using a steadystate algorithm for three geostrophic wind directions for Østerild site (90, 200, and 290^{∘}) and one geostrophic wind direction (210^{∘}) for Rödeser Berg site both with a full description of canopy drag forces and with forest replaced by the effective roughness. Note that wind direction in the surface layer is turned in respect to the geostrophic wind direction by approximately 20^{∘} counter clockwise. Thus, modelled surface wind directions are circa 70, 180, and 270^{∘} for Østerild site and circa 190^{∘} for Rödeser Berg site. In these calculations, G is taken as 15 m s^{−1}, to provide calculated wind speed magnitude in surface boundary layer closer to that observed here.
3.1 Østerild site
The Østerild site is located on a coastal plain, in northern Jutland, Denmark. Figure 2a shows that the terrain is flat, except for a small hill in the southwest direction and some sand dunes to the north. The site has grasslands, and forests with scotch pine, mountain pine, silver fir, Sitka spruce, oak, and birch. Most of the tree heights are between 15 and 20 m. The PAD for this site was derived from ALS data using the method of Boudreault et al. (2015). Figure 2b shows the estimated values of effective roughness over the target domain. The original elevation data set had a horizontal resolution of 1.6 m. To demonstrate the method and to speed up data processing, here we use a resolution of 16 m. There are multiple sets of observational data from the Østerild site, but only data derived from the 250 m tall met mast equipped with sonic anemometers at 5 height levels (for detail see Peña, 2019) are used in this study. The selection of date was done following Dellwik et al. (2014) recommendations.
3.2 Rödeser Berg site
The Rödeser Berg site is located in central Germany (in the northern part of the federal state of Hesse). Figure 3a shows that the terrain in this area is hilly, with the elevation reaching from 200 to 500 m. The forest type at the Rödeser Berg is a mixed forest, consisting mainly of beech trees, pine trees, and spruces, with a few areas of larch and oak. Less than 1 % of the trees in the area are taller than 40 m, and most of the trees have heights between 25 and 32 m. There are several clearings, forest roads and some scattered houses close to the forest edges. Details about ALS data and how PAD values were produced for this area are given in (Freier, 2017). Figure 3b shows the estimated values of effective roughness over the target domain. There are several met measurement facilities on this site (for detail see Mann et al., 2017). Only data from the 200 m tall met mast equipped with sonic anemometers at 9 height levels and located at the top of Rödeser Berg are used in the current study. Relatively few cases were selected, and the sample size was too small for reliable uncertainty estimates.
4.1 Coastal plain
Figure 4 compares the normalized wind speed profiles derived from measurements with numerical simulations, for the mast position at Østerild. Measured data are presented for the directional sectors of 30^{∘} wide closest to modelled wind directions. Figure 4 demonstrates the model's ability to qualitatively reproduce wind speed over the forested area. Good agreements between modelled data and observations are noted for easterly (Fig. 4a) and southerly (Fig. 4b) winds. For westerly winds (Fig. 4c), the RANS simulations provide wind speed profiles laying between those that were observed for two adjacent wind sectors. A possible reason for the scatter in observations with wind direction is the large clearcut, W/SW of the mast (see Fig. 2b). A minimal difference in wind speed profiles is found between the two different canopy approaches for easterly winds; this is probably because the longest forested fetch exists towards the east (see Fig. 2b). Small differences in wind speed profiles derived by the canopy approaches for southerly and westerly winds can be explained by the presence of clearcuts of different sizes, and limited fetch. In both cases, the difference in data derived by two approaches decreases for heights farther above the ground. No difference is observed above a height of ∼100 m.
Figure 5 shows wind speeds (S) and local friction velocity (u_{∗}) along the northsouth transect over the area indicated in Fig. 2b (dashed line), modelled using the two canopy approaches for southerly winds. The results of the numerical experiments presented in Fig. 5 show that models with different parameterization of the vegetation effect produce similar wind speeds (Fig. 5a, c) and local friction velocities (Fig. 5b, d) over the investigated domain. The effects of vegetation on the airflow are clearly seen. When using the canopy drag parameterization, the model provides more details on airflow near the ground, but with increasing the height flow parameters become smoothed. One can see that the large difference in estimated wind speed fields (Fig. 5e) as well as in estimated local friction velocities (Fig. 5f) is observed in the lowest layer of atmosphere and of course inside the layer occupied by vegetation. Data presented in Fig. 5e indicate that the largest difference in wind speed occurs in the zones of forest edge, especially when the wind hit the forest from an open place with large upwind fetch.
To exclude the local effect of the underlying surface on wind profiles we consider average values of flow parameters over Østerild target area at different heights. Figure 6 shows that the model with the different parameterizations of the vegetation effect produces similar results, except the effective roughness approach provides a smaller scatter of modelled variables around its average value.
4.2 Hilly terrain
Figure 7 compares the normalized wind speed profiles derived from measurements and numerical simulations at the mast position at the Rödeser Berg site. Observed data are scattered because they are for different wind directions inside the southwestern sector (180–225^{∘}). The model does a good job in reproducing the wind speed profile with the canopy drag prescription approach, but fails with the effective roughness approach. Figure 8 shows results along the northsouth transect over the area indicated in Fig. 3b (dashed line) modelled using different canopy approaches for the southern wind sector. One can see that the two canopy model implementations give similar airflow only upwind of the hill; on the summit and downwind, the modelled airflow structures are different. The approach based on effective roughness results in higher wind speed (Fig. 8a, see also Fig. 7) above topographic peaks, in comparison with the full drag approach (Fig. 8c, e). The roughnessbased approach overestimates TKE values (consideration of the local friction velocity in our opinion has limited quantitative value regarding complex terrain flow) near the ground (again especially at topography summits), but underestimates TKE values at higher altitudes, and especially downwind of topography summits (Fig. 8b, d and f). Analysis of average values of flow parameters over Rödeser Berg target area at different heights (Fig. 9) retrievals that the roughness model approach in general provide higher wind speed (Fig. 9a, b) and lower TKE levels (Fig. 9c) than provided by the full drag model approach. The roughness model approach also provides a narrow scatter of modelled variables around its average value.
Targeted numerical experiments and their comparison with observations indicate that using a roughness representation of forest within typical RANS flow modelling, one can estimate wind speed and turbulence levels over forested areas, at heights potentially interesting for wind energy applications (∼60 m and higher) – but only above the flat terrain. Not surprisingly, the roughnessbased approach becomes less valid near forest edges, and more valid in directions with longer forest fetch; for the former limitation, there is consistency between our results and the linearized modelling attempts made by Bingöl (2010). Advanced flow models capable of incorporating (distributed) local drag forces are recommended for complex terrain covered by forest, and the RANS turbulence model likely needs to be extended to predict flow over hilly forested terrain.
All simulations used in this paper are available via request to the corresponding author. Datasets for observations are currently not publicly available.
AS and MK conceptualized the research goal and designed the numerical experiments. ED provided forest data and observations for Østerild site, TK and PK provided related data for Rödeser Berg site. DC performed the simulations. AS prepared the manuscript with contributions from all coauthors.
The authors declare that they have no conflict of interest.
We thank Julia Freier from Fraunhofer IEE, Kassel, Germany for providing elevation map for Rödeser Berg site. We thank The New European Wind Atlas (NEWA) project consortium for making their experimental data available at a stage before the data became publicly accessible.
This article is part of the special issue “19th EMS Annual Meeting: European Conference for Applied Meteorology and Climatology 2019”. It is a result of the EMS Annual Meeting: European Conference for Applied Meteorology and Climatology 2019, Lyngby, Denmark, 9–13 September 2019.
This paper was edited by SvenErik Gryning and reviewed by Lars Landberg and Javier Sanz Rodrigo.
Bingöl, F.: Complex Terrain and Wind Lidars. Roskilde, Risø National Laboratory for Sustainable Energy, RisøPhD, No. 52(EN), 60 pp., 2010.
Blackadar, A. K. and Tennekes, H.: Asymptotic similarity in neutral barotropic planetary boundary layers, J. Atmos. Sci., 25, 1015–20, https://doi.org/10.1175/15200469(1968)025<1015:ASINBP>2.0.CO;2, 1968.
Boudreault, L.E.: Reynoldsaveraged NavierStokes and LargeEddy Simulation over and inside inhomogeneous forests, DTU Wind Energy, PhD0042 (EN), 118 pp., 2015.
Boudreault, L.É., Bechmann, A., Tarvainen, L., Klemedtsson, L., Shendryk, I., and Dellwik, E.: A LiDAR method of canopy structure retrieval for wind modeling of heterogeneous forests, Agr. Forest Meteorol., 201, 86–97, https://doi.org/10.1016/j.agrformet.2014.10.014, 2015.
Clarke, R. H. and Hess, G. D.: Geostrophic departure and the functions A and B of Rossbynumber similarity theory, Bound.Lay. Meteorol., 7, 261–287, https://doi.org/10.1007/BF00240832, 1974.
Dellwik, E., Bingöl, F., and Mann, J.: Flow distortion at a dense forest edge, Q. J. Roy. Meteor. Soc., 140, 676–686, https://doi.org/10.1002/qj.2155, 2014.
Finnigan, J. J.: Turbulence in plant canopies, Annu. Rev. Fluid. Mech., 32, 519–571, https://doi.org/10.1146/annurev.fluid.32.1.519, 2000.
Freier, J.: Approaches to characterise forest structures for wind resource assessment using airborne laser scan data, Master thesis, Department of Mechanical Engineering, Kassel University, Kassel, 111 pp., 2017.
HVBG, Hessisches Landesamt für Bodenmanagement und Geoinformation: Landesweites Laserscanning Hessen – Die Basis für ein Solarkataster, with assistance of Dorn, C. and Köhler, G., Hessian Administration Office of Land Management and Geological Information, Wiesbaden (DE), 2010.
Koblitz, T.: CFD modelling of nonneutral atmospheric boundary layer conditions, DTU Wind Energy PhD0019 (EN), 104 pp., 2013.
Lefsky, M. A., Cohen, W. B., Acker, S., Parker, G. G., Spies, T. A., and Harding, D.: Lidar remote sensing of the canopy structure and biophysical properties of Douglasfir western hemlock forests, Remote Sens. Environ., 70, 339–361, https://doi.org/10.1191/0309133303pp360ra, 1999.
Mann, J., Angelou, N., Arnqvist, J., Callies, D., Cantero, E., Arroyo, R. C., Courtney, M., Cuxart, J., Dellwik, E., Gottschall, J., Ivanell, S., Kuhn, P., Lea, G., Matos, C. J. S., Palma, J. M. L. M., Pauscher, L., Peña, A., Rodrigo, J. S., Söderberg, S., Vasiljevic, N., and Rodrigues, C. V.: Complex terrain experiments in the New European Wind Atlas, Philos. T. R. Soc. A, 375, https://doi.org/10.1098/rsta.2016.0101, 2017.
Melgarejo, J. W. and Deardorff, J. W.: Stability functions for the boundary layer resistance laws based upon observed boundary layer height, J. Atmos. Sci., 31, 1324–1333, https://doi.org/10.1175/15200469(1974)031<1324:SFFTBL>2.0.CO;2, 1974.
Michelsen, J. A.: Basis3d – a platform for development of multiblock pde solvers, Technical report AFM9205, Technical University of Denmark, 1992.
Peña, A.: Østerild: A natural laboratory for atmospheric turbulence, J. Renewable Sustainable Energy, 11, 063302, https://doi.org/10.1063/1.5121486, 2019.
Pope, S. B.: Turbulent flows, Cambridge University Press, UK, 771 pp., 2000.
Raupach, M. R. and Shaw, R. H.: Averaging procedures for flow within vegetation canopies, Bound.Lay. Meteorol., 22, 79–90, https://doi.org/10.1007/BF00128057, 1982.
Richardson, J. J., Moskal, L. M., and Kim, S. H.: Modeling approaches to estimate effective leaf area index from aerial discretereturn LIDAR, Agr. Forest Meterol., 149, 1152–1160, https://doi.org/10.1016/j.agrformet.2009.02.007, 2009.
Sogachev, A., Menzhulin, G. V., Heimann M., and Lloyd, J.: A simple threedimensional canopy – planetary boundary layer simulation model for scalar concentrations and fluxes, Tellus B, 54, 784–819, https://doi.org/10.3402/tellusb.v54i5.16729, 2002.
Sogachev, A., Kelly, M., and Leclerc, M. Y.: Consistent twoequation closure modelling for atmospheric research: buoyancy and vegetation implementations, Bound.Lay. Meteorol., 145, 307–327, https://doi.org/10.1007/s1054601297265, 2012.
Sogachev, A., Cavar, D., Bechmann, A., and Jørgensen, H. E.: Assessment of consistent twoequation closure for forest flows, Paper presented at EWEA Annual Conference and Exhibition, Paris, France, 2015.
Sogachev, A., Cavar, D., Kelly, M., and Bechmann, A.: Effective roughness and displacement height over forested areas, via reduceddimension CFD, Report of DTU Wind Energy – E, vol. 0161, 36 pp., 2017.
Sørensen, N. N.: Technical report, RisøR827(EN), Risø National Laboratory, 154 pp., 2003.
Sørensen, N. N., Bechmann, A., Boudreault, LE., Koblitz, T., Sogachev, A., and Laan, M.: CFD for Atmospheric Flows and Wind Engineering, Lecture Series 20162017, von Karman Institute for Fluid Dynamics, Belgium, https://doi.org/10.35294/ls201701, 2017.
Troen, I. and Petersen, E. L.: European Wind Atlas. Roskilde: Risø National Laboratory, Denmark, 656 pp., 1989.