Micro-scale model comparison (benchmark) at the moderately complex forested site Ryningsnäs

This article describes a study in which modellers were challenged to compute the wind field at a forested site with moderately complex topography. The task was to model the wind field in stationary conditions with neutral stratification by using the wind velocity measured at 100 m at a metmast as the only reference. Detailed maps of terrain elevation and forest densities were provided as the only inputs, derived from airborne laser scans (ALSs) with a resolution of 10 m × 10 m covering an area of 50 km × 50 km, that closely match the actual forest and elevation of the site. The participants were free to apply their best practices for the simulation to decide the size of the domain, the value of the geostrophic wind, and every other modelling parameter. The comparison of the results with the measurements is shown for the vertical profiles of wind speed, shear, wind direction, and turbulent kinetic energy. The ALS-based data resulted in reasonable agreement of the wind profile and turbulence magnitude. The best performance was found to be that of large-eddy simulations using a very large domain. For the Reynolds-averaged Navier–Stokes type of models, the constants in the turbulence closure were shown to have a great influence on the yielded turbulence level, but were of much less importance for the wind speed profile. Of the variety of closure constants used by the participating modellers, the closure constants from Sogachev and Panferov (2006) proved to agree best with the measurements. Particularly the use of C_μ≈0.03 in the k–ε model obtained better agreement with turbulence level measurements. All except two participating models used the full detailed ground and forest information to model the forest, which is considered significant progress compared to previous conventional approaches. Overall, the article gives an overview of how well different types of models are able to capture the flow physics at a moderately complex forested site.

1 Introduction

To respond to the increasing demand for wind power, new areas for wind turbine sites are being explored. Large offshore farms further away from the shore are being developed, as are wind farms in more complex onshore areas, such as terrain with a more varied topography and roughness. This is the case in northern countries, such as the Scandinavian region, where large remote forested areas are being explored for wind development. However, when exploring these complex sites it is evident that new challenges arise due to comparatively higher turbulence levels and wind shear. While the magnitude of wind shear and turbulence increase the fatigue load, uncertainties in the estimation of wind shear and turbulence have shown to be problematic in forested areas (Enevoldsen, 2016). Hence, it is important to assess the uncertainties in the modelling process of wind conditions over forests.

In addition to the actual difference in wind climate between traditional wind energy sites and complex forested ones, modelling of the wind conditions is challenging. Trees are elevated sources for both momentum absorption and heat transfer, and thus they differ from traditional surfaces since the exchange may be distributed at several model levels. The degree of physical description is a choice by the modeller, from describing plant area densities (PADs) in each grid cell to representing an entire forest with a single roughness length value. The required numerical demand, however, varies with many orders of magnitude when making that choice.

To the knowledge of the authors, no large-scale studies have been published comparing different micro-scale models over forested terrain with high-quality meteorological data. However, Ayotte (2008) compared models of varying complexity to wind tunnel measurements and concluded that inaccurate representation of all physical scales may result in significant errors. Earlier model intercomparison studies of micro-scale models in non-forested areas have provided insight into the performance of different model types and highlighted important differences in modelling choices such as closure constants (Bechmann et al., 2011; Bosveld et al., 2014). So far there are significant uncertainties in how well wind climate models perform in forested areas, and there are also large differences among model descriptions. Hence, there is a need for more validation studies and a better understanding of how the different modelling choices affect the final result. This study aims to take the first steps to fill that knowledge gap by presenting model performances at a forested wind turbine site.

The progress of forest flow modelling now enables the direct simulation of the tree densities. Such values as PAD may be derived from airborne laser scans (ALSs) that are becoming increasingly available from national mapping services (Boudreault et al., 2015). Using PAD data instead of estimated roughness lengths may be a way to reduce the uncertainties of site assessment. However, the use of drag modelling through PAD is increasingly being adopted in the wind energy community, and how to best make use of the data is still an open question, as is how their use affects the model abilities. In order to the test the performance of the wind simulation determined by models using PAD derived from ALS, data corresponding to this quantity were made available to the modellers taking part in the study.

The study started with a call for a benchmarking model validation study to modellers involved in the European ERANET+ project New European Wind Atlas (NEWA). The aim of the benchmark is to illustrate how well micro-scale models are able to simulate winds above a forest in moderately complex topography. The participating models range from industrial wind models to front-line research models. The modelled case consists of a typical site located in Ryningsnäs in southern Sweden, i.e. a patchy forested site with moderately complex topography (Arnqvist et al., 2015). The benchmark study is in that sense quite broad and aims to answer the question of how successful different modelling programmes, as well as modelling practices, are in matching the measured wind profile in specified conditions. The instructions for the benchmark were deliberately kept rather open in order to study how the different groups made use of the PAD data and how the compromise between resolution and domain size was handled.

The NEWA project includes several large-scale field campaigns designed for flow model validation (Mann et al., 2017); however, the Ryningsnäs measurement campaign was performed prior to the start of NEWA project and was identified as an appropriate data set for a benchmarking study. As such, it also forms a basis for model validation methodology as preparation for upcoming micro-scale benchmarks (Mann et al., 2017) using measurement input from the extensive measurement campaigns performed within the NEWA project.

The paper begins by outlining the following: the benchmark, the validation data, and general modelling. This is followed by a description of first the RANS models and then the LES models. It then continues with the main results and finally concludes with a “Discussion and conclusions” section.

2 Benchmark description

The benchmark task was to model the wind profile at the location $\mathrm{57}{}^{\circ }{\mathrm{16}}^{\prime }{\mathrm{34.26}}^{\prime \prime }$ N, $\mathrm{15}{}^{\circ }{\mathrm{59}}^{\prime }{\mathrm{12.23}}^{\prime \prime }$ E for the wind directions 100, 240, and 290^∘ (directions at 100 m of height). The input provided to the modellers was a target wind speed of 7.4 m s⁻¹ at 100 m of height, neutral atmospheric stratification, and a data set of forest density and topography in a 40 km × 40 km grid. The modellers were asked to provide the wind profile from the ground up to at least 200 m, geostrophic wind speed U_g, wind speed in planes at 40, 100, and 140 m above local ground level (a.g.l.), and information about their model.

The choice behind having a target wind speed at 100 m of height, rather than having a fixed geostrophic wind speed, was partly due to the lack of measurements of the geostrophic wind speed and partly due to the desire to have as similar a wind speed as possible in the lower part of the boundary layer. This relates to the question of whether or not the models can accurately predict the flow footprint given that the ALS data enable the model to have surface conditions very similar to reality. Setting a fixed geostrophic wind speed would risk the modelled wind speed in the surface layer being lower or higher than the measured, with a subsequent uncertainty of the ability to capture the footprint of the flow. In a strictly neutral boundary layer, the ratio of the turbulence level and the wind speed is expected to be constant, and hence the footprint would be the same for different wind speeds. However, the boundary layer height changes with wind speed, as do the gradients of velocity and turbulence level. Thus, in addition to scaling the wind speed with the friction velocity, one would have to scale the height with the boundary layer height in order to make a fair comparison, which cannot be done since the boundary layer height was not measured.

3 Measurements

3.1 Forest characterization by laser scans

In order to characterize the forest ALS data from the Swedish map authority, Lantmäteriet has been utilized (Lantmäteriet, 2016). The data were collected at a flight height of 1700 m, yielding a footprint area at the ground of 0.5 m² for the laser beam. The density is around 1 shot m⁻². The data were processed according to the method described in Boudreault (2015). The method uses the Beer–Lambert law for the attenuation of the laser beam as it travels through the forest canopy, and the plant area density (PAD) can be derived if the extinction coefficient is known. The extinction coefficient was here (Boudreault, 2015) assumed to be 0.5∕cos(θ_ALS), where θ_ALS is the scanning angle corresponding to a spherical distribution of canopy elements. The height of the forest was derived by the maximum return height in a grid box as defined by the distance from the median of the returns in that grid box that had been classified as ground (the ground height of the grid box). The PAD was derived in vertical steps of 1 m from the highest return reflection towards the ground. In order to avoid numerical problems when the beam becomes fully attenuated (in very dense forest patches) PAD estimation was set to a constant of 0.1 m⁻¹ below the level of 95 % attenuation. Since most of the returns normally come from the ground, this filter only affected 0.1 % of the PAD estimates. Two data sets with a resolution of 10 m × 10 m and 50 m × 50 m, respectively, were then prepared to be used as model input. The data sets include horizontal coordinates, ground height, tree height, and PAD in a vertical grid from 0 m above ground to the tree height in steps of 1 m.

Figure 1 (a–c) Forest cover of the Ryningsnäs site with increasing magnification factors. The contour maps are in greyscale and the colour scales illustrate the tree height distribution, both inside and outside the selected area. (d) A topographic map of the Ryningsnäs area.

3.2 Site description

The measurement site is located in Ryningsnäs in south-east Sweden ($\mathrm{57}{}^{\circ }{\mathrm{16}}^{\prime }{\mathrm{34.26}}^{\prime \prime }$ N, $\mathrm{15}{}^{\circ }{\mathrm{59}}^{\prime }{\mathrm{12.23}}^{\prime \prime }$ E). A meteorological measurement tower equipped with cup and sonic anemometers, with the highest measurement level at 138 m above ground, has been used for the validation data. The area around the tower site has moderate complexity in terms of topography, but the forest cover is very heterogeneous with many clearings and forest stands of different ages. Figure 1 shows the forest cover on different scales, the largest being 50 km × 50 km and the smallest 1 km × 1 km. The three sectors chosen for the validation study have been coloured.

The tower is situated in the north-west corner of a 400 m × 400 m clearing. The surrounding forest has a peak in the tree height distribution at approximately 20 m and predominantly consists of Scots pine (Pinus sylvestris). The actual tree height distribution can be seen in Fig. 2a and b in which the distribution is shown for the three sectors within a radius of 10 km (panel a) and 1 km (panel b). In the larger scale, the tree height distributions of sectors 100 and 240^∘ are very similar and have a peak at 21.5 m. The distribution of sector 290^∘ is more flat, with two maximums, one at 7.5 m and one at 19.5 m. In the closer region the 290^∘ distribution is different from the two other sectors in that it does not contain a clearing, but does contain a large patch of young forest at 7.5 m of height. The clearing of the tower itself can be seen in the distribution of the 100^∘ sector, whereas the clearing in the 240^∘ sector is almost 1 km upwind.

Two wind turbines are situated approximately 200 m to the north-east and south of the tower, respectively, but the three sectors used in the validation study exclude directions which these turbines would influence.

Figure 2 (a–c) Distribution of tree heights in the three sectors within a radius of (a) 10 km, (b) 2 km, and (c) 1 km. The colours indicate sectors 100^∘ (blue), 240^∘ (green), and 290^∘ (red). (d) Average PAD profiles from within 10 km are shown as full lines, 2 km as dotted lines, and 1 km as dashed lines. Colouring is as in the histograms.

3.3 Wind measurements

The full measurement set-up and the wind climate has been reported earlier in Arnqvist et al. (2015). The instruments used in this study include six sonic anemometers (Metek Gmbh, USA-1) located at heights of 40, 59, 80, 98, 120, and 137.7 m as well as seven Thies first-class cup anemometers at heights of 25.5, 40.1, 60.5, 80.1, 95.85, 120.75, and 137.6 m.

The sonic anemometers were sampled at 20 Hz and statistics was evaluated by 30 min block averaging and 3-D rotation of the coordinate system, aligning it with the local mean wind direction and yielding the wind vector u, v, w and temperature t. Flow distortion correction and quality checks were applied as in Arnqvist et al. (2015). The stratification was evaluated by means of $(z-d)/L$, where z is the height, d is the displacement height, and L, the Obukhov length, was determined as

$$\begin{array}{}\text{(1)}& L=-{\displaystyle \frac{u_{*}^{\mathrm{3}}T}{\mathit{\kappa }g\stackrel{\mathrm{‾}}{w^{\prime }{t}^{\prime }}}},\end{array}$$

where the velocity fluctuations are defined by the substraction of the mean velocities from the instantaneous values $u^{\prime }=u-\stackrel{\mathrm{‾}}{u}$, $v^{\prime }=v-\stackrel{\mathrm{‾}}{v}$, $w^{\prime }=w-\stackrel{\mathrm{‾}}{w}$; $u_{*}=({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{4}}$ (the friction velocity), κ=0.4 is the von Kármán constant, g=9.81 m s⁻¹ is the gravitational acceleration, and $t^{\prime }=t-\stackrel{\mathrm{‾}}{t}$ is the instantaneous temperature fluctuation.

To select only neutral conditions $(z-d)/L$ was required to be between −0.1 and 0.07 at all heights. The limits were selected based on the shape of the ϕ function for momentum (Högström, 1996) and allows for a roughly ±35 % variation in the wind gradient given a certain u_* value. In addition, the 100 m wind speed was required to be between 7 and 8 m s⁻¹ and the conditions quasi-stationary, as defined by the fact that the wind speed was allowed to vary a maximum of 10 % between adjacent 10 min segments and the wind direction was allowed to vary a maximum of 10^∘ between adjacent 30 min segments. The wind directions were required to be within ±10^∘ of the target wind direction. After applying all the conditions and quality controlling the data, 9, 13, and 9 separate 30 min segments remained in sectors 100, 240, and 290^∘, respectively. To sum up, the data selection consisted of

• a quality check that was passed;

• neutral stratification;

• stationary flow;

• $(\mathrm{7}<{\stackrel{\mathrm{‾}}{u}}_{\mathrm{98}\text{m}}<\mathrm{8})$ m s⁻¹; and

• wind direction within the target sector.

The data that fulfilled the neutral and stationary conditions constitute 10 % of the total data set and the ones also fulfilling $(\mathrm{7}<{\stackrel{\mathrm{‾}}{u}}_{\mathrm{98}\text{m}}<\mathrm{8})$ m s⁻¹ make up 1 % of the entire data set. Selecting for the three specific directions further reduces the number of data points.

Table 1 Overview of the models used and the model family.

^* PALM used mean values of the ALS-derived PAD combined with flat terrain.

4 Modelling

The models that participated in the benchmark were all CFD models using a Reynolds-averaged Navier–Stokes (RANS) or large-eddy simulation (LES) methodology. Table 1 shows an overview of the models used by the respective participants. Some of the models include a full topography and PAD description. All models use a drag formulation to simulate the forest, with a drag coefficient of 0.2 (except Meteodyn which uses another drag formulation; please see Sect. 4.2.2). In the following section a description of each model set-up is given.

4.1 General CFD modelling

Computational modelling of the fluid flow employs a filtered version of the Navier–Stokes equations due to the impracticality of resolving every temporal and spatial scale. The Reynolds-averaged Navier–Stokes (RANS) equations make use of the Reynolds decomposition to divide the velocity field into the time-averaged velocity and the velocity fluctuation around the mean, $u_i=\stackrel{\mathrm{‾}}{u_i}+u_i^{\prime }$. This yields a momentum equation for the mean flow in which the effect of the turbulent motions is represented via the Reynolds stresses $\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}$, requiring a model to represent their effect on the average field. Most approaches (known as turbulence viscosity models) employ the Boussinesq approximation in which the Reynolds stresses are parameterized as a function of an eddy viscosity ν_t and a rate of strain tensor, which assumes that the turbulence fluxes are proportional to the mean velocity gradient. This yields

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial t}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}\stackrel{\mathrm{‾}}{u_j}}{\partial x_j}}{\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mathit{\nu }+{\mathit{\nu }}_t\right)\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right)\right]\\ {\displaystyle }& {\displaystyle }+f_{\text{c}}{\mathit{\epsilon }}_{ij\mathrm{3}}\stackrel{\mathrm{‾}}{u_j}+\stackrel{\mathrm{‾}}{f_i}\end{array}$$

as the momentum equation, where f_i represents an external body force. The Coriolis force is included using the Coriolis parameter f_c=2Ωsinλ, where Ω is the Earth's rotational velocity and λ the latitude of the wind farm. The eddy viscosity is modelled through the introduction of transport equations, such as in the frequently employed k–ε technique.

RANS modelling supposes that the effect of all ranges of fluctuations on the mean flow can be accounted for by the models. Conversely, in the LES approach the energy-containing flow structures are fully resolved, whereas only the effect of the smaller fluctuations is modelled. This is achieved through the decomposition of the velocity field into filtered (or resolved) and residual (or subgrid-scale, SGS) components, $u_i=\langle u{\rangle }_i+u_i^{\text{SGS}}$. Although various types of filters exist, a very common method in wind research is to associate a filter width with the grid spacing. The application of this decomposition on the Navier–Stokes equations leads to the apparition of the SGS stress tensor ${\mathit{\tau }}_{ij}^{\text{SGS}}\equiv -\mathit{\rho }\left(\langle u_i{u}_j\rangle -\langle u{\rangle }_i\langle u{\rangle }_j\right)$, which needs to be modelled. As in RANS, the prevalent strategy is to apply the Boussinesq approximation to introduce a subgrid viscosity ν_SGS to derive the LES momentum equation:

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial \langle u_i\rangle }{\partial t}}+{\displaystyle \frac{\partial \langle u_i\rangle \langle u_j\rangle }{\partial x_j}}{\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \langle p_{\text{m}}\rangle }{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mathit{\nu }+{\mathit{\nu }}_{\text{SGS}}\right)\left({\displaystyle \frac{\partial \langle u_i\rangle }{\partial x_j}}+{\displaystyle \frac{\partial \langle u_j\rangle }{\partial x_i}}\right)\right]\\ {\displaystyle }& {\displaystyle }+f_{\text{c}}{\mathit{\epsilon }}_{ij\mathrm{3}}\langle u_j\rangle +f_i,\end{array}$$

where p_m denotes the modified pressure, which includes the isotropic part of the SGS stress tensor. The simplest approaches to calculate ν_SGS (and amongst the most commonly employed in wind research) make use of the resolved scales.

4.2 RANS

The participating RANS models use one- or two-equation turbulence models, presented in general form in this section. The specific set-up of each model is presented in the following sections. The two-equation turbulence-closure model corresponds to the classical k–ε model (Launder and Spalding, 1974), in which transport equations for the turbulent kinetic energy k and its dissipation rate ε,

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial k}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\stackrel{\mathrm{‾}}{u_j}k\right)-{\displaystyle \frac{\partial }{\partial x_j}}\left(\left(\mathit{\nu }+{\displaystyle \frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_k}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right){\displaystyle }& {\displaystyle }={\mathit{\nu }}_t{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right)-\mathit{\epsilon }+S_k,\end{array}$$

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\epsilon }}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\stackrel{\mathrm{‾}}{u_j}\mathit{\epsilon }\right)-{\displaystyle \frac{\partial }{\partial x_j}}\left(\left(\mathit{\nu }+{\displaystyle \frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_{\mathit{\epsilon }}}}\right){\displaystyle \frac{\partial \mathit{\epsilon }}{\partial x_j}}\right){\displaystyle }& {\displaystyle }=C_{{\mathit{\epsilon }}_{\mathrm{1}}}{\displaystyle \frac{\mathit{\epsilon }}{k}}{\mathit{\nu }}_t{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right)-C_{{\mathit{\epsilon }}_{\mathrm{2}}}{\displaystyle \frac{{\mathit{\epsilon }}^{\mathrm{2}}}{k}}+S_{\mathit{\epsilon }},\end{array}$$

are solved. $C_{{\mathit{\epsilon }}_{\mathrm{1}}}$, $C_{{\mathit{\epsilon }}_{\mathrm{2}}}$, σ_k, and σ_ε are modelling constants, and S_k, S_ε are the source–sink terms representing the drag-based energy loss in the canopy. The eddy viscosity is determined from

$$\begin{array}{}\text{(6)}& {\mathit{\nu }}_t=C_{\mathit{\mu }}{\displaystyle \frac{k^{{}^{\mathrm{2}}}}{\mathit{\epsilon }}},\end{array}$$

where C_μ is another modelling constant. Following Sogachev and Panferov (2006) the source–sink canopy terms are parameterized as

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{S}_k=\mathrm{0},\text{(8)}& {\displaystyle }& {\displaystyle }{S}_{\mathit{\epsilon }}=\mathrm{12}\left(C_{{\mathit{\epsilon }}_{\mathrm{2}}}-C_{{\mathit{\epsilon }}_{\mathrm{1}}}\right)C_{\mathit{\mu }}^{\frac{\mathrm{1}}{\mathrm{2}}}{C}_{\text{d}}a\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|\mathit{\epsilon },\end{array}$$

where C_d is a drag coefficient and a= PAD is a frontal area density.

The length scale in the standard k–ε model is not bound and grows indefinitely with height. In order to adjust the model to ABL-relevant flow cases, a correction suggested by Apsley and Castro (1997) is applied, whereby the $C_{{\mathit{\epsilon }}_{\mathrm{1}}}^{{}^{*}}$ constant is (re)defined in the following manner:

$$\begin{array}{}\text{(9)}& C_{{\mathit{\epsilon }}_{\mathrm{1}}}^{{}^{*}}=C_{{\mathit{\epsilon }}_{\mathrm{1}}}+(C_{{\mathit{\epsilon }}_{\mathrm{2}}}-C_{{\mathit{\epsilon }}_{\mathrm{1}}}){\displaystyle \frac{l}{l_{max}}},\end{array}$$

with the mixing length l defined by

$$\begin{array}{}\text{(10)}& l=C_{\mathit{\mu }}^{{}^{\frac{\mathrm{3}}{\mathrm{4}}}}{\displaystyle \frac{k^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathit{\epsilon }}}.\end{array}$$

The limiting mixing length scale l_max is determined based on the relationship proposed by Blackadar (1962):

$$\begin{array}{}\text{(11)}& l_{max}={\displaystyle \frac{\mathrm{0.00027}G}{f_{\text{c}}}},\end{array}$$

where G is the geostrophic wind and f_c the Coriolis parameter.

The two-equation methodology explained above is used by Ellipsys3D, CFDWind, and Alya. In the case of Meteodyn, a one-equation RANS turbulence model k–l_m is instead used (Delaunay, 2007). This methodology consists of solving the turbulence kinetic energy (TKE) equation, Eq. (4), by replacing ε in terms of k and a parameterized mixing length l_m. Thus, $\mathit{\epsilon }=\mathit{\epsilon }(k,l_{\text{m}})$.

Furthermore, assuming the canopy elements exert a drag force on the flow, the effects of the plant drag inside the canopy on the main flow are parameterized, presuming the form drag dominance in the momentum Eq. (2) as

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{f_i}=\stackrel{\mathrm{‾}}{f_{\text{D},i}}=-C_{\text{d}}a\mid \stackrel{\mathrm{‾}}{\mathit{u}}\mid \stackrel{\mathrm{‾}}{u_i}.\end{array}$$

4.2.1 Meteodyn

Model description

Meteodyn WT is a commercial site-assessment software that models the surface boundary layer (no Coriolis force included) using RANS, in particular a one-equation k–l_m turbulence model and wall functions based on the Monin–Obukhov theory (Delaunay, 2007, 2013; Hurley, 1997). This turbulence-closure scheme uses a prognostic equation for the turbulent kinetic energy and a mixing length approach for estimating the turbulent diffusion. Meteodyn WT version 5.2.1 was used in the present investigation.

The forest model in Meteodyn is based on a mean flow model, which treats the forest as a porous media (Costa, 2007), similar to what is used in other commercial solvers like WindSim (Crasto, 2006). A volumetric sink term is introduced in the momentum equations for all cells lying inside or partially covering the forest volume. The volumetric force depends on the drag coefficient C_d, which is a function of the forest density. In the Meteodyn simulations presented in this paper we used a value of $C_{\text{d}}^{*}=\mathrm{0.005}$.

Table 2 Summary of the model constants used.

In Meteodyn the canopy height h for each cell is directly derived from the local roughness using the relation h=Cz₀, where C is an adjustable constant usually set to 20 and z₀ the roughness length at the cell surface. Cells with roughnesses higher than 0.8 cm are considered forest in the domain. Meteodyn can only handle roughness information provided in the .map format (contour line information as in the WAsP software; WAsP, 2018). For this reason, it was not possible to use the detailed PAD information provided to the modellers. The Meteodyn simulations presented here use roughness maps derived from SRTM data obtained using windPRO (WindPro, 2009).

Two versions of the forest model are available in Meteodyn, which differ in the computation of the mixing length in the one-equation turbulence model, called robust and dissipative in the Meteodyn documentation. The dissipative forest model is used in the Meteodyn simulations presented in this investigation, in which a 15 m extra high dissipation zone is used above the forest.

Numerical set-up

The computations are performed employing a Cartesian structured mesh on a square domain with the dimensions 13.5×13.5 km² and a 2.9 km height. The mesh is refined around the metmast location, with a grid-stretching factor of 1.1 in the horizontal and 1.2 in the vertical directions. The final mesh has 224 000 cells. On the vertical direction the lines are always orthogonal to the topography surface.

Monin–Obukhov inlet profiles for velocity are defined at the inlet of the domain, as is a constant turbulent kinetic energy (Richards and Hoxey, 1993). The sides of the domain are defined as symmetry planes. The top and outlet sides of the domain are set with pressure outflow boundary conditions. At the domain surface wall functions are used based on the local roughness of the cell and thermal stability classes.

4.2.2 EllipSys3D

Model description

EllipSys3D is a CFD solver designed for various wind engineering applications – e.g. atmospheric boundary layer flows, turbine rotor computations – it is a multi-block finite-volume solver of the incompressible Navier–Stokes equations in the general curvilinear coordinates. It uses collocated variable arrangement, employing the revised Rhie–Chow interpolation technique in order to avoid the odd–even pressure coupling. The pressure–velocity coupling in the present study was based on the SIMPLE algorithm. Furthermore, the code is designed based on a non-overlapping domain decomposition technique, which combined with its MPI parallelization enables it to highly efficiently run on distributed and/or shared-memory high-performance computation (HPC) systems.

The standard and modified model constants according to Table 2 are used in the EllipSys3D set-up in the present work. The geostrophic wind chosen is G=13 m s⁻¹, giving a maximum length scale of $l_{max}=\mathrm{28.71}$ m. A 1-D precursor computation has been conducted in order to obtain the suitable inlet profiles applied at all inlet boundaries.

To model the effects of surface roughness on the mean flow and avoid resolving the laminar sub-layer, wall functions as boundary conditions at wall surface boundaries are typically applied. In EllipSys3D, the wall boundary is placed on the top of the roughness elements; this allows large near-surface velocity gradients to be resolved using shallow (high-aspect-ratio) computational cells. The wall shear stress is accordingly used to specify the wall boundary conditions for momentum and ε equations, while a Neumann boundary condition is set for k; for a detailed description see Cavar et al. (2016). A uniform roughness of 0.1 m is applied on the entire wall surface. The laser scan map provided for the present benchmark extending over a 52.5 km × 52.5 km zone centred at the Ryningsnäs site location, in principle covering the whole wall surface area in the present study, was also fully incorporated into the EllipSys3D computations.

Numerical set-up

The computational domain is a circular grid with a radius of 17 km, centred at the Ryningsnäs metmast location. The inner zone surrounding the site has a quadratic form. It is based on equally spaced grid points and covers a zone of 5 km × 5 km. The inner zone domain fully resolves the underlying topography, while the topography in the outer (buffer) zone is gradually smoothed towards the outer boundary. The same computational grid is used for all three investigated cases (flow directions); only the inflow and outflow boundaries on the grid circumference were adjusted for each single run accordingly. Two grid sizes are considered, one using 512 × 512 grid points in the inner zone and 128 points in the outer (buffer) zone and the other one using 128 × 128 grid points in the inner zone and 64 points in the outer zone. The 3-D grid was constructed by using an EllipSys3D default hyperbolic grid generator. 192 points were used in the vertical direction, with the first cell located at 1 cm above the terrain. The vertical hyperbolic mesh growth was controlled, so the zone up to a 50 m height had cells not stretching higher than 1 m. The top boundary was located at a 9 km height. The considered grids had approximately 100 million grid points (3072 blocks of 32³) in the 10 m resolution run and approximately 9.5 million grid points (288 blocks of 32³) in the 50 m resolution run.

4.2.3 Alya

Model description

Alya is an HPC code developed at the Barcelona Supercomputing Centre (BSC) to run large-scale applications. The code was recently tested on 100 000 processors and showed a parallel efficiency above 90 % (Vázquez et al., 2016). The k–ε model has been implemented and used in Alya for the present work.

The k–ε model in Alya uses the modified model constants according to Table 2. Roughness-based wall functions are applied as boundary conditions at the ground to avoid solving the laminar sub-layer. A wall law satisfying the Monin–Obukhov equilibrium is imposed to the momentum and turbulence equations, removing a boundary layer of thickness δ_w. The wall shear stress is imposed on the momentum equation in terms of two velocity scales. Zero diffusion through the wall is imposed for the TKE, and the dissipation ε is imposed in terms of the TKE value at a distance δ_w from the ground. For a detailed description see Avila et al. (2017).

The Navier–Stokes equation (Eq. 2) and turbulence equations (Eqs. 4–5) are discretized using a stabilized finite-element method using equal interpolation for all the unknowns. The algebraical subgrid-scale method (ASGS) was used (Codina, 1998), extended for non-linear equations (Avila et al., 2015), to provide stability for transport and Coriolis-dominating terms in the momentum equation and for transport and production–dissipation terms in the turbulence equations, permitting the removal of spurious oscillations. The ASGS method also gives stability to pressure, allowing equal interpolation spaces for pressure and velocity. The velocity–pressure problem is decoupled using an Orthomin solver (Houzeaux et al., 2011) that converges to the monolithic scheme.

A robust finite-element scheme written in block-triangular form (Codina and Soto, 1999) is obtained for the k–ε equations (Eqs. 4–5). In order to avoid instabilities and numerical convergence issues, the k and ε unknowns are not allowed to drop below a predefined limit by applying a clipping. In addition, the innermost iterative loops of the k and ε equations (Eqs. 4–5) are linearized using a Newton–Raphson scheme for the quadratic terms, considering νt and P_k constants within the innermost loops.

Once the algebraical system of equations is obtained, a deflated conjugate gradient (Lohner et al., 2011) solver with a linelet preconditioner (Soto et al., 2003) is used to solve the pressure, and a generalized minimizing residual (GMRES) solver is used for the velocity and turbulence unknowns.

Numerical set-up

The Ryningsnäs problem was solved using a cylindrical mesh with a radius of 20 km. The mesh is centred on the metmast. Surrounding the metmast the mesh resolution is 10 m over a 4 km × 4 km horizontal square. Farther from the mast the horizontal mesh size grows until it reaches a 500 m horizontal element length. The vertical resolution starts with a 0.5 m element length close to the wall and 1.2 m inside the forest. The computational domain has a vertical extension of 2000 m.

The inflow boundary conditions are defined from a precursor simulation over flat and homogeneous terrain (i.e. single-column model 1-D). The obtained fields are used also as initial conditions. Zero traction is imposed over the outflow boundaries. No velocity penetration and zero tangential stress are imposed over the top boundary.

Three different geostrophic velocities were set to the three different wind directions to match the desired velocity at mast. The geostrophic velocities were set to 12.7, 13.2, and 12.7 m s⁻¹ for the wind directions 100, 240, and 290^∘, respectively.

4.2.4 CFDWind

Model description

CFDWind is a modelling framework developed at CENER on top of the open-source CFD platform OpenFOAM (Weller et al., 1998) version 2.4.0 (CFD Direct, 2015). The model is designed for the simulation of atmospheric boundary layer flows through the solution of the incompressible RANS equations in which turbulence closure is achieved by the eddy viscosity theory and a modified version of the k–ε closure scheme that Apsley and Castro (1997) described.

As only neutral atmospheric stability was considered, the flow is assumed stationary so the SIMPLE algorithm is employed to solve the pressure–velocity coupling, while second-order upwind schemes are used for the discretization of both velocity and turbulence convective terms.

The Coriolis apparent force is added explicitly to the momentum equation together with the horizontal pressure gradient that drives the system, which is derived from the hydrostatic relation for stationary cases.

The perturbations induced by forests are modelled by adding drag and source–sink terms in the momentum and turbulence-closure equations, respectively, as proposed by Sogachev and Panferov (2006). Table 2 (Modified2 k–ε) shows the drag and closure model constants employed for the simulations. Rather than being tuned, these values follow the set employed by Detering and Etling (1985), which are derived from the experiments carried out in Panofsky et al. (1977).

Despite the fact that it is expected that wind flow will be dominated by the effects of forest features near the surface, z₀-based wall functions are implemented as boundary conditions at the ground assuming wall-bounded flow. That is, the applied horizontal kinematic shear stress is set via an effective eddy viscosity ${\mathit{\nu }}_t^{\text{wall}}$ which, together with the dissipation rate and production term of the turbulent kinetic energy equation, is obtained with the local velocity scale $u_{*}^{\text{wall}}$ computed from values of velocity and turbulent kinetic energy of the cells adjacent to the ground (see Chávez-Arroyo et al., 2014, for more details).

Similar to EllipSys3D, the wall functions consider that the computational grid is placed on top of the roughness elements so that restrictions related to the height of the cells adjacent to the ground and z₀ are avoided and high-aspect-ratio cells can be used. Outlet conditions are specified at the sides and at-slip (only tangential velocity and no-gradient) conditions are prescribed at the top of the domain.

Numerical set-up

The numerical grid was created with the meshing software WindMesh. The tool has been developed jointly by BSC and CENER for the automatic and fast generation of grids over terrain. There are currently two different versions further developed by each institution: BSC-WindMesh (Gargallo-Peiró et al., 2015), which was employed in the simulations of the Alya model, and the CENER-WindMesh (Gancarski and Chávez-Arroyo, 2017) version, which was used for the generation of the grids for the CFDWind and UUCG runs.

CENER-WindMesh creates structured terrain-following grids optimizing parameters such as orthogonality and skewness by applying filters to the 2-D (ground) mesh and elliptic smoothing techniques for the final 3-D mesh. The mesh is designed so that terrain is smoothed far from the area of interest, whereas towards the central zone the cells are refined to the maximum resolution established. Only real topography is considered for the grid generation in the centre. The “transition” zone between boundaries and the central zone is a progressive blend between the real terrain and flat boundaries.

Similar to previous approaches, a precursor run is conducted prior to the full-terrain simulation (successor) in order to create the equilibrium profiles that serve as inlet conditions. Precursor simulations are conducted on flat domains with periodic boundary conditions on the sides with the top and wall treatment mentioned above. The PAD is set to a constant value of 19 m² m⁻³ with values of roughness length and forest height of 0.72 and 14 m, respectively. For the successor runs, the heterogeneous roughness values are created based on the canopy height map H using the simple relation of $z_{\mathrm{0}}=H/\mathrm{20}$.

The value of the geostrophic wind is chosen so that the velocity magnitude obtained in the simulations is approximately 10 m s⁻¹ at 100 m at the position of the mast for each of the three flow cases. The values are 14.66, 15, and 15 m s⁻¹, resulting in maximum length scales of $l_{max}=\mathrm{32.2}$, $l_{max}=\mathrm{33}$, and $l_{max}=\mathrm{33}$ m for the inflow directions 100, 240, and 290^∘, respectively.

The computational domain is square-shaped and covers an extension of 18 km × 18 km × 3 km centred on the Ryningsnäs tower. From that, only a 12 km × 12 km region considers real topography in which PAD and z₀ data are interpolated from the input canopy information. The rest of the domain is set as a flat buffer area with the same PAD, z₀, and H as the precursor simulation. For each flow case, the mesh is rotated in order to align the wind direction with the normal vector of the inlet patch at 100 m above ground. The meshes consist of 20×10⁶ cells with 60 vertical levels. The first cell height is set to 1 m and then grows with a geometric function with a constant growth rate of 1.08.

4.3 LES

4.3.1 UUCGWind

Model description

The computations by UUCG were carried out with a solver implementation based on the OpenFOAM platform, version 3.0.1. A neutrally stable wind flow is computed with LES coupled with an SGS model (Yoshizawa, 1986; Yoshizawa and Horiuti, 1985), whereby ν_SGS is estimated from the subgrid turbulence kinetic energy k_SGS, which is in turn computed from a transport equation.

It is assumed that the forest acts as a porous surface exerting a drag on the flow. This is represented in the simulation with the introduction of a source term in the LES momentum equation (Eq. 2):

$$\begin{array}{}\text{(13)}& f_{\text{D},i}=-C_{\text{D}}a\left|\langle \mathit{u}\rangle \right|\langle u_i\rangle ,\end{array}$$

where C_D is the forest drag coefficient, and a is the frontal area density (assumed here to be equal to the PAD). This approach has been successfully used in wind computations over forests with LES, e.g. by Nebenführ (2015) and Boudreault (2015). The employed value of C_D=0.2 throughout the domain is taken from the latter. While the effect of the forest in Eq. (13) is applied over the resolved part of the velocity field, the dissipative effect of TKE caused by the forest is included within the subgrid scales by adding the term

$$\begin{array}{}\text{(14)}& {\mathit{\epsilon }}_{\text{SGS}}=-{\displaystyle \frac{\mathrm{8}}{\mathrm{3}}}{C}_{\text{D}}a\left|\langle \mathit{u}\rangle \right|k\end{array}$$

to the transport equation of k_SGS.

A wall model is also used to account for the roughness of the ground, although it is expected that its influence on the wind flow will be much smaller in comparison to the forest. For this, the wall model implementation available in the OpenFOAM libraries of SOWFA (Churchfield et al., 2014) was employed. The velocity deficit due to the interaction with the ground is introduced indirectly by means of applying a surface stress. For this, the model of Schumann (1975) is used, in which the non-zero components of the stress tensor at the surface are computed as a function of the friction velocity, which in turn is calculated from the logarithmic law with a local time average for the horizontal velocity. Only the modules corresponding to the modelling of the surface stress are used from SOWFA by importing these from OpenFOAM 2.x into the version used for the simulations.

Numerical set-up

The computational domain consists of a box with the dimensions 32 km × 20 km $\times \sim \mathrm{1.2}$ km in the longitudinal, spanwise, and vertical directions (the height varies due to the differences in terrain elevation for each wind direction). The metmast is located at 20 km in the longitudinal direction, in the mid-spanwise crossing. The mesh is created using CENER-WindMesh, described in Sect. 4.2.12, producing a mesh with a varying ground elevation. In this way, the grid consists of zones in the horizontal plane: a farm zone (20 km × 12 km, metmast at 14 km) at the interior, which is then successively surrounded by a transition zone and a buffer zone. The two outermost regions can be described as rectangular edges with a width in the longitudinal and spanwise directions of 3 and 2.5 km for the transition zone and 3 and 1.5 km for the buffer zone. The arrangement of the grid in the horizontal plane is uniform in the farm and buffer regions, while the cells stretch in the transition zone, changing their size from that of the farm to that of the buffer zones. The terrain becomes flat at the buffer edges until, at the outermost boundary (with a width of 500 m), the elevation is equal in all sides with a value of 63.06 m (100^∘), 163.25 m (240^∘), and 137.76 m (290^∘). The horizontal cell resolution is 25 m for the farm and 250 m in the buffer regions. The height of the first cell at the location of the metmast is approximately 3.4 m, while the size increases vertically with a growth rate of ∼1.05. The domain height and the number of cells in the vertical direction in every case is 1.172 km and 84 cells for 100^∘, 1.305 km and 86 cells for 240^∘, and 1.267 km and 85 cells for 290^∘.

The longitudinal axis of the domain is aligned with the wind direction for each case, so the inlet is perpendicular to the inflow. All lateral boundaries are set to periodic boundary conditions. Hence, the inlet flow is recycled from the outlet. The flow is driven by a uniform pressure gradient following the procedure described by Bechmann (2006), which also comprises the introduction of Coriolis forcing (assuming a latitude of 57^∘). In this manner, the pressure gradient is calculated for the desired geostrophic wind, which is set to yield the desired wind velocity at 98 m for each case. The complete height of the ABL is simulated to avoid the parameterization of the components of the shear stress, as they become negligible at this altitude. The ground surface is set to a wall with a uniform roughness of z₀=0.03 m, while the PAD for the cells covering the tree area is extracted (by linear interpolation) from the file at 10 m × 10 m resolution using the same method as for CFDWind in Sect. 4.2.10. For each wind direction, simulations are run during about 400×10³ s to develop the flow and achieve convergence of second-order statistics. Results are obtained from values averaged during subsequent computations lasting 20×10³ s with Δt=0.296 s, yielding a maximum Courant–Friedrichs–Lewy number of CFL ≈0.6 over the whole domain.

4.3.2 PALM

Model description

PALM is a massively parallelized LES solver designed for studies of the atmospheric and oceanic boundary layer. It is an open-source code (Leibniz Universität Hannover, 2018) and has been applied to the simulation of a variety of atmospheric boundary layer studies in the past 20 years. PALM solves the filtered, incompressible, non-hydrostatic Navier–Stokes equations under the Boussinesq approximation on an equidistant Cartesian grid. The subgrid-scale turbulence is parameterized by a 1.5-order closure after Deardorff (1980) solving a prognostic equation for k_SGS. Dirichlet (no-slip) boundary conditions are prescribed at the surface. Further details on the numerics and physics of PALM can be found in Maronga et al. (2015).

The forest effect is modelled by adding a sink term to the momentum equation following Shaw and Schumann (1992) and Watanabe (2004). Furthermore, a sink term is added to the prognostic equation for the k_SGS according to Shaw and Schumann (1992) to ensure a rapid breakdown of turbulence in the canopy. A source term is added to the temperature equation, allowing a heat flux at the canopy top to be prescribed to account for the effect of incoming solar radiation. See Kanani et al. (2014) for equations and further details on the canopy model. A forest canopy can be prescribed by specifying the tree height and a vertical PAD profile. The PAD profile can be prescribed by using a beta probability density function (parameters α, β and leaf area index) or by specifying PAD values at discrete levels and doing a piecewise linear reconstruction. The latter approach has been used for the benchmark simulations. As per default, only a single PAD profile and tree height can be specified, and hence only a homogeneous forest can be simulated. Simulating a heterogeneous forest would have required significant code development, which was not feasible in the scope of the benchmark.

Figure 3 Wind speed profiles.The blue dashed line shows the average from the cups, and the red full line shows the average from the sonics. Error bars indicate the 95 % confidence level for the mean value. The various other markers indicate simulated wind speeds. Please see Table 1 for a legend of marker representations. LES models have been given black filling in order to increase readability. (a, d) Results from 100^∘, (b, e) from 240^∘, and (c, f) from 290^∘. (d–f) are zoomed in and in plotted in log space to improve readability.

Numerical set-up

The benchmark simulations use a model domain of 2304 m × 1152 m × 1867 m with a grid spacing of 4 m. A homogeneous forest canopy is prescribed by setting periodic horizontal boundary conditions and using averaged PAD profiles in which each of the three sectors has been averaged over the innermost 2 km. A conventionally neutral atmospheric boundary is simulated by prescribing a constant potential temperature up to a height of 500 m capped by a stable layer with a gradient of 1 K per 100 m. Coriolis force is considered assuming a latitude of 57^∘ N. The roughness length is set to 0.1 m. A geostrophic wind speed of u_g=13.0 m s⁻¹ and $v_{\text{g}}=-\mathrm{9.3}$ m s⁻¹ had to be prescribed to achieve a mean wind along the x direction at 100 m of height of about 7.4 m s⁻¹ as demanded by the benchmark specification. The simulations were run for 10 h to reach a steady state. The results have been averaged over the entire horizontal model domain and over the last 30 min of the simulation. Additionally, a vertical point profile at the centre of the domain, averaged over 2 h, is provided.

4.4 Numerical set-up overview

To summarize, four different RANS codes and two different LES codes are included in the study. Forest modelling is basically done in the same way in all codes apart from Meteodyn. All models apart from PALM use heterogeneous forest, but Meteodyn is based on a different surface data set. Domain sizes are similar apart from PALM, which uses a significantly smaller domain, but since PALM has homogeneous forest with recirculation the domain size is directly comparable. A summary of some key modelling properties is found in Table 3. The numerical challenge stretches from the use of a commercial code to state-of-the-art research codes using up to about 40 million cells, modelling of 20 000 physical seconds, and the use of approximately 20 000 CPU hours.

Gargallo-Peiró et al. (2015)Gancarski and Chávez-Arroyo (2017)Gancarski and Chávez-Arroyo (2017)

Table 3 Numerical set-up. Cell size refers to the horizontal grid size in the inner domain.

5 Results

One main purpose of RANS and LES modelling within the wind energy community is to extrapolate tower measurements vertically and spatially. In the next section, the vertical extrapolation (vertical profiles) is reported first, followed by the horizontal extrapolation (planes).

5.1 Vertical profiles

Wind speed, wind veer, and turbulence are crucial to power production. These three quantities, in the form of mean wind speed $S=\sqrt{U^{\mathrm{2}}+V^{\mathrm{2}}}$, mean wind direction profile, and turbulent kinetic energy TKE $=\mathrm{0.5}(\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{v^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}})$, are evaluated in Figs. 3, 5, and 6, respectively.

Modelled and measured profiles are shown for the three different wind directions described in Sect. 3.2. As is apparent from studying Fig. 3 most models actually show a slightly lower wind speed than the targeted 7.4 m s⁻¹ at 100 m of height. The wind profiles are also provided in logarithmic height coordinates and it is apparent that the measurements have a deeper log-linear region than most of the modelled curves.

Figure 4 Modelled and observed wind shear. The lower three plots are zoomed in on the instrument heights to increase readability. The blue dashed line shows the average from the cups, and the red full line shows the average from the sonics. Error bars indicate the 95 % confidence level for the mean value. The various other markers indicate simulated wind speeds. Please see Table 1 for the marker representations. LES models have been given black filling in order to increase readability. (a, d) Results from 100^∘, (b, e) from 240^∘, and (c, f) from 290^∘. Panels (a–c) and (d–f) show the same results, but panels (d–f) are zoomed in to improve readability.

Figure 5 Wind direction profiles. The red full line shows the average from the sonics. Error bars indicate the 95 % confidence level for the difference between the wind direction at each height and direction at 40 m. The various other markers indicate simulated wind directions. Please see Table 1 for the legend of marker representations. LES models have been given black filling in order to increase readability. (a) Results from 100^∘, (b) from 240^∘, and (c) from 290^∘.

Figure 6 TKE profiles. The red full line shows the average from the sonics. Error bars indicate the 95 % confidence level for the mean value. The various other markers indicate simulated TKE levels. Please see Table 1 for a legend of the marker representations. LES models have been given black filling in order to increase readability. (a) Results from 100^∘, (b) from 240^∘, and (c) from 290^∘.

Most models overestimate the wind speed gradient, as reported in Fig. 4. The overestimation increases with height and in the upper layers is close to a factor of 2. The fact that most models have a lower 100 m wind speed compared to the target 7.4 m s⁻¹ means that one has to be careful when interpreting how good the models are at estimating the wind speed gradient, but a 0.5 m s⁻¹ difference in 100 m wind speed roughly means a $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ s⁻¹ difference in wind speed gradient (assuming the difference is spread out over 100 m). As can be seen in Fig. 4, the difference for most models is at least twice as large and is consistently over-predicted by all RANS models except Meteodyn. Meteodyn is the only RANS not to use the PAD input, and thus it is difficult to know whether the smaller shear is due to the surface boundary condition or some other modelling aspect. In terms of shear, one can see that OpenFoam UUCG stands out as the best-performing model. The forest parameterization and PAD data were the same for that model as for the RANS models using PAD, so the difference in shear cannot be explained by a difference in the amount and placement of the surface drag elements. The main difference between OpenFoam CFDWind and OpenFoam UUCG is that the latter is run in LES mode, which seems to result in either better estimation of the boundary layer height, more realistic mixing of the velocity deficiencies, or both. The other LES model, PALM, did not run with the detailed PAD input, but did use averaged PAD profiles averaged with the innermost 2 km radius for each of the three directions. As seen in Figs. 3 and 4, PALM has lower shear compared to the RANS models, but it is unclear how much of that is due to the LES effect and how much is due to the constant PAD profile used. One interesting thing to notice is that while all the models using the detailed PAD fields, as well as the measurements, have the largest shear at 100 m in the 240^∘ sector, PALM shows the largest shear in 290^∘ at 100 m despite having much lower PAD in that direction. This is in line with analytical theory for homogeneous forests, which predicts a maximum of the roughness length at moderate PAD after which blocking effects lead to gradually lower roughness lengths with increasing plant area index (PAI) (Jackson, 1981). Most forests considered for wind energy have a substantial heterogeneity, though, and as can be seen in Fig. 1b the lower average PAD within the innermost 10 km does not come from areas with less dense forest, but rather from a large area without forest, over which the flow may be able to adopt to a lower roughness environment.

In an earlier publication (Arnqvist et al., 2015) it was shown that the wind turning with height (veer) was considerable at rotor heights, especially in stable stratification. Most models participating in the study show a veer of 1–3^∘ between 50 and 150 m (Fig. 5); this is about half of the veer found by the measurements, but the kinking of the measured curves also indicates the difficulty in measuring small deviations of the wind direction as wind load on the tower and booms as well as alignment accuracy all add to the uncertainty. The general shape of the wind direction profile in the Ekman layer, however, suggests that the models represent the relevant physics accurately apart from the Meteodyn model in which wind direction strangely does not seem to be coupled to the balance among Coriolis force, stress divergence, and pressure gradient. Another very interesting point is that all the models using detailed PAD input show a reverse in the wind direction turning centred between 20 and 60 m: the lowest in 240^∘ and the highest in 100^∘. This reverse could be an indication that the main driver of the flow is turbulence transport from above, and given the Ekman spiral aloft, the flow carries momentum directed to the right-hand side, leading to a counterclockwise turning with height within the forest. This behaviour has recently been shown to govern the direction of the wake behind a wind farm (van der Laan and Sørensen, 2017). It should be noted, however, that earlier field studies have shown a prevalence of clockwise turning with height also within the forest (Pinker and Holland, 1988; Smith et al., 1972). Unfortunately, the measurement site lacked wind measurements within the forest, and the mechanism behind the forest wind direction profile has not been further investigated in this study.

Figure 7 (a–e) Simulated wind speed at 40 m above the local ground height. (a) EllipSys3D high-resolution 1, (b) Alya, (c) OpenFoam CFDWind, (d) OpenFoam UUCG, and (e) PALM. The size of the boxes is 3 km × 3 km centred around the measurement tower (marked by an ×). (f) The forest height.

Of the three directions, 290^∘ stands out as the one with the lowest overall mean PAD, as seen in Fig. 2. The main contributor to the lower PAD average, however, is mainly low forest or lack of forest in the far upstream region, which is apparent from Fig. 1b in which an area without forest is seen at around 7–12 km from the metmast in the fetch in the 290^∘ direction. It is likely that this low roughness area that causes the TKE level to drop more quickly with height in the 290^∘ direction compared with the other directions (Fig. 6). While all models seem to get the overall level of TKE approximately right, only the LES version of OpenFoam captures the decrease in TKE with height for all three directions. The RANS models show much less variation between the different directions. Interestingly, the PALM LES shows almost no difference in TKE between the different directions even though the average PAD is different in all three directions; see Fig. 2d. It should be noted though, that the wind speed at 100 m in PALM was lower in the 240^∘ direction, which may explain why the turbulence level is not higher for that sector.

The purple lines in Fig. 6 shows the different set-ups of EllipSys3D RANS, and while the effect of resolution does not seem to influence the results much (used values are found in Table 3), the choice of turbulence-closure constants play a huge role (values found in Table 2). The use of standard k–ε values produce only about half of the TKE compared to the values tuned for atmospheric boundary layers despite the shear being virtually the same.

5.2 Horizontal planes

In order to evaluate spatial differences the modellers were instructed to submit horizontal planes surrounding the measurement tower. Planes are shown here for 40 m above the local elevation in Fig. 7 and 140 m above the local elevation in Fig. 8. Also displayed in the figures are tree height and terrain elevation. Although there is some correlation with tree height at 40 m, most of the correlation is with elevation, reflecting the fact that the height is above local terrain and all models can be seen to have wind fields with streamlines that are smoother than the terrain, which results in higher wind speed over high terrain and lower wind speed over low terrain. This feature is common for all models dealing with varying terrain.

Figure 8 (a–e) Simulated wind speed at 140 m above the local ground height. (a) EllipSys3D high-resolution 1, (b) Alya, (c) OpenFoam CFDWind, (d) OpenFoam UUCG, and (e) PALM. The size of the boxes is 3 km × 3 km centred around the measurement tower (marked by an ×). (f) The elevation height.

Although the models all show similar wind speed patterns there is a difference in the amount of wind speed streaks present and in the strength of the streaks. All the models show more intense streaks at higher height. The LES models show more tendency for streaks than the RANS models. EllipSys3D shows almost no streaks, whereas Alya and OpenFoam CFDWind have similar streak patterns as LES OpenFoam UUCG. The streaks correlate with topographical features, but there are also clear streaks in the PALM LES, which ran without topography.

6 Discussion and conclusions

The main aim of the study has been to investigate the state-of-the-art abilities of modelling groups to replicate measurements in neutral conditions over a forested site. Six modelling groups in total completed the whole process and submitted results. The RANS modellers using research codes used a fairly homogeneous approach to the model task, while the LES modellers took quite different approaches. Overall, a variety of options were used and in this section we will try to discuss some of the implications of these different choices.

6.1 Use of surface data

All models except Meteodyn and PALM use ALS input for topography and forest data. The fact that a variety of models (including LES and RANS) were able to use the ALS input was considered a success. The use of PAD data from ALS removes the uncertainty of having to guess the PAD or the roughness length and displacement height, which in practice can be a large source of uncertainty when estimating the wind resource at a potential site. The only model not to use the ALS was Meteodyn (instead using a constant PAD of roughly 0.025 m⁻¹) and that model also clearly stands out in the validation. Yet, some of the differences may also be attributed to the use of first-order turbulence closure.

An interesting observation noted by several of the participants was that the roughness was totally dominated by the drag of the forest and that the value of the ground roughness did not make a visible impact on the results. That could, on the other hand, be expected, since even though the forest characteristic is heterogeneous, the landscape as a whole can be considered forest covered.

6.2 Footprint

An initial question at the start of the study was whether the differences in the measured profiles among the three directions would be captured by the models, given the detailed surface data. In summary, the differences between the directions turned out to be small for the RANS models, but the LES model that used the detailed surface data produced profiles that resembled the measured profiles.

From Fig. 4 it is clear that the majority of models over-predict the shear and one may be led to believe that the forest representation, ALS conversion to PAD, is causing this discrepancy, but OpenFoam LES does in fact match the shear very well using the same forest data and the same C_d value. Especially interesting is the difference between OpenFoam LES and RANS, which are computed with implementations based on the same platform and the same grid generator. The reason for this discrepancy is an interesting point for further study.

The LES version of OpenFoam furthermore showed a much more pronounced difference between the inflow angles both in terms of shear and TKE; a possible explanation may be that the RANS models are over-diffusive, something also indicated by the fact that RANS models show fewer streaks in the horizontal planes.

6.3 Closure constants

One of the most striking outcomes of the study is that the k–ε closure with standard constants produces far too little TKE. The difference is attributed to the value of C_μ: a value of around 0.03 seems to give a reasonable TKE level. The conclusion of the participating modellers is that constants following Sogachev and Panferov (2006) should be used for future studies. Also worth mentioning is the point that all of the RANS models (apart from Meteodyn) show a too-high shear, in fact almost by a factor of 2 in the upper layers, and therefore they would be expected to generate higher levels of TKE than found in the measurements. This is also the case in the 290^∘ direction (Fig. 6c), but not in the 100^∘ direction for which the shear also is too high in the upper parts.

6.4 Resolution and domain size

As seen in Table 3, the two modelling attempts that had the lowest resolution were Ellipsys3D low-resolution and OpenFoam LES at 39.1 and 25 m, respectively. This, however, does not seem to affect the results in a negative way. The difference between the low- and high-resolution Ellipsys is very minor, both in comparison to the measurements and in comparison to the difference of changing the closure constants. The relatively low-resolution OpenFoam LES captured the wind profile well for all heights and all directions, and it seems to be important to accurately model the PAD and topography in the footprint. In order to accurately represent the turbulence moments the LES resolution needs to be such that the resolved part of the flow is large enough. For the OpenFoam LES the ratio of k_SGS to the total k was around 20 % at the tree tops and 10 % at 100 m, while k was rather well predicted and the spatial resolution is probably on the limit of accurately representing the individual components of k, at least in the lower part of the boundary layer.

Following the reasoning in Wyngaard (2010), the upstream domain size needs to be of the order of $zU/u_{*}$ to capture inhomogeneities affecting the wind at height z, which for the current site is around 15 km, exactly the upstream distance for the inner domain of the OpenFoam LES simulations. An interesting point for future studies is to examine how the results depend on upstream domain size for flow modelling with heterogeneous surface data.

Another computational expense in LES modelling is the integration time. This particular study was aimed at simulating a stationary case, and since the Coriolis force may introduce inertial oscillations it is important to make sure that the influence of oscillations does not impact the results. Another important conclusion of the study is that stationary, neutral conditions are very rare in the atmosphere, and hence future studies should investigate naturally occurring transient conditions such as diurnal cycles, evening transitions, and developing unstable turbulence.

The orders of magnitude difference in numerical challenge, both in the set-up and in the computational time used, should be considered in terms of the accuracy of the modelling results.

6.5 Recommendation for future studies

Many modellers expressed the difficulties involved in trying to determine the correct value of the U_g or pressure gradient in order to match the target 100 m wind speed. While it would be possible to instead normalize the results with u_* or some other appropriate quantity, the fact that the upwind topography and PAD seem crucial for good results points to the fact that normalization needs to be done with care since different wind speeds and turbulence levels would imply differences in the fetch, especially in stratified conditions. Furthermore, the use of a turbulence scaling parameter introduces substantial uncertainty to the measurements since the statistical uncertainty is much larger for second-order moments than first-order moments.

Based on the problem of adjusting forcing in order to match a target wind speed, measurement campaigns designed for flow model validation should attempt to measure the boundary conditions and forcing, such as the boundary layer height, vertical and horizontal fluxes, radiation, ground temperature, and geostrophic wind speed. Another option is to use a mesoscale model to compute the boundary forcing for the micro-scale models, but then care has to be taken that no additional uncertainty is introduced due to bias between the mesoscale model results and reality.

Future micro-scale model comparisons at complex forested sites should focus on the modelling of thermal stability effects. The radiative cooling during cold nights strongly affects the wind profile, introducing strong variation of heat flux in the forest. Heat transfer models over forested sites have been already implemented over flat terrain. However, new thermal models need to be developed and validated that account for complex terrain. It is also clear that the given thermally stratified conditions are more difficult to compare to stationary conditions, both due to model drift and to the fact that the atmosphere itself is mostly in transition. Comparing natural cycles would also remove the uncertainty of integration time since the physical time modelled would be the same for all participants. The representativeness of the results would also increase significantly if thermally stratified conditions were included since the strictly neutral conditions used in this study are rare in the atmosphere. Such a study would of course necessitate the use and development of unsteady RANS models.