Numerical study of how stable stratification affects turbulence instabilities above a forest cover: application to wind energy

Forest areas are of increasing interest for the wind energy industry. However, they induce complex flows with strong shear and high turbulence levels. Stably stratified atmospheric conditions, typical during nighttime and especially in winter, add to the challenge of accurately estimating wind resources. Such conditions typically imply strong wind shear and cause larger structural fatigue loads to wind turbines. In this work, large-eddy simulations are performed in neutral and stable conditions over a forest to analyze the influence of the combined effect of forest and thermal stabilities on the unsteady characteristics of the wind flow. Taking advantage of the unsteady resolution provided by the simulations, turbulent characteristics of each thermal stability including the organization of turbulent structures are presented. The resulting comparison between the two cases is put into perspective for wind energy applications.


Introduction
The rapid growth of wind energy during recent years has made it economically feasible to build wind turbines in forested regions. Such locations have a huge benefit from the social acceptance point of view, as many people are sensitive to heavy constructions in their neighborhoods. However, the wind profile over forested terrain is characterized by an increased level of turbulence which may be responsible for a significant degradation of the power curve for the same wind speed at hub height as shown by Clifton et al. [1]. Wind turbines installed within forested areas are also expected to have high aerodynamic loads due to strong vertical wind shear and the associated turbulence production. This would eventually affect their life-time and maintenance costs. In addition, forest covers induce specific turbulent organization unlike a flat terrain that is likely to affect power production [2].
Stably stratified atmosphere and thermal surface inversions are typical of nighttimes and wintertimes and must be accounted in inland wind energy production, especially in high-latitudes locations. Stable boundary layer (SBL) is usually associated with low wind speeds at turbine hub height as well as strong wind shear. These effects induce high deformation strain on the blade as well as greater structural fatigue loads [3]. Another characteristic of stably stratified flows is the low-level jet (LLJ) created in the upper part of SBL. This jet may lead to a high energy potential but can be harmful for wind turbines. Kelley [4] showed that a significant amount of wind turbine damage occurs in weakly stable flow.
These observations emphasize the need for a more detailed understanding of the flow in the lower part of the Atmospheric Boundary Layer (ABL) with thermally stable conditions above a forest canopy. This situation is typical of Northern European countries during winter. A number of studies [5,2,6,7,8,9] exists on the flow over, mainly homogeneous, forest under thermally neutral stratification. However, from a numerical point of view, very little literature [10] exists coupling forest and stable conditions. Therefore, more research is needed to improve our understanding of the flow characteristics and the turbulence organization in the lower part of the SBL, as this regime is highly relevant to wind energy.
This paper proposes to study the case of neutral and stable stratifications over a forest cover. Large-Eddy Simulations (LES) are carried out to understand how stable conditions affect the turbulent characteristics of the wind over a forest cover. In the first part, a LES is performed in neutral conditions and validated with wind tunnel experiments. Then, a simulation is performed with thermally stable conditions representative of typical situations and compared to the neutral case. The result comparison is shown into perspective for wind energy applications. Finally, taking the advantage of the unsteady resolution provided by the LES calculations, turbulent characteristics of each thermal stability are presented and their implications for the wind energy sector are discussed.

Numerical set-up
Following the work of Shaw and Schumann [5], a drag-force approach is used to account the presence of forest canopy into the numerical simulations. The total drag F d i in the x i direction due to forest is parameterized as where a is the leaf biomass area per unit volume, also known as a Leaf Area Density (LAD), c d is the leaf-level canopy drag coefficient taken here 0.15 [5], u i (u 1 = u, u 2 = v, u 3 = w) is the instantaneous velocity components in the x i (x 1 = x, x 2 = y, x 3 = z) direction, and S is the instantaneous velocity magnitude. Using the Reynolds decomposition, instantaneous velocity components u i can be expressed as is the time-averaged velocity and u i the fluctuating part in the respective directions. In this work, the incompressible filtered Navier-Stokes equations are solved including the Coriolis forces and using the Boussinesq approximation for buoyancy forces along with the continuity equation. In addition, a potential-temperature transport equation is solved to take into account the atmospheric stability. Here, the LES calculations are carried out using the Simulator for On/Offshore Wind Farm Applications (SOWFA) package [11]. SOWFA is a collection of OpenFOAM-based incompressible atmospheric/wind-farm LES solvers as well as boundary conditions for ABL flow problems especially in wind energy. Using SOWFA, Churchfield et al. [11] noticed that the standard Smagorinsky sub-grid scale model performs much better than, for example, the Lagrangian-averaged scale-independent (LASI) dynamic Smagorinsky model. Indeed, the standard model with a constant C s =[0.1 0.17] gives similar results to the GABLS inter-comparison case [12,13] for SBL flows, while the results from the LASI dynamic model are poorer. Therefore, following [11], the standard Smagorinsky model with C s = 0.13 is used in this study. As SOWFA does not provide forest related treatment, a canopy-drag model (Eq. 1) is implemented into the LES solver during this work.
For the purpose of validation of the neutral case, the numerical simulations use the forest LAD profile from the recent work by Conan et al. [9]. In this work, wind-tunnel measurements are performed on a 1/400 scale forest canopy model under neutrally stratified conditions and time series are analyzed to study high order turbulence statistics, including the contribution of coherent structures to momentum flux. To consider the flow over an infinitely long forest, the forest canopy is specified as horizontally homogeneous but vertically heterogeneous using the LAD profile used in the wind tunnel study [9]. The Leaf Area Index (LAI), which is the vertical integration of the LAD, is 10.5 and is more than twice the LAIs utilized in the previous studies [8,5,2,6,7,10]. Therefore, it represents a denser forest cover.
In this study, the wind-flow simulations over a homogeneous forest cover are performed using a computational domain of 80H ×40H ×20H in longitudinal (x), transversal (y) and vertical (z) directions, respectively. Here, H = 20 m is the height of the forest. The domain is discretized into 320 × 160 × 126 grid-cells in x, y and z directions respectively. In this way, the grid resolution in both horizontal (longitudinal and transversal) directions is fixed to 5 m. In the vertical direction, the grid resolution is set to 2 m in the forest (z < 20 m), and is then gradually decreased to reach 5 m at the top boundary. Similar grid-resolutions have been used in previous LES studies over forest canopies [8,10]. Periodic boundary conditions are used in both horizontal directions. The flow is driven by fixing the pressure gradient force which adjusts the geostrophic wind to U g = (6.13, 0, 0) m·s −1 at the height of 400 m. The velocity is initialized to the geostrophic wind throughout the domain. The Monin-Obukhov scaling laws are used at the lower boundary. The surface roughness-length z 0 on the lower boundary is set to 0.05 m at full scale, however, this value has very limited influence on the velocity profile over forest.
In the stable case, the potential temperature θ is initialized to raise from 297 K at the ground to 300 K at the top boundary. This leads to a lapse rate of 7.5 K·km −1 . Below 100 m, random fluctuations with amplitudes of 0.1 K and 0.5 m·s −1 are added to the initial potential temperature and horizontal-velocity fields, respectively.
Following Nebenführ and Davidson [10], the stable case is based on the observations in the Ryningsnäs field measurement [14], with a heat flux of Q H = −0.005 K·m·s −1 at the top of the canopy. The latitude is set to 57 • that corresponds to the Ryningsnäs field measurement site. The stable case can be validated using previous work [10,14] (same heat flux and latitude), whereas the wind-tunnel measurements [9] are strictly limited to the neutral-case validation. However, the forest LAD is the same in both cases and based on the wind-tunnel experiment.
In all simulations, the time step is fixed to 0.1 s. This ensures that the maximum Courant number stays below 0.25 in the entire domain. It was anticipated from literature [10] that the flow would develop faster in neutral conditions compared to stable. For this reason, the flow development time was set to 15 hours in the neutral case and to 20 hours in the stable case. Simulations are run for 50 hours of flow-physical time after the flow development. In addition to flow statistics, instantaneous wind and potential temperature fields are recorded for 50 hours of flow time at every height (cell-center) of a vertical line in the middle of the domain.

Validation of the simulations
This subsection discusses the validation of the LES results for both conditions. Following [10], LES results presented in Sections 3.1 and 3.2 are transfered to a local coordinate system in order to account for the wind turning effects induced by the Coriolis forces. Using the transformation, the x -axis of the local coordinate is aligned with the local wind direction at every height. The z-axis is kept vertical. Results in the local coordinate system are denoted with subscript " ".
Further, the results are normalized with the frictional velocity u * = [(u w ) 2 + (v w ) 2 ] 1/4 , where (u w ) and (v w ) are the Reynolds shear stresses. Here, the overline " · " symbol indicates the time-averaged quantities. The notation u * 2H indicates that the frictional velocity is calculated at z = 2H, otherwise the maximum value of u * (typically just above the canopy) is used in the normalization. The time-averaged LES results presented in Section 3 are also space averaged over the horizontal directions. Figure 1 compares the normalized LES results of the horizontal mean wind-speed U =  , and the resolved Reynolds shear stress (c) obtained from the neutral condition and compared with the wind-tunnel measurements (squares) [9].

and the
Reynolds shear stress u w from the neutral condition against the wind-tunnel measurements [9]. From the figure, it can be observed that the simulation successfully reproduces the mean velocity and the Reynolds stress over the forest canopy. However, the predicted velocity fluctuations in the transversal and vertical directions are underestimated, for example, by approximately 11% and 20% at the peaks respectively, except at the highest measurement positions. This might be due to the high turbulence contribution into the longitudinal direction predicted by LES (Figure 1b). In the simulation, = 0.68 and = 0.52 at z ≈ 2H. While in the windtunnel data, the corresponding values are 0.80 and 0.61, respectively. The LES predictions of the horizontal fluctuations, i.e. σ (x, ) and σ (y, ) , agree well with the wind-tunnel data inside the canopy as well as on the edge of it. The longitudinal fluctuations show higher magnitudes but in the same trend as in the wind-tunnel data. These LES results also agree with another wind-tunnel data for neutral condition by Brunet et al. [15] (the comparison is not shown here).
The stable condition in the simulation is estimated by the stability parameter H/L, with L being the Monin-Obukhov length, defined at the canopy-top as where, κ = 0.4 is the von Kármán constant, g is the acceleration due to gravity, and θ 0 is the reference temperature. At z = H, we have H/L ≈ 0.0587, which indicates stably stratified atmospheric condition over forest according to [10,16]. Figure 2 shows the normalized LES results of the mean and turbulence quantities from the stable case and compares them to the reference data [10,14]. It can be observed that the mean flow obtained by our LES agrees greatly with the reference data despite having completely different forest properties, such as LADs or LAIs. The mean velocity profile shows typical features of a stable boundary layer, such as proper representation of the LLJ at 10 < z/H < 18 (Figure 2a). In the reference LES, the LLJ appears weaker, this may be owing to their initial constant potential temperature profile without capping inversion imposed in [10]. (σ x,y,z ) ℓ /u * 2H z/H Figure 2. Vertical profiles of the LES results (solid lines) of the horizontal mean wind-speed (a) and a zoom (b), the RMS of the fluctuations (c), the resolved Reynolds shear stress (d), and the turbulent kinetic energy (e) obtained from the stable condition and compared with the reference LES results (dash-dotted lines) from [10] and the field measurements (circles) [14].
For the fluctuations (σ i, in Figure 2c) some deviations are observed between the two LES cases. However, our LES results at z = 2H show better agreement with the field data [14] compared to the other LES (Figures 2c and 2e). Above this height, the field data show stronger turbulence fluctuations, especially in the transversal direction, while both simulations show a decrease of the fluctuations with height. At z = 2H, the resolved Turbulent Kinetic Energy (TKE) is 90% of the total TKE, which indicates the simulation resolves most of the turbulent motions. More detailed validation on the higher order turbulence statistics is shown in Section 4.

Influence of thermal stratification
In Figure 3, the normalized LES results from the two stability conditions are compared. At first, it can be noticed that the normalized wind-speed increases in the stable condition (Figure 3a). The wind-speed in the stable case is higher by 34% at z = 5H (hub height) and 40% at z = 20H as compared to the neutral case. The maximum is found to be 71% at z = 13.4H due to the LLJ. The velocity increase can certainly lead to high power potential. Further, the longitudinal velocity fluctuation (σ x, ) is smaller in the stable case compared to that in the neutral case at every height above the forest (Figure 3b). The transversal (σ y, ) and vertical (σ z, ) fluctuations are found to be of similar amplitude for 1 ≤ z/H < 5. At z = 5H, the longitudinal fluctuation in the stable condition is weaker by 75% compared to that of the neutral. The overall reduction of the turbulence is clearly depicted in Figure 3(d) presenting the vertical profile of TKE. Moreover, in the stable case, the fluctuations in all three directions decrease more rapidly with increasing altitude above the forest. Inside the canopy, there is no difference for the mean velocity but a very small difference can be seen for the resolved TKE. In general, the TKE is smaller with the stable condition, for example by 71% at z = 5H, with respect to the neutral. Because of the different boundary-layer heights in the two cases, the normalized Reynolds shear stress differs between the two stability conditions as it decreases more rapidly for the stable case (see Figure 3c). The boundary-layer height in the stable case is at z/H = 15, which is about 25% smaller than in the neutral case. The lower boundary-layer height is generally responsible for the decreased level of turbulence such as TKE or turbulent fluctuations in the stable case. In summary, above the forest, SBL provides higher velocity magnitude and a lower level of turbulence, which are favorable for wind energy. (σ x,y,z ) ℓ /u *

Application to wind energy
Here, the results obtained in Sections 3.1 and 3.2 are put into a wind energy perspective. For this, a typical middle size wind turbine (≈ 2.5 MW) with diameter of 90 m and hub height of 100 m is selected. In Figure 4, the altitude is normalized with the hub height, and the light-gray background depicts the rotor diameter. The different stability conditions tested are compared in the frame of the wind turbine, the mean wind direction and speed are taken at the hub height. The LES results are also compared to the International Electrotechnical Commission (IEC) standards [18] and to the norms by VDI guidelines [17]. The IEC standards are used for the selection of the turbine by providing the shear exponent, the reference speed and the reference turbulence intensity specific of each turbine category. VDI guidelines provide empirical formulations of mean velocity and turbulence profiles over a flat terrain (without forest) for various surface roughnesses (z 0 ) under neutral conditions. The comparison of the velocity profiles between the VDI norms [17] (neutral and no forest) and the neutral simulation with forest in Figure 4(a) clearly shows an increased wind shear (due to the forest) affecting the entire rotor area. The shear exponent α = z hub U hub dU dz , presented in Figure 4(d), is calculated as a local value. It is very high in the presence of the forest gives the same conclusion (see Table 1). The increased shear is likely to affect the pitching angle of the turbine. Figure 4(b), shows that the stably stratified condition dramatically affects the change in wind direction. A rotation of 10 • is observed through the rotor in stable conditions whereas only 3 • are observed in neutral conditions. Such changes in wind direction (due to the prevailing of the Coriolis forces compared to the inertial forces and turbulent mixing) are expected to affect power production and to increase the yawing moment with possibly harmful consequences for turbine lifetime.
The global turbulence intensity ( Figure 4(c). As expected, under the neutrally stratified conditions, the forest creates a dramatic increase in the turbulence level and maintains it throughout the entire rotor as compared to the norms [17] without forest. This leads to turbulence levels in the order of 18%. The maximum is reach at the lowest position of the rotor (19%). In the stable case, the negative thermal flux has a tendency to reduce the turbulent level and therefore it damps the turbulence generated by the forest. As a consequence, the maximum turbulent level seen by the rotor is 12% at its lowest location, and decreases rapidly with hight down to 11% at the hub height and 8% at the highest rotor location. This indicates that stable conditions should act favorably in forest areas as the turbulent level is greatly reduced. Table 1 summarizes the results. Compared to the IEC values, the recommendation for the maximum turbulence intensity is not fulfilled in the neutral case. However, the maximum turbulence intensity under the stable conditions matches the IEC recommendations.
It has to be noted that the values in the IEC recommendations are taken only at hub height. However, they may not be representative of the situation along the entire rotor diameter. For example, the shear exponent decreases from 0.8 to 0.3 in the neutral case and from 0.85 to 0.5 in the stable case along the rotor diameter. Similarly, the turbulence intensity in the stable case ranges from 12% to 8%. This issue is important as it may lead to different categories of turbulent characteristics [18].

Analysis of turbulent instabilities
In this section, the focus is placed on the analysis of turbulent fluctuations in the neutral and stable conditions. The study is based on 50-hour time series recorded at 5 Hz after the convergence time of 20 hours. A quadrant-hole analysis (Q-H) is performed in order to understand the organization of turbulent structures and to elucidate dominant mechanisms in the atmospheric conditions modeled here. In this method, commonly used in the atmospheric turbulence analysis [19,20],  time series of the longitudinal and vertical velocity fluctuations (respectively u and w ) are divided in four quadrants: Q 1 (u > 0 and w > 0), Q 2 (u < 0 and w > 0), Q 3 (u < 0 and w < 0), and Q 4 (u > 0 and w < 0). Q 2 events are called ejections and Q 4 sweeps. A filter is applied to account only for the most important fluctuations: |u w | > hσ u σ w with h being the hole size. From this classification, conditional averages can be performed to calculate the contribution of each quadrant to the turbulent shear stress. Q T (= Q 1 +Q 2 +Q 3 +Q 4 ) is defined as the sum of the contribution of all quadrants.
To bring the discussion towards wind energy applications, results are plotted with height normalized by the hub height similarly to Section 3.3. To be comparable to numerical results, the altitude in the wind tunnel experiment [9] is first normalized by the roughness sub-layer height (z RSL ) and then by the hub height. z RSL is determined by the change of sign of the skewness of both velocity components as proposed by Raupach [19].

Neutral conditions and validation
Numerical results for neutral stratification are compared to experimental wind tunnel data from Conan et al. [9]. Figure 5(a) shows the vertical profiles of the velocity skewness u 3 i (thirdorder statistics) in the longitudinal and vertical directions. Similar behavior is found between the experiment and the numerical simulation, however, absolute values are higher in the wind tunnel.
The contribution of each quadrant to the turbulent shear stress u w is calculated using h = 1. It is found that the sum of all quadrants Q T ) represents only 20% of the total time (meaning that 80% of |u w | is below σ u σ w ) but still 80% of the turbulent shear stress. Q 2 and Q 4 (ejections and sweeps, respectively) are found to be of utmost importance compared to quadrant Q 1 and Q 3 which have a negligible contribution (5%). From Figure 5(b), it is observed that the contributions of all quadrants Q T from the numerical modeling are in very good agreement with the wind tunnel data. However, the relative contribution of sweeps and ejections is different between the experiment and the numerical simulation for the neutrally stratified flow. For both cases, sweeps (red) dominate the near-canopy region and ejections (blue) dominate above it, but it is clear from Figure 5(c) that the dominance of one quadrant to the other is weaker in the numerical simulation (ratio Q 2 Q 4 closer to 1) than in the wind-tunnel data. This may be due to a stronger generation of turbulence in the transversal and vertical directions during the experiment (Figure 1). The time series used in the Q-H analysis account only the resolved part of the LES, some structures may have been missed due to the LES spatial filtering.
Relatively to the mean flow, turbulent structures carrying most of the turbulent shear are high speed events entering inside the canopy in the vicinity of the forest (Q 4 ), whereas, at a higher altitude, they are mainly rising low-speed events (Q 2 ). This organization has been well described by a number of authors [2,20] and is well reproduced in the present study. This result is also consistent with the work of Raupach [19] showing the rising influence of near-wall sweeps when the surface roughness is increased.
In the context of wind energy, such organization of turbulent events is likely to increase the pitching moment applied to the rotor (see black arrows in Figure 5(c), representing turbulent fluctuations in a simple way). The unsteady nature of such events may induce instantaneous load peaks on wind turbine. The Q-H analysis performed can be related to the presence of turbulent coherent structures [21] such that the simultaneous action of an ejection at the top of a rotor and a sweep at its bottom is not to be excluded. However, the understanding of the time and space organization of these turbulent structures would require the detection of turbulent coherent structures. Globally, turbulent instabilities act on pitching moment in the opposite direction compared to the mean wind shear (see Figure 4).

Effect of stratification
An analysis of turbulent structures is then performed on the stable stratification case compared to the neutral case. The results show that the stable conditions do not affect the general turbulent organization described in the neutral case (Section 4.1). Ejections and sweeps represent 20% of the time and 80% of the turbulent shear stress (Figure 5b). However, the repartition between ejections and sweeps has changed in the sense of lowering the dominance of one compared to the other. This is clearly seen in Figure 5(c). The stable condition seems to damp the dominance of sweeps near the forest canopy and of ejections above it. This is in line with the lower level of turbulence observed in Section 3.2.

Conclusions
In this study, LES are performed to model the ABL flows above a forest canopy under neutrally and stably stratified atmospheric conditions. The numerical results are compared to the existing literature, a wind tunnel study for the neutral case, and observations and another LES study for the stable case. For this, a drag model was implemented in the OpenFOAM-based SOWFA package. In both cases, the mean velocity and turbulence results are found to be in agreement with the literature.
For the hub height chosen, the stable condition is found to provide higher (34%) wind-speed and lower turbulence level (75%) at hub height compared to the neutral case. Thus for a lonestanding turbine, stably stratified flows are favorable to wind energy as they damp turbulence coming from forest as well as increase wind speed. With respect to a flat terrain, the shear and turbulence level at hub height are increased by the presence of the forest in both conditions ( Figure 4). In the neutral case, a high turbulent level (18%) is observed up to the rotor height, which in practice is harmful to wind-turbine installation. The stable stratification is found to amplify the shear factor (by 56% at hub height compared to the neutral condition) and to increase the angle deviation. These two parameters are prejudicial to wind turbines. In both conditions, due to the forest, the shear factor is out of the IEC limit [18]. However, the stable conditions damp the intensity of the turbulence cause by the forest in an important manner, it decreases from 18% (in neutral) to 11% (in stable) at the hub height. Therefore, neutral stratification is found to be the worst case for turbulence intensity, but stable stratification the worst case for shear factor. A recent field experimental study [22] shows that turbulence intensity is the most influencing parameter to power coefficient compared to wind shear and angle. This indicates that stable conditions are somewhat favorable to wind energy production in forested regions compared to neutral conditions. It is stressed that the IEC evaluation accounts only for hub height values whereas Figure 4 shows that there can be very large variations even within rotor diameter.
The unsteady Q-H analysis shows a global turbulent organization made of ejections at the rotor top and sweeps at its lower position. Both turbulent events act in the same direction and can potentially induce instantaneous peaks in the pitching moment. A time-space study would elucidate whether these events can appear simultaneously and lead to shear factor peaks. In thermally stable conditions, the turbulent organization is similar but the dominance of sweeps on the lower part and ejections on the upper part is less pronounced. This may lead to lower unsteady loads.
Conclusions concern the study of the wind resource without wind turbine modeling, it is therefore expected that results could apply to a stand-alone turbine. In the case of a wind farm with wake interaction, the effect of thermal stability coupled with forest cover and Coriolis effect may lead to different conclusions. For example, the low level of turbulent mixing in stable conditions is known to lead to a poorer wind farm efficiency [23] but the forest on the ground enhances turbulence and may be beneficial to the wind farm.