Numerical modelling of the wind over forests: roughness versus canopy drag

Parameterizing the effect of vertically-distributed vegetation through an effective roughness (z0,eff) – whereby momentum loss through a three-dimensional 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 z0,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.


Introduction
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 (Reynolds-Averaged Navier-Stokes) solvers have become a commonly-used 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 drag-force 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 ver-tical) 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 three-dimensional 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).

Models
Forests significantly reduce wind speeds and increase turbulence intensity. To address this, DTU Wind Energy has been developing RANS parameterizations for drag-induced effects of vegetation for many years. The atmospheric boundary layer (ABL) model SCADIS (Sogachev et al., 2002(Sogachev et al., , 2012, based on a finite-difference numerical approach, is a basic tool for investigation of vegetation effects on atmo- spheric 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 in-house RANS (finite-volume) 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 length-scale determining variable ϕ (which can be e.g. ε or ω, where ε is the dissipation rate of TKE and ω ≡ ε/k) (e.g. Pope, 2000).

Vegetation description in models
In the models with canopy drag-force 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: 1. 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) Here S = ( U i U i ) 1/2 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.
2. 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 vegetation-induced 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. (2012Sogachev et al. ( , 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 : 1. using a statistical approach to rank frontal density data into land-use classes; 2. applying empirical methods to describe the mean properties (per each class): 3. utilizing the RANS solver to estimate the aerodynamic roughness, z 0,eff and as such to effectively describe each land-use class with a single parameter.
The SCADIS model using k-ε closure is employed in the last step here. SCADIS provides steady-state 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 * ,0 is the friction velocity at the "virtual" surface (generally the mean tree-top 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).

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)

Ø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.

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. 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 clear-cut, 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 clear-cuts 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 north-south 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 param-eters 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  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. 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 north-south 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 roughness-based 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.

Conclusions
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 sur-prisingly, the roughness-based 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