Modelling Uncertainty in t-RANS Simulations of Thermally Stratified Forest Canopy Flows for Wind Energy Studies

Abstract

The flow over densely forested terrain under neutral and non-neutral conditions is considered using commercially available computational fluid dynamics (CFD) software. Results are validated against data from a site in Northeastern France. It is shown that the effects of both neutral and stable atmospheric stratifications can be modelled numerically using state of the art methodologies whilst unstable stratifications will require further consideration. The sensitivity of the numerical model to parameters such as canopy height and canopy density is assessed and it is shown that atmospheric stability is the prevailing source of modelling uncertainty for the study.

1. Introduction

The motivation of this paper is to assess the use of computational fluid dynamics (CFD) modelling to consider forest canopy flows where there is uncertainty in the canopy density and the level of atmospheric stability. The accuracy of the model predictions within levels of uncertainty of these two parameters is assessed in order to highlight areas where further validation data and research are required.

The computational power required to run full CFD simulations on the scale of a typical wind farm is now accessible and as a result, CFD is beginning to see greater adoption by industry for the purposes of wind resource assessment [1]. Following this trend, research activities have increased into the flow dynamics generated by non-trivial terrain and atmospheric features in order to fully realise the capabilities of CFD to describe the atmospheric boundary layer (ABL) and to meet the demanded uncertainty standards.

One element of terrain complexity which has been found to significantly increase flow modelling uncertainty is the presence of forestry. It was shown in [2] that forestry increases modelling uncertainty in terms of root mean square error by a factor of 4–5 when modelling the flow between meteorological mast pairs using a variety of industry standard modelling software packages. In [3] it was suggested that one reason for these elevated levels of uncertainty may be the regular occurrence of non-neutral atmospheric stability events in forested terrain. The buoyancy forces associated with non-neutral events are generally neglected in industry standard modelling software packages. However, they have been shown to have a significant impact on how the wind interacts with obstacles such as forestry [4,5,6].

In [7] the possibility of including the joint effects of atmospheric stability and forest canopy drag within a CFD domain was examined through the use of validation data from stratified ABL wind tunnel experiments. Whilst the results achieved in [7] were promising, the analysis was limited by a lack of availability of experimental data for an unstably stratified ABL and also possible Reynolds number scaling problems when using architectural model trees to represent a forest canopy.

For this paper, non-neutral Reynolds Averaged Navier Stokes (RANS) CFD simulations have been validated against field data from a heavily forested site in Northeastern France. Firstly, sets of stable, neutral and unstable events are identified. The neutral events are then numerically modelled in order to identify the appropriate terrain, canopy, mesh and atmospheric configurations to successfully model flow over the site. The effects of atmospheric stability are then introduced in an attempt to replicate the non-neutral events observed in the dataset.

All CFD simulations in this paper have been configured using the WindModeller (WM) software (ANSYS UK Ltd., Abingdon, Oxfordshire, UK) package which is a front end for the ANSYS CFX flow solver (ANSYS UK Ltd., Abingdon, Oxfordshire, UK). WM has been specifically designed to meet the needs of the wind energy industry and it includes the ability to simulate the effects of non-neutral stability.

The novelty of this work lies in the use of WindModeller, a reasonable approximation of the industrial state of the art in flow modelling software, to consider the extremely complex flows real world flows generated by the combined effect of thermal stratification and canopy drag.

2. Validation Data

This study uses data from a meteorological mast near Vaudeville-le-Haut which is located adjacent to a wind farm in Northeastern France (46°26′58″ N, 05°35′02″ E). There is an extensive mixed forest located to the west at a distance of c. 170 m, as shown in Figure 1.

An Institut National de l’Information Géographique et Forestière (IGN) map of the area under consideration is given in Figure 2. Four operating turbines are marked on this map; the two closest turbines to the meteorological mast are located at a bearing of 85° and a distance of 400 m and 25° at a distance of 600 m.

Data were provided between 1 January 2010 and 31 December 2011 as 10 min averages from a series of sonic anemometers (METEK USA-1, METEK, Elmshorn, Germany), temperature sensors and wind vanes on a 100 m meteorological mast as summarised in

Table 1. Meteorological sensors present on the Vaudeville meteorological mast. Instrumentation model numbers are given in parenthesis where available. Sonic anemometer were orientated into the prevailing wind from the south west and thus were not affected by tower shadow for the director sector considered, shown in Figure 3.

Height (m)	Sensor 1	Sensor 2	Sensor 3
80	Temperature sensor (PT 100, SKS Sensors, Vantaa, Finland)	3D Sonic anemometer (Metek USA-1)	Cup Anemometer (Thies First class, Thies, Göttingen, Germany)
70	Wind vane (Thies compact)	-	-
60	Temperature sensor (PT 100)	3D Sonic anemometer (Metek USA-1)	-
40	Temperature sensor (PT 100)	3D Sonic anemometer (Metek USA-1)	-
10	Temperature sensor (PT 100)	3D Sonic anemometer (Metek USA-1)	-
3	Temperature & Humidity (CS215)	Pyranometer (CMP6, Kipp & Zonen, Delft, The Neterlands)	-
1	Pluviometer	-	-
-1	Temperature sensor (PT 100)	-	-

. Solar irradiance data were provided for the same period from a pyranometer on site. The wind turbines on site were in operation during the measurement period.

Access to the full 3D sonic datasets were not available with only 10 min mean wind speed and standard deviation of wind speed provided for this research, thus it was not possible to calculate the Obukhov Length directly. In order to isolate non-neutral data with which to validate the CFD simulations, the steps described in Section 2.1 were taken. This methodology was previously applied to four sites, including Vaudeville, to isolate non neutral events and was found to provide an accurate demarcation of stability class when compared with more conventional measures of stability such as the Obukhov Length and the Richardson number [3].

2.1. Wind Speed Data

The 250–260° direction sector was examined in order to limit variations in the proximity of the meteorological mast to the forest edge. This range can be seen in Figure 3.

The effect of the forest canopy on the wind resource will vary seasonally and annually as the trees grow and develop. Such variations in the data will complicate the validation process, thus it was deemed necessary to focus analysis on data relating to a single season. The maximum possible number of observations were required for the selected season in order to provide sufficient data for validation. Also, as the analysis in [3] showed that unstable events are the least common in the Vaudeville site, a season was selected in which high irradiance levels were recorded in order that sufficient validation data would be available for all three stability classes. A summary of the available data is given in Figure 4.

Months 7 and 8 were selected as highlighted by the yellow shading in Figure 4. These data relate to July and August 2010 and allow the analysis to avoid complications due to seasonal variance in canopy density whilst providing an adequate spread of irradiance values and number of observations.

The next step was to apply the methodology outlined in [3] to identify non-neutral events. Thus turbulence intensity (TI) at 80 m and wind shear between 40 m and 80 m were calculated using Equations (1) and (2), respectively:

$$\mathrm{TI}=\frac{{\sigma }_u}{\overline{U}}$$

(1)

$$\alpha =\frac{\ln \left(U_{80}/U_{40}\right)}{\ln \left(80/40\right)}$$

(2)

As can be seen in Figure 5, the values of the observed wind shear and turbulence intensity become less sensitive to solar irradiance levels at higher wind speeds. Following [3] we assume that the narrower range of values of wind shear and turbulence intensity at higher wind speeds are characteristic of neutral stratification for the 250–260° direction sector.

The estimated neutral threshold values for the considered data are indicated as red lines in Figure 5. These are 0.15–0.28 for turbulence intensity and 0.32–0.52 for wind shear. These thresholds are then applied to the selected data set in order to identify stable, neutral and unstable events as shown in Figure 6.

The stability demarcation displayed in Figure 6 is used for qualitative purposes in this paper to assess the performance of the CFD model. The quantitative results achieved in determining stability class using this method are discussed in [3] where it was shown that up to 90% agreement was achieved when compared with demarcation achieved using direct measures of the Obukhov Length.

As can be seen in Figure 6, there are a limited number of observations which display both wind shear and turbulence intensity values which would be indicative of unstable stratification for the given site. This is despite the fact that the selected analysis period is one in which high levels of irradiance were observed, and is a limitation of the Vaudeville dataset. Regardless of this statistical constraint, the effects of stability can be clearly seen in the sample profiles presented in Figure 7. These sample profiles relate to the oversized data points in Figure 6. The time and date at which each event was measured are provided in

Table 2. Time and date at which each of the profiles in Figure 7 were recorded.

Stability Class	Time & Date
Stable	19:40 13 July 2010
Neutral	23:40 17 August 2010
Unstable	12:00 10 August 2010

.

2.2. Canopy Height Data

Data were provided for a 5 km radius around the Vaudeville meteorological mast by Intermap Technologies Ltd. (Denver, CO, USA). These data were measured using Interferometric Synthetic Aperture Radar which combines aerial, satellite and ground measurements to gather x, y, z coordinates for the ground surface surveyed. These data are then analysed using a canopy height model to derive vegetation height. The resolution of the supplied data is 5 m with an approximate accuracy in terms of measured canopy height of 2 m [8] Information on the canopy measurement techniques can be found in [9]. The distribution of canopy heights in the examined region is shown in Figure 8.

The mean canopy height in Figure 8 is 10.7 m with a standard deviation of 5.62 m. Unfortunately, no data relating to the density of the forest canopy and variation of this parameter with height were available. Thus a constant canopy density is assumed and this parameter is tuned during the neutral simulations (Section 4) to identify the appropriate value for use in the non-neutral simulations (Section 5 and Section 6). The benefit of using a constant rather than a variable canopy density profile for CFD simulations where accurate site canopy density data are not available was discussed in [10].

3. CFD Modelling

As stated previously, all CFD simulations in this paper were configured using the WM front-end to the CFX software. WM solves the Navier-Stokes equations (mass and momentum conservation) in a RANS mode. Following the analysis in [3] the Shear Stress Transport (SST) turbulence closure [11] was used for all simulations. The effects of atmospheric stability are accounted for by solving an additional transport equation for the potential temperature $\theta$, and by including stability effects in the vertical momentum equation (term $F_{B,i}$) and in the turbulence model (buoyancy turbulence production $P_{kB}$, defined below). The model also has the option to include the effect of the Coriolis force, implemented as a difference to the geostrophic balance, to capture effects associated with the development of an Ekman spiral in the boundary layer. The effect of the forestry on the flow is modelled via a quadratic resistance term in the momentum equations, as well as sources and sinks in the turbulence model to account for turbulence production and length scale redistribution. The forestry drag sources and sinks are applied to all control volumes which are identified to be located below the top of the forest canopy. Specific details of the configuration used are given below.

3.1. Model Equations

In all simulations the flow is treated as incompressible. The effect associated with non-constant density is modelled in the buoyancy force using the Boussinesq approximation. This effect is accounted for in the vertical velocity equation and in the turbulence model. The model solves the following equations:

Continuity:

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho U_i\right)=0$$

(3)

Momentum:

$$\frac{\partial }{\partial t}\left(\rho U_i\right)+\frac{\partial }{\partial x_j}\left(\rho U_j{U}_i\right)=-\frac{\partial }{\partial x_i}p+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_T}{{\sigma }_{}}\right)\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)\right]+F_{B,i}+F_{Cor,i}+F_{D,i}$$

(4)

with the body forces:

$$F_{B,i}=g\beta {\rho }_{ref}\left(\theta -{\theta }_{ref}\right){\delta }_{i3}\mathrm{Buoyancy}$$

(5)

$$F_{Cor,i}=\rho f\left[\left(U_i-U_{i,geo}\right){\delta }_{i1}-\left(U_i-U_{i,geo}\right){\delta }_{i2}\right]\mathrm{Coriolis}$$

(6)

$$F_{D,i}=-\frac{1}{2}\times \rho C_{d}A\left(z\right)\left|U\right|U_i\mathrm{Forestry}\mathrm{drag}$$

(7)

Energy conservation equation via a transport equation for the potential temperature $\theta$:

$$\frac{\partial }{\partial t}\left(\rho \theta \right)+\frac{\partial }{\partial x_j}\left(\rho U_j\theta \right)=\frac{\partial }{\partial x_j}\left[\left(\frac{\lambda }{C_p}+\frac{{\mu }_T}{{\sigma }_{\theta }}\right)\left(\frac{\partial \theta }{\partial x_j}\right)\right]$$

(8)

Turbulence closure is provided by the SST 2-equation turbulence model [11,12]:

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_j}\left(\rho U_{j}k\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_T}{{\sigma }_k}\right)\left(\frac{\partial k}{\partial x_j}\right)\right]+P_k+P_{kB}-\rho C_{\mu }\omega k+S_k$$

(9)

$$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_j}\left(\rho U_j\omega \right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_T}{{\sigma }_{\omega 3}}\right)\left(\frac{\partial \omega }{\partial x_j}\right)\right]+\left(1-F_1\right)2\rho \frac{1}{{\sigma }_{\omega 2}\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}+\rho \frac{{\alpha }_3}{{\mu }_T}{P}_k+P_{\omega B}-{\beta }_3\rho {\omega }^2+S_{\omega }$$

(10)

The effect of buoyancy on the turbulence kinetic energy is included via the source term $P_{kB}$:

$$P_{kB}=-\frac{{\mu }_T}{{\sigma }_{\theta }}g\beta \frac{\partial \theta }{\partial z}\mathrm{Buoyancy}\mathrm{source}\mathrm{for}k$$

(11)

For the eddy frequency equation, the effect of buoyancy is included with:

$$P_{\omega B}=\frac{\omega }{k}\left[\left({\alpha }_3+1\right)C_{3}max\left(P_{kB},0\right)-P_{kB}\right]\mathrm{Buoyancy}\mathrm{source}\mathrm{for}\omega$$

(12)

The level of turbulent mixing in the model is modelled via the eddy viscosity ${\mu }_T$, which in the SST model is calculated via:

$${\mu }_T=\rho \frac{a_{1}k}{max\left(a_1\omega ,S{F}_2\right)}$$

(13)

where the viscosity limiter is activated by the function $F_2$ near the wall only. $S$ is an invariant measure of the strain rate:

$$S=\sqrt{2S_{ij}{S}_{ij}},S_{ij}=\frac{1}{2}\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)$$

(14)

Another feature of the SST turbulence model is the use of a shear production limiter to avoid over production of turbulence kinetic energy in stagnation regions. The turbulence production term $P_k={\mu }_T{S}^2$ is implemented with the limiter:

$$P_k=min\left(P_k,C_{lim}\rho \epsilon \right)$$

(15)

with a value of 10 for $C_{lim}$.

The effect of the forestry drag on turbulence quantities has been modelled and discussed by various authors e.g., [13,14,15,16], often in the context of $k-\epsilon$ models. For such models, sources and sinks are added to the turbulence kinetic energy $k$ and turbulence dissipation $\epsilon$ equations, to model the added turbulence and redistributed length scales as follows:

For the turbulence kinetic energy equation

$$S_k=F_F\frac{1}{2}\rho C_{d}A\left(z\right)\left|U\right|[{\beta }_p{\left|U\right|}^2-{\beta }_{d}k]$$

(16)

where ${\beta }_p$ and ${\beta }_d$ are constants, the values of which are given in

Table 3. Modelling constants used for the canopy model [14].

Constant	Value
${\beta }_p$	0.17
${\beta }_d$	3.37
$C_{\epsilon 4}$	0.9
$C_{\epsilon 5}$	0.9

.

If using the $k-\epsilon$ turbulence mode the source terms of the dissipation rate equation:

$$S_{\epsilon }=F_F\frac{1}{2}\rho C_{d}A\left(z\right)\left|U\right|\epsilon \left[\frac{C_{\epsilon 4}{\beta }_p{\left|U\right|}^2}{k}-C_{\epsilon 5}{\beta }_d\right]$$

(17)

where $C_{\epsilon 4}$ and $C_{\epsilon 5}$ are constants, the values of which are also given in Table 3.

However, for the work presented in this paper the SST turbulence model is used. The required source terms are as follows:

For the turbulence frequency equation:

$$S_{\omega }=F_F\frac{1}{2}\rho C_{d}A\left(z\right)\left|U\right|\omega \left[\frac{(C_{\epsilon 4}-1){\beta }_p{\left|U\right|}^2}{k}-\left(C_{\epsilon 5}-1\right){\beta }_d\right]$$

(18)

This source term for the $\omega$ equation is derived from the generic relationship:

$$S_{\omega }=-\frac{\omega }{k}{S}_k+\frac{1}{C_{\mu }k}{S}_{\epsilon }$$

(19)

which itself results from the transformation of the $\epsilon$ equation into the equation for $\omega$, via the identity $\epsilon =C_{\mu }\omega k$.

A discussion on the formulation of these equations for a $k-\epsilon$ model can be found in [14,16]. The appropriate value for the modelling constants in the above equations has been an area of some research. For the current work, the values as recommended by [14] are used and these are summarised in Table 3.

The porosity of the canopy was defined by a loss coefficient, L_x, which is the product of the canopy drag, C_d, and the Leaf Area Density, A(z). In WM, this loss coefficient can be set to a constant value or can vary with height. As no data relating to the vertical structure of the canopy were available, a constant value was used for all simulations. The specific values used for each simulation will be given in the appropriate section. Note that the WM implementation of the drag force and turbulence source are based on a definition of a drag force including a factor ½. In [14] the factor of ½ is omitted in the definition of this parameter. As a consequence, the loss coefficient in [14] is to be interpreted as half the loss coefficient in WindModeller.

3.2. Boundary Conditions

At the ground, a no-slip boundary condition is used for the velocity, where the momentum fluxes through the ground are evaluated with a wall treatment, using the automatic wall function for rough walls developed by [17]. The roughness at the ground is implemented in terms of an equivalent sand roughness $h_s$. In the rough limit, the rough wall treatment implements:

$$U_{+}=\frac{U}{u_{\tau }}=\frac{1}{\kappa }\ln \left(y_{+}\right)+C-\frac{1}{\kappa }\ln \left(1+0.3h_{s+}\right)$$

(20)

with $\kappa =0.41$ and $C=5.2$.

When the ground is fully rough regime ($h_{s+}=\frac{h_s{u}_{*}}{\nu }>100)$, this can be approximated as:

$$\begin{gathered} U_{+}=\frac{U}{u_{\tau }}=\frac{1}{\kappa }\ln \left(y_{+}\right)-\frac{1}{\kappa }\ln \left(0.3h_{s+}\right)+C \\ =\frac{1}{\kappa }\ln \left(\frac{y_{+}}{h_{s+}}\right)-\frac{1}{\kappa }\ln \left(0.3\right)+C \\ =\frac{1}{\kappa }\ln \left(\frac{y}{h_s}\right)-\frac{1}{\kappa }\ln \left(0.3\right)+C \\ =\frac{1}{\kappa }\ln \left(\frac{y}{h_s}\right)+8.14 \end{gathered}$$

(21)

The above returns a $\frac{1}{\kappa }\ln \left(\frac{y}{z_0}\right)$ profile as a function of the aerodynamic roughness $z_0$ if the sand roughness $h_s$ is prescribed as:

$$h_s=z_0\exp \left(8.14\kappa \right)$$

(22)

In the log limit for the automatic rough wall treatment, the friction velocity $u_{*}$ is calculated as:

$$u_{*}=C_{\mu }^{1/4}\sqrt{k}$$

(23)

The momentum flux through the wall is calculated as:

$$F_U=-\rho u_{*}{u}_{\tau }=-\rho u_{*}\frac{U_1}{U_{+}}$$

(24)

where $U_1$ is the velocity just above the ground. For the turbulence model, the wall treatment for the turbulence kinetic energy is adiabatic (i.e., zero flux), while for the $\omega$ equation, an algebraic closure is imposed. In the rough case, with an assumed log limit, the wall value of $\omega$ is set with:

$${\omega }_{wall}=\frac{\rho u_{*}^2}{\mu }\frac{1}{\kappa \sqrt{C_{\mu }}}\frac{1}{y_{+}}$$

(25)

For the heat transfer at the wall, the boundary condition on the potential temperature is either adiabatic (zero flux) when modelling neutral surface stability conditions, or a ground temperature is prescribed from a temperature offset with respect to the advected neutral surface layer prescribed at the inflow. A negative temperature offset leads to the development of a stable surface condition downstream of the inflow, while a positive offset leads to unstable surface conditions. A wall treatment for the potential temperature is used to relate the ground heat flux $q_{wall}$ to the difference in potential temperature between the ground (${\theta }_w$) and the air (${\theta }_f)$ just above the ground. The wall function for this is a modification to the Kader wall treatment, to account for roughness effects as described [17]. It implements the following relationship:

$$q_{wall}=\frac{\rho C_p{u}_{*}}{{\theta }^{+}}\left({\theta }_w-{\theta }_f\right)$$

(26)

$${\theta }^{+}=2.12\ln \left(Pr{y}_{+}\right)+{(3.85P{r}^{1/3}-1.3)}^2-\Delta B_{th}$$

(27)

with:

$$\Delta B_{th}=\frac{1}{0.41}\ln \left(1+C0.3Pr{h}_{s+}\right)$$

(28)

$$Pr=\frac{\mu C_p}{\lambda }$$

(29)

$$C=0.2$$

At the inflow, profiles are prescribed for the velocity, for the turbulence kinetic energy and dissipation rate, while Neumann boundary conditions are applied to the pressure. When modelling stability, a profile for the potential temperature is also prescribed. By default, the latter assumes adiabatic conditions in the boundary layer at the inflow, capped by stable conditions above the boundary layer, with a potential temperature gradient that can be prescribed by the user (default value of 3.3 K/km as per the standard US atmosphere [18]. When non-adiabatic conditions are imposed at the ground for the temperature, the model then develops a stable or unstable surface layer which grows downstream of the inflow. The conditions applied at the inflow for the velocity and turbulence quantities depend on the selected physics. When modelling purely neutral flow, without the Coriolis force or Ekman spiral, the profiles applied at the inflow are set up with the standard Equations (30) and (31) as defined by [19]. These profiles are calculated using a default $C_{\mu }$ value of 0.09:

$$U(\tilde{z})=\frac{u_{*}}{\kappa }\ln (\frac{\tilde{z}}{z_0})$$

(30)

$$k\left(\tilde{z}\right)=\frac{u_{*}^2}{\sqrt{C_{\mu }}}$$

(31)

$$\epsilon \left(\tilde{z}\right)=\frac{u_{*}^3}{\kappa \tilde{z}}$$

(32)

where $\tilde{z}$ is the height above the ground. In terms of user input, the profiles are prescribed from a reference mean horizontal wind speed, U_ref, and the height above ground level at which it occurs Z_ref along with the surface roughness $z_0$. From these user-defined criteria, WM then calculates a value of $u_{*}$ using a form of the log law as shown in Equation (33):

$$u_{*}=\frac{\kappa \times U_{ref}}{\ln \left(\frac{Z_{ref}}{z_0}\right)}$$

(33)

When modelling stability effects, and including Coriolis, the associated Ekman spiral is prescribed following a formulation proposed by [20] This provides profiles for the horizontal wind speed components which follow the Monin-Obukhov similarity theory in the surface layer, and adapt the atmospheric length scales in the upper part of the boundary layer accounting for static stability effects above the boundary layer and effects associated with the earth rotation. The boundary layer height required to specify the velocity profile at the inlet is obtained from a multi-limit diagnostic method proposed by [21]. More details on the implementation of this approach in an earlier version of CFX can be found in [22]. For the cases simulated here, the Obukhov length was set to 10,000 m (essentially neutral surface layer), to be consistent with the assumed neutral conditions applied at the inflow for the potential temperature. For the turbulence quantities, the following profiles are imposed at the inflow, in conjunction with the Zilitinkevich et al. velocity profiles:

$$k\left(\tilde{z}\right)=\frac{u_{*}^2}{\sqrt{C_{\mu }}}{\left(1-\eta \right)}^{1.68}$$

(34)

$$\epsilon \left(\tilde{z}\right)=\frac{u_{*}^3}{\kappa \tilde{z}}\times 1.03\times \left[1+\frac{0.015}{{\tilde{z}}^{0.9}}\max (\ln \left(\frac{\tilde{z}}{z_0}\right),0)\right]\exp \left(-2.8{\eta }^2\right)$$

(35)

where $\eta =\tilde{z}/h$, and $h$ is the boundary layer height. These profiles were fitted from equilibrium profiles resulting from 1D simulations of a vertical columns with homogeneous flow conditions in the horizontal directions, obtained for adiabatic ground conditions.

At the domain outflow and at the top boundary an opening type of boundary is used. Von Neumann boundary conditions are applied to the velocity components, the turbulence quantities and the potential temperature, as long as the flow at this location is out of the domain. In case of flow entering the domain at the outflow location, a Dirichlet boundary conditions is applied to the turbulence variables and potential temperature, using the same profiles as used at the inflow. The pressure at the outflow and top boundary is prescribed with a profile balancing the hydrostatic conditions associated with the buoyancy term in the vertical velocity equation for the temperature conditions applied at the inflow. Hydrostatic equilibrium implies:

$$\frac{\partial }{\partial z}p=g\beta {\rho }_{ref}\left(\theta -{\theta }_{ref}\right)$$

(36)

For an inflow potential temperature profile given with:

$$\begin{gathered} \theta \left(z\right)={\theta }_{ref}forz<Z_{ref} \\ \theta \left(z\right)={\theta }_{ref}+\gamma \left(z-Z_{ref}\right)forz\ge Z_{ref} \end{gathered}$$

(37)

The integrated pressure profile balancing the hydrostatic is then:

$$\begin{gathered} p\left(z\right)=p_{ref}forz<Z_{ref} \\ p\left(z\right)=p_{ref}+g\beta {\rho }_{ref}\gamma \frac{{\left(z-Z_{ref}\right)}^2}{2}forz\ge Z_{ref} \end{gathered}$$

(38)

In Equations (37) and (38), the parameter $\gamma$ is the temperature lapse rate. When not modelling stability, the pressure profile at the outflow and top boundary is simply set to a constant value of 0 Pa. The initial conditions prescribed within the computational domain for all t-RANS simulations were as per the height dependent boundary conditions set at the inlet.

3.3. Domain Description

Circular computational domains are generated in WM where the outer edges are divided into 24 surfaces, as shown in Figure 9, which allows various wind directions to be considered using a single domain configuration. Twelve of the outer surfaces are used for the inflow condition, and the other twelve represent the outflow. For the present study, the radius of the domain was set to 7.5 km and the domain was centred on the meteorological mast (46°26′58″ N, 05°35′02″ E). The wind direction was set to 255° in order to coincide with the centre of the direction sector investigated, shown in Figure 3.

Topographical details from the 90 m resolution Shuttle Radar Topography Mission (SRTM) [23], dataset were used to generate a tessellated surface of triangular elements which captured the undulations in the terrain. This resolution was considered satisfactory given the domination of canopy effects and the simple terrain in the 250–260 ° direction sector. This terrain detail was limited to 5 km from the mast in all directions with the outer most 2.5 km extended radially at constant local elevation. This configuration is used to allow the wind characteristics to adjust to the applied surface roughness height before encountering the topography.

The aerodynamic surface roughness length applied in all simulations was z₀ = 0.04 m which is what would be expected for a site containing low grass [24]. Values of 0.1 m and 0.001 m were also tried, however the impact on results was negligible. This is due to the fact that much of the fetch along the 255° direction is occupied by forestry and thus the surface roughness itself will have a reduced role in dictating the wind characteristics.

The height and extent of the Vaudeville forest was described by a set of x, y, z coordinates derived from the Intermap data described in Section 2.

The height of the domain was set to 2 km for the majority of simulations. Any alterations to this will be discussed where applicable. A description of the mesh used will be given in Section 3.4.

In order to capture the additional flow detail introduced by the buoyancy effects, all WM simulations which include stability are investigated as transient RANS simulations. The overall physical simulation time is calculated using:

$$\mathrm{Overall}\mathrm{time}=\frac{2.5\times \mathrm{Domain}\mathrm{diameter}}{U_{geo}}$$

(39)

The initial time step is set to 10 s and increases to a maximum value of 30 s depending on how quickly the simulation converges.

3.4. Mesh Sensitivity

A mesh sensitivity study was conducted using a neutral configuration. A constant canopy loss coefficient of L_x = 0.05 m⁻¹ was used for all simulations along with U_ref = 6.5 m/s at Z_ref = 40 m.

The circular domain generated by WM is divided into nine zones for the purposes of meshing as shown in Figure 9. In each of these zones, a block structured hexahedral mesh is generated in accordance with user-defined criteria. This configuration allows all direction sectors to be considered using a single mesh which considerably reduces the time required to set up simulations for the purpose of a resource assessment.

In Figure 9a, the critical dimensions which define the mesh are shown. For all simulations the following values were used: the edge length of the centre block, L = 2.33 km, the radius of the inner zones R1 = 5 km and the radius of the outer zones R2 = 7.5 km. The height of the domain was set to 2 km for the mesh sensitivity study. The structure of the mesh itself is defined by setting a maximum horizontal, Hz, and vertical, Vt, mesh resolution for the centre block.

For all simulations, a 10 cell inflation layer of 2 m high cells was applied to the floor boundary throughout the domain with a vertical expansion factor of 1.15 thereafter. A horizontal expansion factor of 1.1 was used for the both the inner and outer zones. The maximum horizontal and vertical cell size within the central block was then adjusted in order to produce three different meshes; details of which are given in

Table 4. Mesh resolutions used for the mesh sensitivity analysis.

Mesh	Maximum Cell Size		Control Volumes	Nodes	CPU Time
	Hz	Vt
Coarse	100 m	100 m	87,696	93,478	5 min
Medium	20 m	50 m	2,149,056	2,215,626	60 min
Fine	10 m	25 m	13,418,460	13,638,322	480 min

. All simulations were conducted on a High Performance Computing (HPC) cluster which consists of 161 nodes, each having two six-core Intel Westmere Xeon X5650 Central Processing Units and 24 GB of memory. Each simulation was divided among twelve cores in order to avoid problems which may occur from segmenting the domain into an excessive number of parallel computations.

In order to compare the quality of the results achieved using the three levels of mesh, values for the mean horizontal wind speed, U, and turbulent kinetic energy, k, at the meteorological mast location up to a height of 200 m were determined. The results of the mesh sensitivity study are shown in Figure 10. Simulated values have been normalised to the reference velocity U_ref = 6.5 m/s.

As can be seen from the results presented in Figure 10, there is a significant alteration to the magnitude of the simulated U and k profiles at the meteorological mast location for the coarse and medium mesh. However, the effect of further refining the mesh to the fine configuration is only very slight whilst a significant computational expense was incurred as shown in Table 4. As we will only be examining a single direction sector and in order to preserve the academic relevance of the presented analysis, the fine mesh was used for all simulations.

4. Neutral Simulations

The first step in this analysis is to understand the neutral flows before we consider the more complicated events in which stability effects are present. As it was not possible to arrange access to the full set of sonic anemometer data from the Vaudeville site, it was necessary to convert the CFD results for turbulent kinetic energy, k, to Turbulence Intensity, TI, in order to provide a direct comparison to the field dataset. This conversion was achieved by assuming that the flow is fully isotropic and thus:

$$\mathrm{TI}\approx \frac{\sqrt{\frac{2}{3}k}}{\overline{U}}$$

(40)

This calculation was performed for k values at 80 m in the converged CFD simulations in order to provide a comparison with the validation dataset. Values for shear exponent factor, α, were also calculated from the converged CFD simulations between 40 m and 80 m in order to provide a direct comparison with the validation dataset.

4.1. Process

The neutral simulations were configured as described in Section 3.2.

Due to a lack of canopy structural data or a detailed description of the atmospheric boundary layer characteristics, it was necessary to adjust various parameters in the CFD model in order to identify the appropriate settings to simulate the neutral events observed in the validation dataset. Thus, the following variables were adjusted iteratively:

• Reference height, Z_ref
• Reference velocity, U_ref
• Canopy loss coefficient, L_x: Variable h_c
• Canopy loss coefficient, L_x: Constant h_c

When the term ‘Variable h_c’ is used, simulations have been conducted using the canopy height data discussed in Section 2.2. When the term ‘Constant h_c’ is used, simulations have been conducted using a constant canopy height for the forested area. The results of this analysis are presented in the following section.

4.2. Results

4.2.1. Reference Height, Z_ref and Reference Velocity, U_ref

The values set for Z_ref and U_ref are used by WM to calculate the value of $U_{*}$ and also to define the inlet velocity profile. The simulations summarised in

Table 5. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of Z_ref. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
1	40	6.5	0.05	0.09	Variable	0.415	0.142
2	60	6.5	0.05	0.09	Variable	0.415	0.142
3	80	6.5	0.05	0.09	Variable	0.415	0.142
4	100	6.5	0.05	0.09	Variable	0.415	0.142
5	500	6.5	0.05	0.09	Variable	0.415	0.142

and in

Table 6. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of U_ref. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
6	100	5	0.05	0.09	Variable	0.418	0.141
7	100	5.5	0.05	0.09	Variable	0.418	0.141
8	100	6	0.05	0.09	Variable	0.417	0.141
9	100	7	0.05	0.09	Variable	0.417	0.141
10	100	13	0.05	0.09	Variable	0.418	0.142
11	100	20	0.05	0.09	Variable	0.418	0.142

were conducted in order to assess the sensitivity of the model to the prescribed value of Z_ref and U_ref respectively. The default WM value for the canopy loss coefficient, 0.05 m⁻¹, has been used for all simulations. The results of these simulations are also displayed in Figure 11 where they are compared to the validation dataset. The target neutral range is highlighted in green. In all tabular results, the adjusted parameter is highlighted in bold for clarity.

As can be seen from the results presented in Table 5 and Table 6 and in Figure 11, the values of α and TI simulated at the location of the meteorological mast are insensitive to the prescribed value of Z_ref and U_ref. In Figure 11 we see the locus of results in this section indicated as a purple oversized data point, the simulated value of α is in line with the observed value for the neutral events whilst the values of TI are significantly lower. Due to the insensitivity of the model to the prescribed values, it is not possible to correct this discrepancy by adjusting Z_ref or U_ref. This confirms a lack of sensitivity to a change in Reynold’s number when operating at high Reynolds number values in the absence of stability effects or significant separation due to complex terrain downstream of the obstruction.

4.2.2. Canopy Loss Coefficient, L_x: Variable h_c

In these simulations, the sensitivity of the CFD simulation to the prescribed value of the canopy loss coefficient, L_x was assessed using the simulations summarised in

Table 7. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of L_x with a variable canopy height. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
12	100	6.5	0.001	0.09	Variable	0.223	0.095
13	100	6.5	0.01	0.09	Variable	0.373	0.129
14	100	6.5	0.02	0.09	Variable	0.397	0.135
15	100	6.5	0.03	0.09	Variable	0.405	0.138
16	100	6.5	0.04	0.09	Variable	0.411	0.140
17	100	6.5	0.045	0.09	Variable	0.413	0.141
18	100	6.5	0.06	0.09	Variable	0.420	0.144
19	100	6.5	0.07	0.09	Variable	0.423	0.145
20	100	6.5	0.08	0.09	Variable	0.426	0.146
21	100	6.5	0.09	0.09	Variable	0.430	0.148
22	100	6.5	0.5	0.09	Variable	0.484	0.169

. The canopy height was allowed to vary as described by the canopy height data outlined in Section 2.2. The CFD outputs for α and TI at the meteorological mast location are summarised in Figure 12 where they are compared to the validation data.

As can be seen in Figure 12, the CFD simulation is significantly more sensitive to the prescribed value of the canopy loss coefficient. It is possible to bring the simulated value of both α and TI into the desired neutral range by applying a canopy loss coefficient of 0.5 m⁻¹ as used in simulation No. 22.

4.2.3. Canopy Loss Coefficient, L_x: Constant h_c

We now examine sensitivity of the CFD simulations to the prescribed value of L_x when using a constant rather than a variable canopy height. Firstly, the canopy height was set to 11 m which is the average of the canopy height data summarised in Figure 8. The simulations conducted using this height are summarised in

Table 8. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of L_x with a constant canopy height of 11 m. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
23	100	6.5	0.02	0.09	11	0.360	0.130
24	100	6.5	0.03	0.09	11	0.360	0.130
25	100	6.5	0.04	0.09	11	0.363	0.130
26	100	6.5	0.05	0.09	11	0.365	0.131
27	100	6.5	0.06	0.09	11	0.368	0.133
28	100	6.5	0.09	0.09	11	0.374	0.136
29	100	6.5	0.12	0.09	11	0.379	0.138
30	100	6.5	0.15	0.09	11	0.383	0.141
31	100	6.5	0.2	0.09	11	0.389	0.144
32	100	6.5	0.3	0.09	11	0.397	0.148
33	100	6.5	0.4	0.09	11	0.404	0.151
34	100	6.5	0.6	0.09	11	0.412	0.156
35	100	6.5	0.7	0.09	11	0.415	0.158
36	100	6.5	0.8	0.09	11	0.414	0.158

.

The canopy height was then gradually increased to the average value of 30 m stated in [8]. These simulations are summarised in

Table 9. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of L_x with a constant canopy height of 20 m. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
37	100	6.5	0.05	0.09	20	0.458	0.154
38	100	6.5	0.7	0.09	20	0.462	0.176
39	100	6.5	0.9	0.09	20	0.465	0.179

,

Table 10. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of L_x with a constant canopy height of 25 m. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
40	100	6.5	0.05	0.09	25	0.544	0.193
41	100	6.5	0.9	0.09	25	0.570	0.238

and

Table 11. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of L_x with a constant canopy height of 30 m. The adjusted parameter is in italics.

Simulation No.	CFD Settings					CFD Output
	Zref (m)	Uref (m/s)	Lx (m−1)	Cμ	hc (m)	α	TI
42	100	6.5	0.05	0.09	30	0.572	0.174
43	100	6.5	0.7	0.09	30	0.514	0.193
44	100	6.5	0.9	0.09	30	0.515	0.197

. As before, all simulations are compared to the validation dataset in Figure 13.

As can be seen in Figure 13 that the effect of varying the canopy loss coefficient is heavily dependent on the average canopy height used. It is again possible to simulate α and TI values which fall within the desired using certain configurations.

4.3. Discussion

It can be seen in the analysis presented above that the CFD simulation is most sensitive to the prescribed value of the canopy loss coefficient. By tuning this variable it is possible to bring both simulated wind shear and turbulence intensity values in line with values observed during neutral events in the validation dataset. In order to visualise the effect of this tuning on the simulated wind characteristics, profiles are extracted at the meteorological mast location for simulation No. 4 where the default value of L_x is used and simulation No. 38 where the value of L_x has been tuned. In Figure 14, these simulated profiles are presented along with the average profiles of all neutral events in the validation dataset.

As can been seen in Figure 14 that there is little difference between the velocity profile simulated in Nos. 4 and 38. Both simulations show good agreement with the normalised mean velocity profiles in the validation dataset for measurement points above 10 m.

The effect of tuning the prescribed value of the canopy loss coefficient is more clearly evident in the profiles for turbulence intensity where we see that the values simulated in No. 38 are more in line with values in the validation dataset. This is with the exception of measurements at 10 m where the simulated values of turbulence intensity in No. 4 are closer to the mean value observed during the neutral events in the validation dataset. However, values of turbulence intensity simulated in No. 38 fall within the expected range.

Whilst the wind characteristics simulated using the configuration in No. 38 are similar to those observed in the validation dataset, the required value of the canopy loss coefficient is 10 times the default value in WM. Thus, it is prudent to investigate whether the required value has any basis in reality. As mentioned in Section 3, the canopy loss coefficient, L_x, is the product of the canopy drag, C_d, and the Leaf Area Density (LAD), A(z). A value of C_d = 0.15 has been suggested by [25] as being appropriate for a variety of forest canopy types. This would indicate that the average LAD for the Vaudeville forest is approximately 4.6 m⁻¹ if use the value of L_x from No. 38.

In order to set this average LAD value in context, we can examine published values for LAD such as those found in [26]. In this paper, the authors provide a selection of LAD profiles for dense canopies. Whilst peak LAD values of up to 8 m⁻¹ were suggested, values of 0.5–3 m⁻¹ were more common.

Thus, it would appear that an average LAD value of 4.6 m⁻¹ for the Vaudeville forest is high but realistic. Given that we are considering a mixed canopy and that the validation dataset relates to the summer months, this value is plausible.

As shown in [3], the ideal situation when modelling a forest within a CFD domain is to include both realistic canopy height and height dependant LAD data. When such a level of detail is not available, the best option is simply to utilise a constant canopy height and a mean value of LAD. As we were unable to gain access to any level of LAD data for the Vaudeville site, and given the quality of the profiles in Figure 14, the configuration used in simulation No. 38 will be taken as the best option to simulate the neutral events for the Vaudeville site.

5. Stable Simulations

In the previous section, we systematically adjusted the CFD simulation settings in order to model the neutral events observed in the validation dataset for the Vaudeville site. Having simulated the effect of the forest canopy on the wind resource, we now include buoyancy forces in the CFD simulations and attempt to model the stable events.

5.1. Process

The simulations were configured as for simulation No. 38, described in Section 4, with the addition of the physics required to model buoyancy effects as outlined in Section 3 with a domain height of 2 km. The floor temperature was gradually adjusted in order to induce stable stratification of the surface layer. The resulting wind characteristics were then compared to the validation dataset.

5.2. Results

The considered simulations in which stable stratification of the boundary layer was induced are summarised in

Table 12. Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed temperature of the domain floor. The time given is the physical computational time that the simulation required to reach a converged solution. The adjusted parameter is in italics.

Simulation No.	Floor Temperature Difference from Ambient (Kelvin)	CFD Output
		α	TI	Time (min)
48	−0.5	0.594	0.120	800
49	−1	0.626	0.117	828
50	−5	0.720	0.107	1088
51	−10	0.764	0.099	2204
52	−25	0.833	0.092	2434
53	−50	0.864	0.068	2574

. The resulting wind characteristics are compared to the validation dataset in Figure 15 where the target stable range is highlighted in blue. In Table 12, the floor temperature is defined in terms of deviation from the ambient air temperature of 288 K.

5.3. Discussion

As can be seen in Figure 15, decreasing the temperature of the floor in the CFD domain has a profound effect on the wind characteristics in the CFD simulation. The resulting values of α and TI simulated at 80 m are in line with those observed during stable events in the validation data set.

In order to validate the simulated wind profile, values were extracted at the meteorological mast location for simulation No. 51. In Figure 16, these values are compared with the average profile of the stable events in the validation dataset.

As can be seen in Figure 16, the simulated stable wind characteristics at the meteorological mast location are well within the range of values observed during stable events in the validation dataset. However, it is clear that the required temperature differential on the floor surface for the latter simulations, up to 50 K less than the ambient air temperature, is far from what could be reasonably be expected in reality. However, the value of 10 K to achieve the results for No. 51 as presented in Figure 16 is in line with expectations.

6. Unstable Simulations

The next step in this analysis was to attempt to simulate the joint effects of forestry and unstable stratification within the considered domain. The simulations were configured as for simulation No. 38, with the addition of the physics required to model buoyancy effects as outlined in Section 3. The floor temperature was then gradually increased from +0.5 K to +10 K in order to induce unstable stratification of the surface layer. Also, the height of the domain was increased to 3000 m and the inversion height was increased to 1250 m. The obtained results and their general trends were not in line with expectations. The reasons for this shortcoming are unclear, however, the fact the unstable simulations required up 27 times the CPU (Central Processing Unit) time of the neutral equivalent indicate that the turbulence model struggles to capture the joint effects of the forest canopy and non-neutral stability.

As the state of the art develops, the formulation used by WM may be found to be inadequate when considering these complex flow regimes and it may be necessary to use Large Eddy Simulation (LES), Direct Numerical Simulation (DNS) or modifications to the K-Epsilon model such as those proposed in [27,28]. Research is progressing [29] on the use of LES to simulate non-neutral canopy flows. Whilst success has been achieved simulating stable events using these more advanced models, the simulation of unstable events also requires further consideration.

7. Conclusions

It has been shown in this paper that is it possible using a t-RANS CFD model to simulate the joint effects of canopy drag and atmospheric stability when considering stable stratification for a site in North-Eastern France. However, it was not possible to simulate the unstable events in the validation dataset despite modification of boundary layer parameters, further analysis will be required.

The study was limited to a 10-degree direction sector, a 2-month period of reasonably constant canopy density and only considered events where the mean wind speed at 40 m was above 3 m/s. The remaining data displayed a variation in turbulence intensity from 0.01 to 0.45 and from 0.1 to 0.9 for wind shear. By considering neutral stratification and making reasonable assumptions for canopy density and canopy height, a range of 0.1 to 0.25 for turbulence intensity and 0.36 to 0.57 for wind shear could be achieved. This range could be extended to the lower values of turbulence intensity (minimum of 0.068) and higher values of wind shear (maximum of 0.864) by numerically simulating the effects of stable stratification. However, the increased CPU time required for convergence (up to 2574 min for stable simulation compared to c. 480 min for neutral equivalent) and the fact that the unstable simulations were unsuccessful, show that atmospheric stability is the prevailing source of modelling uncertainty for this study.

Whilst the neutral and stable simulations appear to match observations, it would be beneficial to have access to additional validation data in order to assess if the real world flow conditions are being accurately captured numerically at locations other than at the meteorological mast. For example, we have not been able to assess the ability of the numerical simulation to capture the recovery of the flow following the obstruction presented by the forest. This highlights the need for a more comprehensive measurement campaign and detailed characterisation of both the canopy height and LAD variation in order to allow assessment of the quality of numerical simulations when such complex dynamics are present.