Numerical simulation of non-neutral forest canopy flows at a site in North-Eastern France

Abstract. The flow over densely forested terrain under neutral and non-neutral conditions is considered using 10 commercially available Computational Fluid Dynamics software. Results are validated against data from a site in North-Eastern 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 are more difficult to simulate accurately. The sensitivity of the numerical model to parameters such as canopy height, canopy density and the turbulence modelling constant Cμ is also assessed. 15


Introduction
The computational power required to run full Computational Fluid Dynamics (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 [Mortensen & Jørgensen (2013)].Following this trend, research activities have increased into the flow dynamics generated by non-trivial terrain and atmospheric features in order to fully 20 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 Brower et al. (2014) that forestry increases modelling uncertainty in terms of root mean square error by a factor of 4-5 when modelling the flow between mast pairs using a variety of 25 industry standard modelling software packages.In Desmond & Watson (2014) 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 [Brunet & Irvine (2000), Morse et al. (2002) and Chaudhari et al. 30 (2016)].
In Desmond et al. (2017) 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 Desmond et al. (2017) 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 35 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 North-Eastern 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, Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.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 package which is a front end for the ANSYS CFX flow solver.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.5 In Section 2 of this paper the validation data are presented and in Section 3 the CFD configurations are discussed.
Sections 4, 5 and 6 present the validation results and discussion for the neutral, stable and unstable intercomparisons respectively.Section 7 presents overall conclusions.

Validation data
This study uses data from a mast near Vaudeville which is located adjacent to a wind farm in North-Eastern France 10 (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.Data were provided between 01/01/10 -31/12/11 as 10 min averages from a series of sonic anemometers (Metek 5 USA-1), temperature sensors and wind vanes on a 100 m mast as summarised in Table 1.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.In order to isolate non-neutral data with which to validate the CFD simulations, the steps described in Section 2.1 were taken.

Wind speed data
The 250 -260˚ direction sector was examined in order to limit variations in the proximity of the mast to the forest edge.This range can be seen in Figure 3. 5 The effect of the forest canopy on the wind resource will vary seasonally and annually as the trees grow and 10 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 Desmond and Watson (2014) 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 15 classes.A summary of the available data is given in Figure 4.The next step was to apply the methodology outlined in Desmond and Watson (2014) to identify non-neutral 10 events.Thus turbulence intensity (TI) at 80 m and wind shear between 40 m and 80 m were calculated using Eq. 1 and Eq. 2 respectively.

𝑇𝐼 = 𝜎 𝑢 𝑈 ̅
(1) As can be seen in Figure 5, the values of the observed wind shear and Turbulence Intensity become less sensitive 20 to solar irradiance levels at higher wind speeds.Following Desmond and Watson (2014) 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 5 selected data set in order to identify stable, neutral and unstable events as shown in Figure 6.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.

Canopy height data
Data were provided for a 5 km radius around the Vaudeville meteorological mast by Intermap Technologies Ltd.
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 10 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 [Texier et al. (2010)].Information on the canopy measurement techniques can be found in Andersen et al. (2008).A comparison between a sub set of the forested area determined using the Intermap data and an aerial photograph of the area is provided in Figure 8.The distribution of canopy heights in the examined region is shown in Figure 9. 5 The mean canopy height in Figure 9 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 10 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 Desmond et al. (2014).

CFD configuration 15
As stated previously, all CFD simulations in this paper were configured using the WM front-end to the CFX software.All simulations were run using the RANS configuration of the CFD software with the Shear Stress Transport (SST) turbulence closure (Menter, 1994).Specific details of the configuration used are given below.

Domain description
Circular computational domains are generated in WM where the outer edges are divided into 24 surfaces, as shown in Figure 10, which allows various wind directions to be considered using a single domain configuration.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 5 direction sector investigated, shown in Figure 3.
Topographical details from the 90 m resolution Shuttle Radar Topography Mission (SRTM) [Werner (2001)], 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 10 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 zo= 0.04 m which is what would be expected for a site containing low grass [Burton et al. (2001)].Values of 0.1 m and 0.001 m were also trialled, however the impact on results was negligible.This is due to the fact that much of the fetch along the 255˚ direction 15 is occupied by forestry and thus the surface roughness itself will have a reduced role in dictating the wind characteristics.
Following the analysis in Desmond et al., (2014), the SST turbulence closure model was used for all simulations (Menter, 1994).
The height and extent of the Vaudeville forest was described by a set of x, y, z coordinates derived from the 20 Intermap data described in Section 2. 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 sources and sinks are applied to all control volumes which are identified to be located below the top of the forest canopy.
In the momentum equations, the drag term is: 25 Where, Fi Drag force per unit volume in the i-direction (kg/m 2 .s 2 )

A(z)
Leaf area density at height z (m -1 )

|𝑈|
The modulus of the windspeed (m/s)

𝑈 𝑖
The wind speed in the i-direction (m/s)

𝐶 𝑑
The canopy drag coefficient (dimensionless) 30 The effect of the forestry drag on turbulence quantities has been modelled and discussed by various authors [e.g.Svensson & Haggkvist (1990) , Lopes da Costa (2007), Katul et al. (2003) , Sogachev (2009) For the turbulence kinetic energy equation Where, β p and β d are constants, the values of which are given in Table 3.
For the dissipation rate equation 5 Where, C ε4 and C ε5 are constants, the values of which are also given in Table 3.
When using the SST turbulence model, source/sink terms are required for the equation for the turbulence eddy frequency instead of .The required source terms are as follows: For the turbulence frequency equation This source term for the  equation is derived from the relationship 10 which itself results from the transformation of the  equation into the equation for , via the identity  =   .
A discussion on the formulation of these equations for a  −  model can be found in Lopes da Costa (2007) and Sogachev (2009).The appropriate value for the modelling constants in the above equations has been an area of some research.For this research, the values as recommended by Lopes da Costa (2007), are used and these are summarised in Table 3. 15

Constant Value
β p 0.17 The porosity of the canopy was defined by a loss coefficient, Lx, which is the product of the canopy drag, Cd, 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 20 WM implementation of the drag force and turbulence source are based on a definition of a drag force including a 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.5.

Neutral 5
The inlet wind profile is defined by prescribing a reference mean horizontal wind speed, Uref, and the height above ground level at which it occurs zref 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 ( 8. The mean horizontal wind speed profile is then calculated using the standard form of the log law along with profiles for Turbulent Kinetic Energy k(z) and Turbulent Eddy Dissipation rate ε(z) using the standard form of 10 the Richards and Hoxey equations [Richards & Hoxey (1993)].These profiles are calculated using a default C μ value of 0.09.
The calculated U(z), k(z) and ε(z) profiles are then applied as boundary conditions to the 12 outer surfaces which constitute the inflow region at the outer edge of the domain.The other 12 outer surfaces constituting the outflow along with the top of the domain are set up as entrainment openings.For these boundaries, a 0 Pa relative pressure 15 is defined rather than a wind speed.
Finally, the floor of the domain is set as a no-slip wall with roughness equal to the equivalent sand grain roughness,   , which is derived from z 0 using Eq. 8. [McCormick et al. (2012)]:

Stable
When the effects of stability are included in a CFD simulation configured using WM, a stable capping layer, 20 which begins at the top of the boundary layer, is introduced.This condition is representative of the situation in the real atmosphere and has been found to be important for the convergence and relevance of simulations in which the effects of stability are included [Montavon (1998)].Above the boundary layer, the potential temperature gradient is set to 3.3 x 10 -3 K/m by default as per the standard US atmosphere [U. S. Standard Atmosphere (1962)].
The height of this inversion, which is equal to the boundary layer height h, is calculated as a function of the 25 Coriolis parameter, inversion gradient, U * , and the Obukhov length [Garratt (1992)].
The value of U * is estimated using a look up table [Garratt (1990)] which provides reference values for a given Surface Rossby number [Stull (1988)], which itself is a function of the free-stream wind speed, the local Coriolis force, f, and the surface roughness as shown in (10: Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.

𝑅𝑜𝑠𝑠𝑏𝑦 =
×  0 (10) Where,   is the free steam of geostrophic wind speed.The value of the Obukhov Length is set to a default value of 1000 m in WM.Once values for U * and h have been calculated, it is possible to describe profiles for the mean horizontal and vertical wind speeds, the Turbulent Kinetic Energy and Turbulent Eddy Dissipation rate using the stability dependant form of the relevant equations as defined in Zilitinkevich et al. (1998).
Buoyancy forces are also introduced into the domain in order to allow stability effects to develop within the 5 domain.These forces are related to a reference temperature of 288 K and a reference location which coincides with the start of the capping stable layer.
Once again, a pressure condition rather than a velocity is set at the boundaries which the software defines as openings.The pressure profile at the outlet is set up to represent an atmosphere in hydrostatic equilibrium with a temperature profile as prescribed at the inlet.10 Once these calculations have been performed, the appropriate velocity, temperature, k and ε values are set for each boundary and are also used as the initial conditions for the entire domain.In addition, it is possible to set a temperature condition on the floor in order to introduce stability effects in the surface layer.
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 15 using: The initial time step is set to 10 seconds and increases to a maximum value of 30 seconds depending on how quickly the simulation converges.

Unstable
The configuration for modelling unstable events is as summarised in Section 3.3 with the exception that the 20 temperature set on the floor surface is greater rather the ambient air temperature in order to simulate the effect of solar heating of the earth's surface.
In theory, this configuration and the adjustment of the surface boundary temperature should allow the full range of atmospheric stabilities within the surface layer to be investigated.

Mesh sensitivity 25
A mesh sensitivity study was conducted using a neutral configuration.A constant canopy loss coefficient of Lx = 0.05 m -1 was used for all simulations along with Uref = 6.5 m/s at zref = 40 m.The circular domain generated by WM is divided into nine zones for the purposes of meshing as shown in Figure 10.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 10 (a), the critical dimensions which define the mesh are shown.For all simulations the following values 10 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.

15
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.All simulations were conducted on a High Performance Computing (HPC) cluster which consists of 161 nodes, each having two six-20 core Intel Westmere Xeon X5650 CPUs and 24GB 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.As can be seen from the results presented in Figure 11, there is a significant alteration to the magnitude of the simulated U and k profiles at the mast location for the coarse and medium mesh.However, the effect of further 10 refining the mesh to the fine configuration is 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.

Neutral simulations 15
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: 20 This calculation was performed for simulated k values at 80m and values for α were calculated between 40 and 80 m in order to provide a direct comparison with the validation dataset.

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 5 iteratively: When the term 'Variable hc' is used, simulations have been conducted using the canopy height data discussed in Section 2.2.When the term 'Constant hc' is used, simulations have been conducted using a constant canopy height for the forested area identified in Figure 8. 15 The results of this analysis are presented in the following section.

Reference height, zref and Reference velocity, Uref
The values set for zref and Uref are used by WM to calculate the value of U * and also to define the inlet velocity profile.The simulations summarised in Table 5 and in Table 6 were conducted in order to assess the sensitivity 20 of the model to the prescribed value of zref and Uref respectively.The default WM value for the canopy loss coefficient, 0.05 m -1 , has been used for all simulations.
The results of these simulations are also displayed in Figure 12 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, Table 6 and in Figure 12 that the values of α and TI simulated 5 at the location of the mast are insensitive to the prescribed value of zref and Uref.In Figure 12 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 zref or Uref.This confirms a lack of sensitivity to a change in Reynold's number when operating at high Reynolds number values in the absence of 10 stability effects or significant separation due to complex terrain downstream of the obstruction.

Canopy loss coefficient, Lx : Variable hc
In these simulations, the sensitivity of the CFD simulation to the prescribed value of the canopy loss coefficient, Lx was assessed using the simulations summarised in Table 7.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 mast location are 15 summarised in Figure 13 where they are compared to the validation data.7.
As can be seen in Figure 13, 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 -1 as used in simulation No. 22.

Canopy loss coefficient, Lx : Constant hc 10
We now examine sensitivity of the CFD simulations to the prescribed value of Lx 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 9.The simulations conducted using this height are summarised in Table 8.The canopy height was then gradually increased to the average value of 30 m stated in Texier et al. (2010).These simulations are summarised in Table 9 to Table 11.As before, all simulations are compared to the validation dataset in Figure 14.   8 to Table 11.
As can be seen in Figure 14 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.

Cμ 10
It can be seen from the preceding simulations that the simulated TI tends frequently to be lower than that observed.
In an attempt to increase levels to the desired range, simulations were conducted in which the value of Cμ was reduced, an approach which has yielded satisfactory results in Gorle et al. (2009) and Richards and Hoxey (1993).
The considered simulations are summarised in Table 12 and are compared with the validation data in Figure 15   As can be seen in Figure 15 it is possible to increase the simulated values of Turbulence Intensity by adjusting the value of Cμ.However, the required value, 0.0005, is well below the range considered in Gorle et al. (2009) and 5 Richards and Hoxey (1993).

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 10 the effect of this tuning on the simulated wind characteristics, profiles are extracted at the mast location for simulation No. 4 where the default value of Lx is used and simulation No. 38 where the value of Lx has been tuned.
In Figure 16, these simulated profiles are presented along with the average profiles of all neutral events in the validation dataset.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, Lx, is the product of the canopy drag, Cd, and the Leaf 15 Area Density (LAD), A(z).A value of Cd = 0.15 has been suggested by Amiro (1990) as being appropriate for a variety of forest canopy types.This would indicate that the average LAD for the Vaudeville forest is approximately 3 m -1 .
In order to set this average LAD value in context, we can examine published values for LAD such as those found in Banerjee et al. (2013).In this paper, the authors provide a selection of LAD profiles for dense canopies.Whilst 20 peak LAD values of up to 8 m -1 were suggested, values of 0.5 -3 m -1 were more common.Thus, it would appear that an average LAD value of 3 m -1 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 in line with expectations.As shown in Desmond et al. (2014), 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 16, 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.10

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

Results
The considered simulations in which stable stratification of the boundary layer was induced are summarised in Table 13.The resulting wind characteristics are compared to the validation dataset in Figure 17 where the target stable range is highlighted in blue.In Table 13, the floor temperature is defined in terms of deviation from the ambient air temperature of 288 K. 20

Simulation
No.

Discussion
As can be seen in Figure 17 that decreasing the temperature of the floor in the CFD domain has a profound effect 5 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 mast location for simulation No. 51.
In Figure 18, these values are compared with the average profile of the stable events in the validation dataset.As can be seen in Figure 18, the simulated stable wind characteristics at the mast location are well within the range of values observed during stable events in the validation dataset.However, it is clear that the required

Unstable simulations
The next step in this analysis is to attempt to simulate the joint effects of forestry and unstable stratification within the considered domain.5

Process
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 in order to induce unstable stratification of the surface layer.The resulting wind characteristics were then compared to the validation dataset.10

Results
For the unstable simulations, the CFD domain was initially sized as described in Section 3 with an overall domain height of 2000 m.The height of the stable capping inversion was set to the height of the boundary layer h = 350 m.Initial simulations were conducted using this standard WM configuration and are summarised in Table 14.The 15 corresponding simulated values of α and TI at the mast location are compared with the validation data in Figure 19 where the target range is highlighted in red.

Simulation
No.  As can be seen in Figure 19, it was not possible to simulate the effect of unstable stratification for the Vaudeville site using the standard WindModeller configuration.Whilst the results of the simulations presented in Section 4 5 and Section 5 showed very defined trends, which were in line with expectations, the unstable simulations results

Domain height (m)
shown in Figure 19 are highly erratic.These inconsistent results may be due to the fact that the free stream stable temperature gradient, imposed above 350 m within the domain, is preventing unstable flow characteristics from emerging in the surface layer below.In order to remove this possible conflict, the height of the capping layer was moved to a height of 1250 m.The height of the simulation domain was also increased to 3000 m in order to allow 10 a sufficiently deep stable capping layer to exist in the new configuration.
Firstly, the height of the domain was increased to 3000 m and simulations were conducted using an identical configuration to simulations No. 54 -57 above.This first step was performed in order to ascertain to what degree the results were modified by the increased domain height.The simulations conducted using this configuration are summarised in Table 15 and the resulting α and TI values are compared to the validation data in Figure 20.15

Simulation
No.   As can be seen in Figure 20, increasing the domain height to 3000 m has a moderate effect on the simulated values of α and TI at the mast location.Simulations were now conducted in which the height of the capping inversion 5 was moved to 1250 m within the 3000 m simulation domain.Simulations conducted using this configuration are summarised in Table 16 and the results are summarised in Figure 21.

Discussion
As can be seen in Figure 19, Figure 20 and Figure 21 it was not possible to simulate the effect of unstable 5 stratification for the Vaudeville site using the WM configuration despite reasonable modifications.
As previously outlined, when the magnitude of the unstable stratification is increased we would expect a reduction in the value of α along with increased levels of TI.This trend is only observed when the floor temperature is increased by 0.5 to 1 K as shown in simulation No. 54 -55 and No. 58 -59.For higher levels of temperature differential between the floor and the ambient air, simulated values of α and TI are inconsistent with these 10 expectations.
It was not possible to bring the unstable simulation results in line with the validation data by adjusting the height of the stable capping inversion with the results presented in Figure 21 being broadly similar to those in Figure 20.

7 Conclusions
It has been shown in this paper that is it possible using a 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.It may be possible to achieve the desired results by adjusting the boundary conditions 20 to better represent the effects of unstable stratification and/or tuning the modelling constants.However, due to the fact that validation data is limited to a single measurement location, it will not be possible to fully appreciate the ramifications of such alterations on the overall quality of the simulation.
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 25 numerically at locations other than at the 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.
Whilst considering non-neutral events in micro-meteorological studies will allow a greater appreciation of the prevailing wind characteristics at a given site, it is clear that this knowledge will come at a significant 5 computational expense.The presented stable and unstable simulations in this paper required approximately 27 times the computational time in order to achieve a converged solution compared to analogous purely neutral simulaitons.Whilst this is a considerable expense, it may be unavoidable if we are to reduce the uncertainty generated by non-neutral events in forested terrain.

Figure 1 .
Figure 1.Location of the Vaudeville met mast is indicated by the red marker.[Picture credit: www.maps.google.com]

Figure 2 .
Figure 2. IGN map of the Vaudeville region.The meteorological mast location (46° 26' 58''N , 05° 35' 02''E) is marked with by the red X circumscribed by a red circle.Turbine locations are indicated by a red inverted Y. [Texier et al. (2010)].

Figure 3 .
Figure 3. Aerial photograph of the Vaudeville site showing the 250-260˚ direction sector.Distances to the forest edge are indicated by the red arrows.[Picture credit: www.maps.google.com]

Figure 4 .
Figure 4.A summary of the available data showing the number of observations and the maximum recorded irradiance level for each month.Month 1 relates to January 2010.The yellow shading identifies the months selected for analysis.5

Figure 5 .
Figure 5. Observed wind shear and Turbulence Intensity at the Vaudeville site for the 250-260˚ direction sectors for July and August 2010.The red lines indicated the applied neutral threshold values.Turbulence intensity values are calculated at 80 m.

Figure 6 .
Figure 6.Neutral thresholds are applied to the selected data.Points in the sector with the green background are considered to be neutral, blue are stable and red unstable.Profiles for the oversized data points in each of these sectors are given in Figure 7.Only events with wind speeds of > 3m/s at 40 m are displayed in this figure.

5Figure 7 .
Figure 7. Sample profiles for the oversized data points in Figure 6.

15 WindFigure 8 .
Figure 8.The red dot indicates the position of the meteorological mast in each image and dashed line indicates the centre of the 250 -260˚ direction sector.The scales used in these images are similar.Figure (a) shows only trees of heights between 2.5 -5 m.This setting identifies the perimeter of the forestry.

Figure 9 .
Figure 9.The height distribution of trees for the region outlined in Figure 8.
], often in the context Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License. of  −  models.For such models, sources and sinks are added to the turbulence kinetic energy  and turbulence dissipation  equations, to model the added turbulence and redistributed lengthscales as follows: Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.factor ½. Lopes da Costa (2007) omits the factor ½ in his definition.As a consequence, the loss coefficient in Lopes da Costa (2007) is to be interpreted as half the loss coefficient in WindModeller.

Figure 10 .
Figure 10.(a) Mesh zones created by WM.The red dot in (b) indicates the mast location.The same domain is detailed in (a) and (b).

Figure 11 .
Figure 11.Results of the mesh sensitivity study.

Figure 12 .
Figure 12.The locus of the results of Simulations 1-11 are represented by the purple oversized data point.

Figure 15 .
Figure 15.The results of Simulations 45-47 are represented by the purple oversized data points.The reference numbers shown correspond to the simulation numbers given in Table12.

Figure 16 .
Figure 16.Graphs showing the simulated normalised velocity and Turbulence Intensity profiles at the mast location for simulations No. 4 & No. 38.The field data points represent the average value at that height for all neutral events whilst the horizontal bars indicate the range of recorded values at each height in terms of 2 x Standard Deviation.As can been seen in Figure 16 that there is little difference between the velocity profile simulated in No. 4 and 38.Both simulations show good agreement with the normalised mean velocity profiles in the validation dataset for 5 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 10 of Turbulence Intensity simulated in No. 38 fall within the expected range.
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 time that the simulation required to reach a converged solution.Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 17 .
Figure 17.The results of Simulations 48-53 are represented by the blue oversized data points.The reference numbers shown correspond to the simulation numbers given in Table13.

Figure 18 .
Figure 18.Graphs showing the simulated normalised velocity and Turbulence Intensity profiles at the mast location 10 Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.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.

Figure 19 .
Figure 19.The results of Simulations 54-57 are represented by the red oversized data points.The reference numbers shown correspond to the simulation numbers given in Table14.
Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed temperature of the domain floor with a domain height of 3000m.∆ Temp is the domain floor temperature minus the ambient air temperature.Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 20 .
Figure 20.The results of Simulations 58-61 are represented by the red oversized data points.The reference numbers shown correspond to the simulation numbers given in Table15.
Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed temperature of the domain floor with a domain height of 3000m and altered capping inversion height.∆ Temp is the domain floor temperature minus the ambient air temperature.10 Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 21 .
Figure 21.The results of Simulations 62-68 are represented by the red oversized data points.The reference numbers shown correspond to the simulation numbers given in Table16.

Table 5 . Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of zref.
25 Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.

Table 7 . Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of Lx with a variable canopy height.
Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.

Table 12 . Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed value of Cμ. 15
.

Table 14 . Summary of simulations run to investigate the sensitivity of the CFD model to the prescribed temperature of the domain floor with a domain height of 2000m. 20
Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2017-34Manuscript under review for journal Wind Energ.Sci. Discussion started: 22 August 2017 c Author(s) 2017.CC BY 4.0 License.