Modelling of wind turbine wakes over forests along the diurnal cycle

This work presents a methodology for the Large-Eddy Simulation (LES) of the continuous transition of atmospheric stability over forests along the diurnal cycle and its effect on the turbulence characteristics of wind turbine wakes. The forest is modelled as a porous surface where temperature changes, transferred to the air via sensible transport, are caused by the variation of net radiation and in proportion to the tree height and leaf density. The flow is driven by a pressure gradient including Coriolis forcing to allow for the development of nocturnal inertial oscillations. An actuator disk is employed to model the wake of a wind turbine located in Ryningsnäs, Sweden, for which metmast measurements are available to carry out a comparison. Results show a good prediction of the inflow and wake characteristics during daytime whereas turbulence fluctuations seem to be overestimated during night periods, attributed to a combination of an excess in geostrophic velocity and coarse mesh resolution. Observations of velocity, heat flux, potential temperature, velocity spectra and other higher order statistics are used to characterize the diurnal variations both in the inflow and across the wake. The results show that the model is capable of representing the turbulence flow dynamics during the diurnal stability transition, hence laying the ground to future studies to assess the performance of wind parks over forested areas.


Introduction
The character of wind turbine wakes is crucial to the power output of wind farms.Wakes from upstream turbines reduce the wind speed and increase the incoming turbulence for those further downstream, which reduces the power output and increases the mechanical loads over the rotors.The recovery of wind speed within the wakes is largely controlled by the turbulence transport of momentum, which in turn is governed by background turbulence and the turbulence generated by the wake itself.An essential fact is that the background turbulence is strongly influenced by the changes in stratification in the Atmospheric Boundary Layer (ABL) along the day.Compared to neutral stratification, an increase in turbulence occurs during the day due to the ascending heat flux from the ground, enhancing mixing and thus reducing the wind shear.During nights, the negative buoyancy suppresses turbulence resulting in a much increased shear relative to neutral or unstable conditions.Broadly, periods of unstable and stable stratification induce a faster wake recovery during the day and longer, less turbulent wakes during nights, respectively [1,2].These variations lead to fluctuations in power production and rotor loads which are needed to be estimated to ensure the optimal performance of the wind park.
Forested regions constitute an area of increasing interest for the development of wind projects due to increasing hub heights, less problems with social opposition in forested areas and synergy effects with forest management, such as shared costs for access roads.The increase in shear and mechanical turbulence modify the wake characteristics compared to those over sites with low or absent vegetation.Several studies point at the distinct effect that a fetch of forested terrain imprints to the wind at a given location, making possible to relate the footprint of the forest distribution and terrain variation at different heights [3], as well as their effect on the wake characteristics [4,5] and rotor loads [6].Recent investigations have been dedicated to combine the representation of forest with the numerical simulation of stratified flows to observe the integrated effects on the development of wakes.In many cases, these works focus on the representation of discrete, stationary stability regimes by employing a fixed heat-flux to compare the stable and/or convective wind conditions above uniform forests [7,8], finding a comparable effect of the stability on the wind flow over a background of increased canopy-produced wind shear.However, by only considering stationary periods of stratification, some features of the ABL highly relevant to the assessment of wind resources are missed, such as nocturnal low level jets and turbulence in the residual layer.Thus, to better estimate the impact of atmospheric stratification, it is needed that numerical models incorporate a realistic representation of the transient nature of the stability changes.
Over flat terrain, the impact of an ideal (symmetrical) variation of heat-flux over a single wake has been studied by [9], with [10] adding obstacles over the ground to represent an urban landscape whereas [11] investigated the diurnal characteristics of wakes within a wind park.Based on methods from [7], a modelling approach for the simulation of diurnal cycles over forests has been presented by [12] for a 1D column model using RANS, laying the ground for the development of numerical approaches to represent the diurnal changes of the ABL over forests.The impact of diurnally varying heat flux on wind turbine wakes is particularly interesting over forests, since compared to other surfaces there is both more surface area and increased flow around the surface elements which leads to very rapid transitions in stratification due to the effective heat transfer.
This work presents a methodology to simulate the transient characteristics of ABL turbulence over forests along the diurnal cycle and their impact over the wind turbine wake.By making use of simple assumptions, it presents a method to drive the diurnal changes in the thermal stratification by means of the heating and cooling of the forest to represent the ensuing changes in the wake turbulence field as recorded in a measurement campaign in Ryningsnäs, Sweden.The methodology employs various techniques found in the literature but presents for the first time -to the knowledge of these authors-an integrated numerical approach for microscale simulations that permits to represent the turbulence dynamics of wake turbulence resulting from the combined backdrop of forests and the transient diurnal stratification.

Measurements
The comparison with the numerical simulation is based on data from the Ryningsnäs site in Southern Sweden (57 • 16 ′ 34.40" N, 15 • 59 ′ 11.69" S).The site is characterized by a heterogeneous forest cover predominately consisting of pine and includes a 140 m high meteorological tower and two wind turbines, one 203 m south of the measurement tower and one 234 m north east of the tower.The southern turbine has a hub height of 80 m and the north eastern one a hub height of 100 m, both with a rotor diameter of 90 m.Only the south turbine is considered in the LES so only measurements from that wake are employed.The height difference of the ground between the southern turbine location and the tower location is 2 m.
The comparison is based on 30 minute block averages which have been subsequently grouped according to the time of the day.The measurements consist of 6 levels of sonic anemometers (METEK USA-1) at the heights 40, 59, 80, 98, 120 and 137.7 m above ground level.The post processing of the 20 Hz sonic raw data is the same as in [13].In order to map the complete shape of the wake, different inflow angles were used from the sonic at 80 m height so data was conditionally sampled based on the direction between the mean flow and the upstream turbine following [4].In other words, it was approximated that the mean wind was accurately measured at the met tower and that the mean wind was directly perpendicular to the turbine rotor plane.By systematically varying the selection for the incoming direction, the whole wake character was mapped.This procedure results in a grid equivalent to an arched cross-section parallel, but situated downstream, to the rotor plane, see Fig. 1-(b).It is important to notice that unlike with the LES, each point in this grid is made up from the average of several, temporally separated and independent 30 minute periods.
In order to assess differences and similarities between the modelled and measured turbulence, spectral densities were computed.This was done in the same way for both simulated time series and measurements.First, the data was split into 30 minute segments, after that a Fast Fourier Transform was computed and spectral densities were constructed by multiplication of complex conjugates of the velocity components.To reduce scatter, the individual spectra were first averaged in logarithmically spaced frequency bins.To facilitate comparison, all individual spectra were normalized by the integral of the concurrent power spectrum of the streamwise velocity component, while the frequency was normalized with the height and the individual 30 min wind speed.Finally all 30 minute periods satisfying the selection criteria were averaged.
To filter out data with similar diurnal cycles in the radiative forcing, measurements of the net radiation balance were used.The data was provided by a net radiation sensor at 40 m above ground level (Kipp and Zonen, NR-Lite).In order to determine corresponding radiative forcing for the LES, the diurnal cycle of the net radiation was determined by averaging over all the occasions that satisfied the quality assessment of [13] and and had a wind speed larger than 4 m/s at hub height.In order to avoid too large deviations in the diurnal cycles only data within one month before and one month after the equinoxes were used.The averaged diurnal cycle of net radiation was then parameterized by an exponential function.Both the parameterized curve and measurements can be seen in Fig. 1-(c).To contrast the character of the flow between day and night, the periods of 21:00 to 03:00 to represent nocturnal conditions and 10:00 to 16:00 to represent daytime conditions were chosen.

Flow model 3.1. Model description
The flow model is based on Large-Eddy Simulations (LES) and it is implemented on the platform OpenFOAM 3.0.1 [14] (the OpenFOAM Foundation Ltd, London, UK).This consists of an incompressible solver that makes use of the Boussinesq approximation to simulate the stratification of the ABL flow.The principal approach to represent the effect that the various components of the simulation have on the flow is by introducing source terms in the momentum and temperature equations [15].Firstly, the forest is conceived as a porous surface extracting momentum from the flow, represented as where a corresponds to the frontal leaf area density (considered here equal to Plant Area Density, PAD), C D = 0.2 is the forest drag coefficient.Make note that the LES filtered/resolved quantities are represented here by ' • ', temporal averages by ' • ', spatial means by '⟨ • ⟩' and magnitudes by ' | • | '.The heating and cooling of the forest is simulated with a source term in the temperature equation [7]: where q r is the heat flux by emitted or absorbed by the forest, calculated as a function of the net radiative flux at the top of the canopy Q h and its downwards penetration within the forest.This, in turn, is a function of the extinction coefficient of light taken as γ = 0.6 and the downward cumulative leaf area index A c , obtained from the downwards vertical integration of the forest density: where h f is the forest height.The Sub-Grid Scale (SGS) model of Deardorff [16] as revised by Moeng [17] is used as it presents a practical way to relate the calculation of the turbulent viscosity ν t = 0.1 √ k SGS l to the stratification regime accounted for by the lengthscale: where g is the gravitational constant, k SGS is the subgrid turbulence kinetic energy and θ, θ 0 are the potential and reference temperatures, respectively.This choice makes l equal to ∆ = (∆ x ∆ y ∆ z ) 1/3 (coinciding with the implicit LES filter) except for only the most stable stratification.This can be problematic when grid resolutions are relatively fine (in the order of a few metres) where the value of l can remain equal to ∆, as discussed by [18].However, the presence of the forest enhances the creation of turbulence even during nights, maintaining a larger size of eddies compared to a terrain without vegetation.Therefore, this problem is not expected in the present setup.The subgrid dissipation ε SGS -also dependent on the stratification-is supplemented by another term that accounts for the increased dissipation due to the forest [19] ε SGS,f : where the value between parenthesis of eq. ( 5) is made equal to 3.9 only at the first grid above the wall.Subgrid turbulence production, on the other hand, is also buoyancy dependent next to the shear generation.The revised Deardorff model of Moeng [17] has been previously implemented in OpenFOAM [20] within the SOWFA package [1].The corresponding libraries have been ported into the current implementation and modified to account for the forest modelling, as presented above.The same applies to the wall model based on the Schumann surface stress concept [21], used in combination with the forest drag model.The calculation of the heat-flux at the wall is also based on a library from SOWFA.

Rotor modelling
The rotor is modelled as Actuator Disk (AD), with a momentum source acting solely on the axial direction.The simplest conception is used here, where the thrust is uniformly distributed over the rotor area A AD according to where U 0,ref corresponds to the undisturbed (in the absence of the rotor) reference velocity at the location of the rotor hub.This velocity is highly turbulent due to the forest, besides varying during the diurnal cycle.An estimation is given by the expression following the work of [22], where ⟨u s ⟩ plane corresponds to the horizontal velocity magnitude averaged over the rotor plane.Since the previous equation is only applicable while C T < 1, the expression is complemented by (Jens N. Sørensen, personal communication), which is used for C T > 0.9.Eq. 9 is derived assuming that C T follows a 2nd-order polynomial in terms of the axial induction .The wind turbine located in Ryningnäs during the measurement campaign corresponds to a 2.5 MW Nordex N90.Information of C T vs incoming velocity is, at the time of writing, available online.Since this velocity corresponds in turn to U 0,ref , this and C T are calculated based on an iterative process of eq. ( 8) or (9).The convergence criterion is set as 001.The cut-in speed is set at 3.25 m/s so for wind speeds below this value, C T is set to 0.15.Hub and tip corrections are also applied using the formulations from [22].The tip correction applied (Prandtl) requires information of the Tip-Speed Ratio (TSR) which could not be found, so its value was estimated as 10.To avoid spatial oscillations in the vicinity of the AD that appear in collocated grids due to the pressure-velocity decoupling, the force in eq. ( 7) is distributed in the axisymmetric direction by taking its convolution with a Gaussian distribution within ±3σ as in [23], with an standard deviation of σ = 2∆ n, where ∆ n is the cell size projected in the axisymmetric direction of the rotor.A simple steering mechanism is devised to maintain the rotor plane perpendicular to the incoming flow.This is based on the determination of the vector normal to the AD plane nAD , computed by means of sampling the horizontal velocity vector u s i = (u x , u y ) at a fixed location p, located at 3D upstream from the AD centre and at hub height (see Fig. 1-(b)).Instead of allowing for an instantaneous alignment of the rotor yaw, a more realistic representation is obtained by employing a moving average in the calculation of u s i , with a window of 10 s, which is then used to calculate nAD .The thrust calculated from eq. ( 7) is applied by introducing a corresponding source term in the LES momentum equation, at the cells centres within the rotor perimeter, extracting momentum only along the direction nAD .Next to the tip correction, no other treatment is applied to the edges of the AD volume.

Numerical setup
LES simulations are carried in box with dimensions L x ×L y ×L z = 7 km × 23.5 km× ∼ 3.88 km.
In the horizontal plane, the domain consists of 3 areas with different cell resolution, as seen in Fig. 1-(a).In the innermost region L ′ x × L ′ y = 1 km × 1.5 km with ∆ x = ∆ y = 5 m.The outer edges (outwards from the dashed lines in Fig. 1-(a) to the domain perimeter) have a width of 4 km in the x−direction and 8 km in the y−direction with ∆ x = ∆ y = 50 m.In between these regions, the cells are uniformly stretched in along x and y.In the vertical direction, the mesh is composed by two sections, with the bottom part defined by a very small stretching of ∼ 1.0054 in order to maintain the cell height approximately uniform within the rotor area.This lower region extends from the flat ground set at the height of the rotor location, i.e. z asl = 122.67m up to z asl = 505.49m with the first cell height ∆ z 1 = 3.93 m.Above this region, cells stretch vertically at uniform ratio of 1.1.This yields ∆ z N /∆ z 1 = 82.41 for the last to first cell size ratio and N z = 115 cells in the vertical direction.The total number of cells in the mesh is N ≃ 40.68 × 10 6 .
The AD is located at 10.5 km from the inlet, at the origin of coordinates in Fig. 1-(a), at 500 m from the start of the fine resolution region and in the middle of the spanwise dimension.The AD centre is at z agl = 80 m.The domain is aligned in the direction of the incoming wind at φ = 181 deg as shown in Fig. 1-(b), producing a difference of only 1 deg with respect to the y−axis.Due to the small stretching in the vertical direction, the cells are maintained A forest of uniform leaf density and constant height is represented in the LES using a PAD field of 0.15 1/m from the bottom surface up to a height of h f = 20 m.These values correspond to approximately the average quantities within the incoming wind direction.At the edges of the domain, PAD is increased to 2.0 1/m within a span of 1 km on the inlet and outlet edges and 0.5 km at each side.This large density value has been observed to aid to erase the forest footprint and streaks in the recirculating flow when periodic conditions are employed [3,23] -potential issues with an internal boundary layer are avoided due to the large separation between this edge and the rotor and wake regions.The development of streaks is also hindered by the long streamwise extension of the domain.The bottom surface is set as a no-slip wall with a uniform roughness of z 0 = 0.1 whereas periodic boundary conditions are set at the sides.
The flow is driven by a constant pressure gradient, calculated from a given geostrophic velocity vector u g .The pressure gradient includes a Coriolis forcing calculated for the latitude of 57 deg of Ryningsnäs, corresponding to a Coriolis parameter of p c = 1.22 × 10 −4 1/s.From preliminary simulations, it was found that |u g | = 16 m/s at φ = 66 deg would yield the desired main incoming flow direction of φ = 181 deg at hub height (given by the relative position of the metmast with respect to the rotor).This previous work also showed that large inertial oscillations occur in the upper part of the domain due to the introduction of the Coriolis force.These are seen as fluctuations in the velocity magnitude as well as 360 deg spins in the flow direction at the uppermost parts of the domain [12].To remove these, Rayleigh damping is applied following suggestions by [24].The damping follows a sinusoidal increase in the vertical direction and linear for the horizontal terms, starting at z asl = 1500 m (z agl = 1377.Next to damping, a preliminary simulation was needed to achieve a balance in the flow solution for the given pressure gradient.This was carried out in a domain of the same dimensions but using coarser grid due to the long period required to achieve the equilibrium in regions above the ABL height, where turbulence is essentially null.The horizontal cell size was 50 m for the innermost region and 100 m for the edges with the same vertical resolution as described above.For this simulation, the initial velocity is obtained from an Ekman profile, calculated for an eddy viscosity of 164 m 2 /s.The simulation is set as neutral, with a constant temperature profile of 282 K up to z agl = 800 m and a capping inversion from there upwards of 8 K/km.Random fluctuations are applied to the initial profiles, with an amplitude of 0.5 K for temperature and 0.5 m/s for the velocity components u x and u y .The preliminary simulation ran for 12.42 days, at the point where the flow was considered nearly devoid of oscillations around the geostrophic velocity magnitude and direction.This flow solution was mapped onto the domain with the finer mesh described above.The simulation was then run for spinup period of 8 h before starting the diurnal cycle. To produce the stratification changes along the day, a transient net radiative flux Q h (t) is imposed on the domain.This value follows an Gaussian-like curve based on measurements from Ryningsnäs obtained during diurnal periods when daylight lasted for about 13 h (see Sec. 2) as shown in Fig. 1-(c).The value of Q h is spatially uniform above the canopy top, producing no changes in temperature following eq.( 2).The same equation shows that the net radiation heats or cools the forest downwards from the canopy top, as a function of the light penetration and A c .The ground temperature is maintained at a constant 282 K, so a heat flux is produced at this boundary to maintain a temperature balance with the forest layer above it.Heat exchange from evaporation is ignored in the present framework.Since the domain surface is completely covered by forest, the source term eq. ( 2) is the only supply of variations in heat flux and hence, the only driver of the stratification changes along the diurnal cycle.
After the spinup time, the LES is run to compute a 24 h diurnal cycle starting at 5:00 h.A varying time-step that maintains a CFL number smaller than 0.65 is used, which in average corresponds to a ∆t ≃ 0.134 s (with a std of 0.012 s).

Results and discussion
In the first set of results, flow fields are time averaged in periods of 30 min and sampled in the z−direction at the locations p and at the metmast (Fig. 1-(b)) to capture the inflow conditions of the rotor as well as the wake.These are in turn averaged within day and night periods (analogously to the measurement data), chosen according to the stretches of positive and negative net radiation in Fig. 1-(c).The corresponding periods are 10:00 h to 16:00 h for day and 21:00 to 03:00 for night, identified by coloured shades in that figure.Next to these, the initial condition at 5:00 h, just before the start of the diurnal cycle and therefore referred to as neutral, is also included.All variables in this section represent resolved quantities, so the ' • ' is dropped to simplify notation.The profiles of Fig. 2-(a) show features typical of the unstable and stable regimes, such as the reduced or extended shear during the day or night, respectively.A Low-Level Jet (LLJ) is clearly visible during the night, peaking at around 500 m.At 5:00 h, a comparatively small LLJ is visible at around 1000 m.The velocity is appreciably slower during the day, but only above ∼ 153 m where it intersects with the night profile, below this height -including the rotor area-the velocity is higher during the day.The velocity in the neutral condition is slower than in the day below ∼ 300 m and slower than the night essentially for all the ABL extension.Fig. 2-(b) shows a closer look at these results, permitting to see the differences between the inflow and the wake flow for the different periods and how they compare with the metmast measurements.Note that LES results are normalized with the inflow velocity at hub height whereas the measurements use the inflow velocity estimated by the nacelle anemometer which employs a transfer function to correct for blockage.The LES points towards a slower wake recovery during the night as observed in the measurements (LES yields faster inflow during the day) although not as markedly.In comparison, the daytime wake deficit is very well predicted by the LES.Fig. 2-(c) shows a noticeable difference in the potential temperature between the day and night periods, both at rotor heights and within the forest.The temperature remains largely equal in the wake compared to the inflow.The turbulent heat flux θ ′ w ′ in Fig. 2-(d) reaches a larger magnitude during the day, due to the larger influx of positive net radiation, as well as spreading over higher regions (it reaches zero at ∼ 1700 m) compared to the night values.Notably, the wake has the effect of increasing the heat flux magnitude during both periods.Due to the constant temperature set at the ground, the heat flux curves change their slope within the forest.
The velocity decay along the wake can be seen in Fig. 3 (left) for the different stratification periods.The results support the hypothesis that wake recovery is faster in the unstable conditions during the day compared to neutral and nighttime.A turning of the wake can also be discerned.Whereas the neutral conditions display very little turning at all, in stable conditions the wake initially exits the turbine to the left of the incoming wind direction but then becomes centralized relative to the inflow direction further downstream.In unstable conditions the results show a very weak turning to the right.It is interesting to note that these results contrast those of [25,26] which both saw turns to the right for actuator disk simulations.Further investigation into the relative strengths of shear stress divergence and the Coriolis effect at different stratifications could help to clarify the situation regarding wake turning in forested regions.The change of temperature and turbulent heat flux along the day and at hub height is displayed in Fig. 3 (right), allowing a direct comparison with the net radiation Q h from Fig. 1-(c).While θ ′ w ′ is the largest at around 11:30 h, matching the maximum of Q h , the temperature peaks only until 14:30 h.It is worth noticing that the maximum/minima of Q h used in this work are appreciably higher/lower that what has been employed in other works that simulate a diurnal cycle without a forest [7,11,10].

Spectral densities
Consistent with the previous results, the spectral densities are presented as a daytime and a nighttime average.The power spectral densities of all three velocity components as well as that of the uw-cospectra are shown in Fig. 4. For the measurements curves, data are taken from the wake free western wind direction whereas the LES data were sampled over probes located over the arc of circle of radius 3D, extending from -40 to 40 deg from the position p in Fig. 1-(b) with probes located every 10 deg and at the height of the sonic anemometers.Overall, the agreement between the simulated and measured spectra is better during daytime than nighttime.As expected, the length scales of turbulence are smaller during nighttime (as indicated by the peak frequency of the curves).It is clear from the comparison that the differences between night and daytime turbulence are larger in the measurements than in the LES.This is attributed to a longer lengthscale during nighttime in the simulations than in the measurements, particularly for the vertical component (but also clearly noticeable in the lateral).One potential reason is the fairly high wind speed in the LES, leading to a high shear production even though the net radiation is negative.One may also observe that while the measurements show a broad frequency range following the f −5/3 decay slope of spectral density, the simulations clearly indicate an increased loss of energy due to unresolved turbulence as normalized frequencies exceed 0.3.The effect of the resolution is larger during nighttime because of the higher relative importance of the smaller turbulence scales.

Wake cross sections
To map the character of the wake from the metmast, the data was conditionally sampled for different incoming flow directions and combined to map the 2-dimensional cross section of the wake.The resulting wake structure is shown in Fig. 5  nighttime conditions.The wake is clearly more diffused during daytime and for both LES and measurements, the nighttime wake is characterized by an asymmetry consistent with a veering wind profile.
Even though the planes are sampled from the same spatial positions in both LES and measurements, there is an important distinction: The LES had a homogeneous forest cover whereas the measured wake is likely influenced by the forest edge which in reality was present along the line between the turbine and the metmast, missing in the simulations.The influence of the homogeneous forest cover may be noted by the elevated turbulence intensities at 40 m height relative to higher heights in the LES, appearing as yellow in the bottom of Fig. 5, a feature that is much less marked in the measurements.Another reason for discrepancy between the modelled and measured wake is the close distance between the wake cross section and the turbine, 2.23D.Since the turbine was simulated as a uniformly loaded actuator disk without tangential forces, the wake structure is not expected to compare very well in the near wake region.It is also clear from Fig. 5 that the structure of the wake is less diffused in the measurements relative to the LES.Other reasons for discrepancies could be random error and difference in stratification.Even though the model was forced with measured radiation balance, latent heat and inertial effects are expected to lead to the stratification being slightly different in the LES and the measurements.

Conclusion
An LES-based methodology has been presented to reproduce the continuous changes in atmospheric stratification along the diurnal cycle over forests and their influence on the turbulence characteristics of wind turbine wakes.The methodology employs a net heat flux radiation input as the main driver of the stability changes, heating and cooling a homogeneous forest that covers all the domain surface.An actuator disk is used to simulate the wake of a wind turbine located in Ryningsnäs, Sweden, for which metmast measurements were available to carry out a comparison.Results show that the methodology is capable of reproducing the main flow features characterizing night and day periods in the inflow as well as in the wake of the rotor.Spectral analysis shows that while the daytime distribution of turbulence energies are well predicted, the nighttime period displays distinct off-peak values in the LES compared to the measurements, indicating the usual challenge of reproducing the small scale turbulence of the stable periods.In general, results during the daytime period are better matched than the night by the LES, which could be due to an excess in the magnitude of the geostrophic velocity that drives the flow, producing an increased forest-induced shear that leads to larger length scales despite the negative heat flux, in comparison with the experimental conditions.Second order statistical momenta support this hypothesis, as well as comparing the wake structure that reveals the need to account for the heterogeneities of the forest and perhaps a more detailed modelling of the rotor.Both of these are tasks that will be addressed in upcoming studies.This work presents an integrated methodology to simulate the turbulence characteristics in the wake across the full diurnal cycle, which to the knowledge of the authors had not been attempted before over forested conditions at microscale levels.

Figure 1 .
Figure 1.(a) Horizontal view of the computational domain and its mesh regions; (b) relative positions of metmast MM, the AD and the sampling location p, additional sampling locations used for the statistical momenta calculations in Secs.4.1 and 4.2 are outlined over the arcs with centre in p and MM; (c) net radiative heat flux showing the measured and the prescribed value used in the LES, where the time periods chosen for night and day are shaded in red and blue, respectively.
3 m).For each term, the damping is proportional to c |u i − U ref |, with c = 0.1p c , U ref = 0 in the vertical direction and c = p c with U ref given by u g for the horizontal components.

Figure 2 .
Figure 2. Vertical profiles of the horizontal velocity (a)-(b), temperature (c) and turbulent heat flux (d) during the day and night periods and in the initial neutral period.Solid lines denote the inflow at p and dashed lines in the wake, at the metmast location.Small circle markers in the neutral curves are added at every cell centre to depict vertical mesh resolution.In (b) measurements are indicated by colour markers.The forested area is coloured in green and the vertical region behind the rotor is shaded in grey.

Figure 3 .
Figure 3. Left: Velocity deficit along the wake, sampled spanwisely at hub height, normalized with the mean inflow value.The grey region corresponds to the area behind the rotor.Right: from top to bottom, variation of temperature, turbulent heat flux and Obukhov lengthscale (calculated at hub height) along the diurnal cycle.Red and blue shades correspond to the day and night periods, respectively.

Figure 4 .
Figure 4. Normalized power spectra at 100 m height of the three velocity components as well as cospectra of uw.The left figure is the average for daytime conditions and the right figure for the nighttime conditions.Measurements are represented by the uniformly coloured curves and the LES results by the curves with black edges on the markers.

Figure 5 .
Figure 5. Turbulence intensity as measured in the tower during (a) day and (b) night and as simulated by the LES during (c) day and (d) night.