The Explicit Wake Parametrisation V1.0: a wind farm parametrisation in the mesoscale model WRF

We describe the theoretical basis, implementation and validation of a new parametrisation that accounts for the e ﬀ ect of large o ﬀ shore wind farms on the atmosphere and can be used in mesoscale and large-scale atmospheric models. This new parametrisation, referred to as the Explicit Wake Parametrisation (EWP), uses classical wake theory to 5 describe the unresolved wake expansion. The EWP scheme is validated against ﬁl-tered in situ measurements from two meteorological masts situated a few kilometres away from the Danish o ﬀ shore wind farm Horns Rev I. The simulated velocity deﬁcit in the wake of the wind farm compares well to that observed in the measurements and the velocity proﬁle is qualitatively similar to that simulated with large eddy simulation mod- 10 els and from wind tunnel studies. At the same time, the validation process highlights the challenges in verifying such models with real observations. a slower modelled recovery compared to that measured, whereas a positive slope shows a faster modelled recovery. There is no vertical resolution dependency for a wind farm parametrisation when the or triangles for a given wind speed bin lie 20 on top of each other. The results show that the bias in velocity between the measurements and the EWP simulations is small ( < 0.15 ms − 1 ) for all wind speeds ( except for the 7 m s − 1 bin at M7). We ﬁnd a positive velocity di ﬀ erence at M6 and a negative one at M7. Hence, the modelled wake recovery oscillates between being slightly slower from M6 to M7 and 25 being slightly faster from the end of the wind farm to M6, as well as from M7 to the end of the wake. The velocity of the EWP scheme at M6 is nearly resolution independent and at M7 the dependency is very weak.


Introduction
Wind turbines capture the kinetic energy of the wind with their turning blades, which transfer the energy to a transmission system that drives an electric generator. In this 15 process the flow in front and behind a wind turbine is decelerated by the forces acting on the rotating blades. In large wind farms, the interaction of the flow and the wind turbines is further complicated by the interaction of the wake of one wind turbine with neighbouring turbines. Besides the changed velocity field around the turbines, there is also evidence that wind turbines affect planetary boundary layer (PBL) processes due We implemented the EWP scheme in the open source Weather Research and Forecast (WRF) model (Skamarock et al., 2008). The WRF model already includes a wind farm parametrisation option, WRF-WF (Fitch et al., 2012;Jiménez et al., 2014). We validated the WRF-WF and EWP parametrisations against long-term meteorological (met) mast measurements in the wake of an offshore wind farm. To our knowledge, 15 measurements in the wake of a wind farm have not been used for the validation of a wind farm parametrisation model in previous literature.
In Sect. 2, we explain the theoretical basis of the EWP scheme and its implementation in the WRF model. Section 3 describes the measurements, whereas Sect. 4 introduces the WRF model set-up and the configuration of the EWP and WRF-WF 20 scheme. In Sect. 5 both wind farm parametrisations are validated against the met mast measurements in the wake of the wind farm. A discussion of the results finalises the article in Sect. 6, followed by the conclusions in Sect. 7.

The Explicit Wake Parametrisation
We start by introducing the relevant mesoscale model equations. Thereafter, the addi- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | added to the model equations. At the end of the section their implementation in the mesoscale model is described.

The mesoscale model framework
Wind turbines are well described by drag devices that slow down the wind velocity from a free-stream value U. Downstream, due to mixing of fluid particles from within 5 and outside the wake, the velocity deficit is gradually reduced to the point at which the background conditions are restored.
We use a mesoscale model for the simulation of the wind farm wake and its recovery. It uses the Reynolds Averaged Navier-Stokes (RANS) equations, 10 to describe the flow evolution. We will use the notation and symbols of Wyngaard (2010). The upper-case letters refer to ensemble-averaged quantities, whereas lowercase letters refer to fluctuations. The exception is the average air density ρ(x, t), where t denotes the time and x the position. In Eq. (1), U i (x, t) and P (x, t) represent the mean velocity components and the pressure, whereas f and g are the Coriolis frequency and 15 the gravitational acceleration constant. Furthermore, the right most term F D i (x, t) is the ensemble-averaged horizontal forcing due to the action of wind turbines (F D 3 = 0). Mesoscale models generally simulate the effects of turbulence in the vertical direction only. The components of the Reynolds stress are parametrised in a 1.5-order PBL scheme as where the turbulence diffusion coefficient for momentum, K m (x, t) = S m √ 2 e, depends on the stability function S m (x, t), the turbulence length scale (x, t) unit mass e = u i u i /2. We write the most general form of the TKE equation: where on the left hand side (l.h.s.) E denotes the TKE and T the transport, which includes the advection by the mean flow, turbulence transport of TKE and the divergence of the pressure correlation. On the right hand side (r.h.s.), P s represents the turbulence 5 production from the vertical shear in the horizontal velocity (shear production), P b the turbulence production or destruction related to buoyancy forces, P t the turbulence induced by the turbine rotor, and the dissipation.

The mesoscale model grid
The previously defined variables in the Eqs. (1) and (3) have to be redefined on the 10 three-dimensional model grid. A new volume-averaged velocity equation is derived by integrating Eq. (1) over the grid-cell volume. For Eq.
(3), we define a (random) velocity fluctuation at any hypothetical measurement point within a grid-cell volume to be the difference between the ensemble-averaged and the instantaneous velocity. This is illustrated schematically in Fig. 1a. Then, the grid-cell averaged TKE can be thought of 15 as the variance of the velocity fluctuations around the ensemble-averaged velocity and not around the grid-cell averaged velocity as shown in Fig. 1b. The aim is to obtain expressions for the terms P t (Eq. 3) and F D (Eq. 1) that are consistent with the mesoscale model flow equations. The expression for P t can be found by multiplying the Navier-Stokes equations with the velocity fluctuation and then 20 applying Reynolds averaging. This gives for the additional source term: where ρ A r c T (U i + u i ) 2 is the instantaneous forcing due to the action of the wind turbine, the lower case c T the instantaneous thrust coefficient, the upper case C T the 3486 Introduction averaged thrust coefficient, and A r the rotor area. This additional turbulence occurs directly behind the turbine blades. Its turbulence length scale is expected to scale with the blade's cord length, which is on the order of a few metres and it is therefore much smaller than that of the atmospheric flow. The smaller length scale implies a significant higher dissipation rate ∝ (1/ )u i u i 3/2 . Therefore, we assume the source term 5 in Eq. (4) to dissipate within a mesoscale model grid-cell and neglect P t of Eq. (4) on the grid-cell average. Additional turbulence is generated by shear production, which we assume to be the dominant mechanism on the grid-cell average. In the following sections, we derive the expression for the grid-cell averaged F D .

10
The velocity deficit expansion within one grid-cell is not negligible, which leads to flow decelerations that extend beyond the rotor swept area. This part of the wake expansion is not accounted for in the mesoscale model, hence we estimate it explicitly with a sub-grid scale (turbulence) diffusion equation. Then, a grid-cell averaged force is determined and added to the model RANS equations.

15
The sub-grid-scale model describes the unresolved turbulence diffusion process that results from turbulence shear production. In this model, we assume the horizontal advection of velocity and the turbulence diffusion to dominate in Eq. (1). Considering first only the flow behind the turbine rotor, we obtain from Eqs. (1) and (2) the diffusion equation: that describes the expansion of the velocity deficit, U 0 − U, behind the turbine. Here, we denote the advection velocity at hub-height h by U 0 = |U(h, t)| and the unresolved velocity in the stream-wise direction x in the wake of the turbine by U(x). The turbu- lence diffusion, which causes the wake expansion, is described by a single turbulence diffusion coefficient K = K m (h, t).

The velocity deficit profile
We define the vertical structure of the velocity deficit as U d = U s ξ, where U s (x) is the maximum velocity deficit at the centre of the wake and ξ(x, y, z) a function that deter-5 mines the vertical wake extension (Tennekes and Lumley, 1972). Equation (5) can be solved for the velocity deficit profile: where the length scale that determines the vertical wake extension is 10 In Eq. (7), σ 0 denotes the initial length scale. Equation (7) represents the vertical wake extension, resulting from turbulent diffusion of momentum and it is similar to the solution of Eq. (4.29) in Wyngaard (2010) for the dispersion of plumes. Equation (6) describes the ensemble-averaged profile of the velocity deficit around hub-height at a given point in the far turbine wake (Tennekes and Lumley, 1972). The turbine wake is normally 15 divided in a near and far wake, where the far wake begins between one rotor diameter D 0 (Vermeer et al., 2003) and 3 D 0 (Crespo and Hernández, 1996). We can find the velocity deficit profile for wind turbines by equating the total thrust to the momentum removed by the action of the wind turbine, i.e., where C T (U 0 ) is the turbine thrust coefficient and R 0 the radius of the rotor. In Eq. (8), the l.h.s. represents the local forces at the rotor swept area, whereas the r.h.s. describes the equivalent distributed force for the expanded wake at any x. From the Eqs. (8) and (6), we find the velocity deficit profile for a specific thrust force, When we insert the velocity deficit of Eq. (9) in the second term of Eq. (8) and integrate in the cross-stream direction y, we have the integrated thrust profile The term on the r.h.s. within the integral describes the equivalent distributed thrust force in the vertical direction at any distance x. Next, we derive from Eq. (10) a single 10 effective thrust force, which represents the average wake expansion within a grid-cell.

Turbine forcing in the mesoscale model
For the mesoscale model we derive an effective thrust force, which describes the average wake expansion within a grid-cell. Therefore, we first determine the effective velocity deficit profile U e by averaging the velocity deficit of Eq. (9) over the cross-15 stream direction y and over a downstream distance L that the wake travelled within the grid-cell. It is convenient to approximate this area-averaged velocity deficit profile by a Gaussian-shaped profile: Here σ e is the effective length scale that is related to the model grid-size, From the definition of the effective length scale, we can obtain the total effective thrust force F e by substituting the length scale σ in Eq. (10), with the effective length scale σ e . This gives The grid-cell averaged acceleration for the model RANS equations, Eq. (1), is now obtained, when we divide Eq. (13) by the mass and apply the grid-cell volume. For every vertical model layer k, we obtain then from Eq. (13) the grid-cell averaged acceleration components, and in the x and y direction, respectively. In Eqs. (14) and (15), we use ∆ x and ∆ y for the horizontal grid-spacing in the x and y direction, ϕ(k) for the wind direction, and the angular brackets for the spatial averaging. For the height z(k), we use the height at the centre of layer k.

Implementation in the WRF model
We assume every turbine within a grid-cell to be positioned at its centre and use L = ∆ x/2 in Eq. (12). For each turbine a thrust force is calculated with Eqs. (14) and (15). The total thrust force for a turbine containing grid-cell is then obtained from a superposition of the single turbine thrust forces and is added to the numerical approxi-5 mation of Eq. (1). We use the grid-cell averaged velocity as the upstream velocity U 0 to the wind turbine (see Sect. 6). The turbulence diffusion coefficient for momentum, K in Eq. (12), comes from the selected PBL scheme. The parametrisation, therefore, can be used with any PBL scheme. However, to be consistent with the assumptions used in the derivation 10 of the wind farm parametrisation, a second order closure scheme with a turbulence shear production term is recommended. A practical description of how to use the EWP scheme in the WRF model is given in Sect. . 15 The parallelogram-shaped Horns Rev I wind farm is situated in the North Sea about 15 km West of the Danish coast (Fig. 2). It is made up of 80 Vestas V80 (2 MW) pitch controlled, variable speed turbines, resulting in a total rated wind farm capacity of 160 MW. The turbines have a swept area diameter of 80 m with the hub mounted at 70 m a.m.s.l. The equally-spaced turbines are placed in 10 columns from West to East 20 and 8 rows from North to South, labelled C1-C10 and R1-R8 in Fig. 2b. The spacing between columns and rows is 560 m, which is equivalent to 7 turbine diameters. Hansen et al. (2012) describes in detail the various production data available at the wind farm, their quality filtering, and processing. Turbine (C1,R7) is used as a reference turbine, since it is not affected by the wind farm wake under westerly flows; in Fig. 2b it is marked with the solid bullet. The wind speed for each turbine is estimated from the power production data, using the turbine specific power curve and the wind direction is obtained from the adjusted yaw angle (Hansen et al., 2012).

GMDD
For the validation, we use data from two met masts, M6 and M7, whose positions are shown in Fig. 2b. The masts are located 2 and 6 km to the east of the eastern edge of the wind farm and are thus directly in its wake for westerly winds. The 70 m tall masts are identically equipped and their instrumentation includes high quality Risø cup anemometers for measuring wind speeds. The data from the two masts, which 10 are independent of the wind farm data, is used to validate the modelled wind speed reductions downstream from the wind farm.

Data selection and averaging for the validation
For the validation of the results of the model simulations, and especially because of the relatively large area-averaged fields in the mesoscale model, it is very important to 15 ensure that measurements and model fields are comparable. This can be achieved by properly selecting the wind direction and wind speed interval at the reference turbine.
Regarding the wind direction selection, Gaumond et al. (2013) found at Horns Rev I, for narrow wind direction sectors, a higher measured turbine production compared to that estimated by wake models. Most likely this is related to wind direction variations, 20 which, for small wind direction bins, exceed the bin size within the 10 min averaging period. Since these variations are not accounted for in the model, we follow the recommendation to average over a relatively large wind direction interval. Furthermore, it is important that the flow reaching the mast anemometers has passed through the wind farm, otherwise it does not characterise the wind farm wake. Therefore, we se- In the period 2005 to 2009, we select wind speeds from 6.5 to 11.5 m s −1 at the reference turbine and bin them in 1 m s −1 intervals. In this range, we guarantee to be above the cut-in wind speed of the turbines (4 m s −1 ) and below the wind speed at which the control system starts to pitch the blades around (13 m s −1 ). To obtain as much data as possible, we do not filter the measurements on stability. For strong westerly winds, 5 we expect the lower atmospheric stability on average to be neutral (Sathe et al., 2011). For every instance that the wind direction and wind speed at the reference turbine was within the above described range, the wind speeds at M6 and M7 were accepted. Afterwards, the selected wind speeds were averaged for each wind speed bin and normalised by the average wind speed at the reference turbine U ref (Table 1).

WRF model configuration and averaging
In the simulations we used the WRF model V3.4. However, the WRF-WF parametrisation has been updated to the version available in WRFV3.6. The model domain is setup with 80 × 40 grid-cells, with a horizontal grid-spacing ∆ x = ∆ y = 1120 m. Equally to the Horns Rev I wind farm, the model wind farm contains 80 V80 turbines and extends 15 over 5 grid-cells in the West-East and North-South direction (Fig. 3a). The Vestas V80 thrust and power curves were used for the turbine parameters.
The model was run in an idealised case mode with open lateral boundary conditions (Skamarock et al., 2008, p. 51), Coriolis forcing and zero heat fluxes from the lower boundary. At the surface, the no slip condition holds. The surface roughness was con-20 stant in time and set to z 0 = 2 × 10 −4 m for the entire domain. The model simulations were run with the MYNN 1.5 order Level 2.5 PBL scheme (Nakanishi and Niino, 2009), on which the WRF-WF scheme is dependent. These and other details of the model configuration are summarised in Table 2.
We ran simulations for 5 wind speed intervals and for 9 wind directions, ranging from Introduction With these vertical resolutions there are 5, 7, and 10 grid-cell volumes intersecting with the turbine rotor as shown in Fig. 3b. All 135 simulations (5 wind speeds, 9 wind directions, and 3 vertical resolutions) were initialised with a constant geostrophic wind in height in a dry and slightly stable atmosphere. After a four day integration period, the wind converged in the whole domain to 5 a logarithmic neutral profile within a 650 m deep boundary layer and remained independent of height above the inversion layer. The atmospheric state of each of these 135 simulations was used to drive: a control simulation without wind farm parametrisation, a simulation with the WRF-WF scheme, and a simulation with the EWP scheme. Each control or wind farm simulation lasts one day, resulting in a total simulation length of 10 five days. The wind speeds in the control simulations were 7, 8, 9, 10, and 11 m s −1 at 70 m (hub-height) after 5 days of simulation.
To summarise, we performed simulations for 9 wind directions, 5 wind speeds, and 3 vertical resolutions, with and without wind farm parametrisation, which gives a total of 405 (3 times 135) simulations as outlined in Table 3. For the validation against the 15 mast measurements, we averaged the model wind speeds over the 9 wind directions. We used the instantaneous wind speeds at the end of the simulation period. Similarly to the observations, we normalise the wind direction averaged wind speeds and hereafter simply refer to them as "velocity". The model wind speeds from the simulations without the wind farm were used to normalise the wind farm flow. Introduction For this calibration, we selected turbine production data at time stamps when the derived upstream wind speed and wind direction at the reference turbine ranged from 8.5 to 9.5 m s −1 and from 255 to 285 • . This wind speed bin was selected, since it contains the most turbine observations. It has a minimum of 850 observations for a turbine at the eastern edge of the wind farm and a maximum of 1612 observations for a turbine 5 in the first wind farm column. The difference in the selected number of observations results from the additional requirement that for a given turbine the local upstream turbines have to be operational, in order to guarantee wake losses. For the power production data obtained in this way, we used the power curve to derive the wind speed. The column averaged velocity was afterwards obtained by averaging over the inner 6 10 rows (R2-R7 for all turbines from σ 0 = R 0 to σ 0 = 2 R 0 , stepping with ∆R 0 = 0.1, to determine the best fitting value. For the comparison to the measurements, the column-averaged wind speed was obtained by averaging over the 3 central wind farm grid-cells in the crossstream direction. Figure 4 shows the results from the simulations with σ 0 = 1.7 R 0 , which had the small-20 est overall bias compared to the measurements. It shows the same velocity reduction for all 3 vertical resolutions. Consequently, the amount of energy extracted by the turbines is independent of the vertical resolution. The simulations with an initial length scale σ 0 = 1.7R 0 follow the measured velocity reduction fairly well. Therefore, we conclude that for neutral conditions the initial length scale can assumed to be independent Introduction

The WRF-WF scheme
In the validation, we also include results from the wind farm parametrisation WRF-WF from WRF-V3.6. This parametrisation was introduced in WRFV3.4 by Fitch et al. (2012). In this approach, a forcing term F D is applied to every model level that intersects with the wind turbine blade sweep area and thus the scheme simulates the local drag 5 forces over the turbine rotor. An additional TKE source term is parametrised as here C T and C P are the thrust and power coefficients of the turbine, whereas U is the grid-cell velocity. The application of one-dimensional theory (Hansen, 2003) on Eq. (16), gives P t,WRF-WF = ρ A r C T a U 3 /2 for the additional TKE term. Here, a denotes 10 the induction factor. The same result, P t,WRF-WF = ρ A r C T a U 3 /2, can also be obtained by defining a velocity fluctuation as the difference between the grid-cell averaged velocity and the instantaneous velocity as illustrated in Fig. 1b. This has been analytically derived by Abkar and Porté-Agel (2015) (their Eq. 21, with their ξ = 1). In this definition of TKE, 15 the turbine induced velocity reduction is counted as TKE. The TKE source term from the WRF-WF scheme in Eq. (16) is much larger than the source term in Eq. (4), which comes from the different definition of TKE in the two schemes. In WRFV3.4, the WRF-WF parametrisation the power coefficient has been obtained from the power curve information. The thrust coefficient was then determined by an To obtain a complete picture of the modelled velocity field within and downstream of the Horns Rev I wind farm, we compare the hub-height velocity simulated by the WRF model using the two wind farm parametrisations for one wind speed bin. We used the second most frequently observed wind speed bin (10 m s −1 ), which is different from that 5 used in the calibration (Sect. 4.1.1), to be as calibration independent as possible. Afterwards, we compare the modelled velocities for all wind speed to the measurements at M6 (2 km downstream) and M7 (6 km downstream). In a qualitative validation, we examine the simulated wind farm wake characteristics, as well as the vertical profile of the wake deficits, using results from LES simulations and wind tunnel experiments as 10 a reference.
In the validation, we use the instantaneous model outputs after the 5 days integration period. Furthermore, the velocities are normalised as described in the Sects. 3 and 4. Figure 5 shows the wind direction averaged hub-height velocity (10 m s −1 bin) as a func-15 tion of the downstream distance for the EWP and WRF-WF schemes and all vertical resolutions. For the EWP scheme ( Fig. 5a) there is no identifiable vertical resolution dependency on the velocity, while for the WRF-WF (Fig. 5b) this is small. In the EWP scheme the velocity within the wind farm decreases almost linearly with distance, whereas for the WRF-WF scheme, after a more rapid initial decrease, the velocity be- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | indicated the velocity from the EWP scheme at the most easterly grid-cell of the wind farm with a dotted horizontal line in Fig. 5. This agreement is noteworthy, given the two different methods used. However, the velocity 2 km downstream of the wind farm, at M6, differs between the schemes by 4.7 %. This means that downstream of the wind farm grid-cells, where the wind farm parametrisations are not active, the velocity di-5 verges for the two schemes and the difference becomes significant. This is important, especially if it was used to estimate the power production on a neighbouring wind farm that was located at this distance from the original wind farm. For example, the Rødsand 2 and Nysted offshore wind farm in Southern Denmark are separated by a comparable distance.

Velocity recovery at turbine hub-height
10 Figure 6 depicts wind direction averaged velocity recovery at the two met masts for all 5 wind speed bins and all vertical resolutions. It shows the differences in the WRFmodelled and measured wind speed at M6 (circles) and M7 (triangles), against the measured wind speed at the same masts. The WRF-modelled wind speed is obtained from linear interpolation of the wind speeds in the nearest grid-cells. The modelled 15 recovery rate can be deduced from the slopes of the solid (between the values at M6 and M7) and dashed lines (between M7 and the free-stream velocity). A negative slope is linked to a slower modelled recovery compared to that measured, whereas a positive slope shows a faster modelled recovery. There is no vertical resolution dependency for a wind farm parametrisation when the circles or triangles for a given wind speed bin lie 20 on top of each other.
The results show that the bias in velocity between the measurements and the EWP simulations is small (< 0.15 m s −1 ) for all wind speeds ( except for the 7 m s −1 bin at M7). We find a positive velocity difference at M6 and a negative one at M7. Hence, the modelled wake recovery oscillates between being slightly slower from M6 to M7 and 25 being slightly faster from the end of the wind farm to M6, as well as from M7 to the end of the wake. The velocity of the EWP scheme at M6 is nearly resolution independent and at M7 the dependency is very weak. The WRF-WF scheme shows a positive bias in velocity of up-to 0.5 m s −1 at M6. This positive bias, between the WRF-WF scheme and the measurements, becomes larger with increasing wind speed. The bias at M6 is a consequence of the too fast modelled wake recovery from the end of the wind farm to M6. Between M6 and the point at which the free-stream velocity is reached again, the modelled wake recovery is slower than 5 that measured. This overall positive bias implies that the modelled velocity with the WRF-WF scheme is overestimated throughout the whole wake. Figure 7 shows the spatial structure of the modelled velocity (10 m s −1 bin) within and in the wake of the wind farm for the L60 simulations in the 270 • wind direction. A comparison of the results from the two schemes confirms the faster initial wind farm 10 wake recovery in the WRF-WF scheme. The 10 % velocity deficit contour, for example, extends for the EWP scheme to around x = 15 km (Fig. 7a), while for the WRF-WF scheme to x = 8 km (Fig. 7b). The difference in the distance at which a 7.5 % velocity deficit is reached is even larger: after x = 21 and x = 11 km for the EWP and WRF-WF scheme, respectively. Further downstream, after around 30 km, the velocity fields 15 from the two parametrisations become similar. Possible reasons for the difference in the initial wake recovery are discussed in the next section.

GMDD
Finally, Fig. 7 shows a difference in the orientation of the axis of the velocity deficit downstream from the wind farm. In both simulations the steady-state wind direction of the free-stream flow was 270 • . As the velocity decreases within the wind farm, the 20 velocity is expected to turn to the North, due to the changed Coriolis force. In the wind farm wake, where the flow accelerates again, the velocity should turn back again to the background direction. This effect is seen for the EWP scheme (Fig. 7a). Whereas, for the WRF-WF scheme (Fig. 7b) the wake turns southward behind the wind farm. A possible reason for this unexpected behaviour is the turbulence transport of momentum 25 from aloft (Ekman spiral) within the wind farm. In the wind farm wake the flow keeps turning to the South because of the flow acceleration from the momentum transport.

Vertical profile of TKE and velocity
To obtain a broader understanding of the mechanisms acting in the two schemes, we compare the simulated TKE (per unit mass) and the velocity deficit profiles for the 10 m s −1 wind speed bin in the 270 • wind direction. Figure 8 shows the difference in TKE (wind farm minus control simulation) for the 5 L60 simulation. The cross-sections in the x, z plane are in the West-East direction and pass through the centre of the wind farm. We used different colour scales in the two plots, due to the relatively large differences in TKE production between the two schemes. However, we have kept the outermost contour (0.02 m 2 s −2 ) the same. The maximum TKE difference was 0.30 and 1.9 m 2 s −2 for the EWP and WRF-WF scheme, 10 respectively. Compared to the environmental TKE of the pure shear flow (not shown), the TKE increases in the EWP scheme (Fig. 8a) at the end of the wind farm by a factor of 2, whereas in the WRF-WF scheme (Fig. 8b) it increases at hub-height by a factor of 5.5. The EWP scheme shows an increased and decreased TKE above and below hub-15 height compared to the reference simulation. The maximum increase occurs behind the wind farm, where the velocity gradients are the largest. In the WRF-WF scheme, the maximum TKE increase happens at hub-height within the wind farm. Recalling that in the WRF-WF scheme, turbulence is generated by the source term P t and by the PBL scheme from turbulence shear production (P s ), the intensity of P t dominates over 20 that of the shear production. The sum of the additional source term and the turbulence shear production, causes within the wind farms an increased turbulence from the lowest model level upwards.
We use the results from the actuator disk approach without rotation from Wu and Porté-Agel (2013), to obtain a qualitative impression of the structure of the turbulence 25 field from a LES model within a wind farm. The actuator disk approach from their LES model is most similar to the drag approach in the mesoscale model. Their Fig. 5c shows that higher and lower turbulence intensities dominate around the upper and lower tur- bine blade tip. Also, a positive and negative shear stress occur above and below hubheight (their Fig. 7c). This indicates that the shear in velocity dominates the turbulence production. Similar features are present in the TKE field from the EWP scheme, where the increased and decreased TKE above and below hub-height are in like manner caused by an enhanced and reduced turbulence shear production with respect to the 5 neutral logarithmic velocity profile. Furthermore, Wu and Porté-Agel (2013) shows an upper wake edge at around 4.5 turbine hub-heights for 10 aligned wind turbines after 60 turbine diameters (their Fig. 12). They defined the wake edge at the point where the velocity reduction for a given height was 1 %. Similarly, the edge of the wake can be found from an increased TKE due to shear production. For the outermost contour 10 in Fig. 8, we obtain a vertical wake extension of around 5 turbine diameters for the EWP scheme at 4.8 km downstream (equivalent to 60 turbine diameters). At the same distance, it is around 7 turbine hub-heights for the WRF-WF scheme. The influence of the TKE field to the velocity profile is analysed in Fig. 9. It shows the velocity deficit profile ∆U x /U 0 h for both schemes and all vertical resolutions. Here 15 U 0 h denotes the free-stream velocity at hub-height. The velocity deficit is defined as ∆U x (z) = U(z) − U free (z), where U free (z) is the free-stream velocity profile from the reference simulation and U(z) the velocity profile from the wind farm simulation. We choose the second and third grid-cell within the wind farm (C2 and C3) and the second point behind the wind farm (C7) that corresponds approximately to the location of M6. 20 Figure 9a shows the velocity deficit profiles from the EWP and WRF-WF scheme from the L60 simulation. The profiles indicate a stronger diffusion in the WRF-WF scheme compared to that in the EWP scheme. This can be recognised by the vertical extension of the vertical deficit profile at C7, behind the wind farm.
For the EWP scheme (Fig 9b), we find a maximum velocity deficit at hub-height and 25 symmetric features around the maximum for the L40 and L80 simulations. Also, results from LES simulations, wind tunnel experiments, and measurements (Vermeer et al., 2003;Wu and Porté-Agel, 2013;Iungo et al., 2013) have a maximum velocity deficit at hub-height in the far turbine wake for neutral conditions. The profiles in Vermeer et al.  . 36) and in Wu and Porté-Agel (2013) (their Fig. 4) additionally show symmetric features around the maximum with a shape similar to the velocity profile of the EWP scheme. Figure 9b demonstrates the vertical resolution independence of the EWP scheme within the wind farm: the velocity deficits of the L40 and L80 simulations are almost identical. A small difference is found below hub-height in the wind farm 5 wake.

GMDD
With the WRF-WF scheme (Fig. 9c) the maximum velocity deficit is displaced vertically above hub-height, which reaches the upper wind turbine blade tip at mast M6 (C7). The dependency of WRF-WF scheme to the chosen vertical resolutions is weak; differences are found within the wind farm (C2, C3) above hub-height. Also, the profiles from the WRF-WF scheme show increased (area-averaged) velocities inside the wind farm at the lowest model level. The wind farm simulations from Wu and Porté-Agel (2013) (their Fig. 13) do not support this behaviour.

Discussion
We use wind farm parametrisations implemented in a mesoscale models to simulate 15 the effect of wind farm wakes in areas on the order of hundreds of kilometres. However, the models have a limited horizontal resolution and hence the local processes within a wind farm remain unresolved. Wu and Porté-Agel (2013) found in their LES model results a sensitivity of the velocity reduction to the wind farm layout. In current implementations of wind farm parametrisations, all turbines within a grid-cell experience the 20 same upstream velocity. Although these parametrisations are not meant to estimate the local velocity field within the wind farm, differences in velocity reduction within the wind farm could affect the velocity in the wake of it.
In future implementations it can be possible to account for the local flow within the wind farm by using data from high resolution models, which can be input to the 25 mesoscale model via a look-up table as suggested by Badger et al. (2013) and Abkar and Porté-Agel (2015). However, we have at the moment no measurement data in the wind farm wake to validate these approaches for different wind farm layouts. At the Horns Rev I wind farm, we were restricted to the geometry of met masts positions, which did not allow to study the sensitivity of the velocity field in the wind farm wake for additional wind direction sectors. A fair comparison between the mesoscale model and long-term measurements can 5 be realised in several ways. One method is to match the simulation period to that observed. By selecting corresponding time frames, one can then compare the fields of interest. The main disadvantage of this method is that errors in the background flow simulated by the mesoscale model also exist, which complicates the analysis. The second method is to sample the data and the model simulations in rather strict idealised 10 conditions of wind speed and direction. We have chosen this second method using the WRF model in idealised case mode and compare its results to properly averaged measurements. Here the equilibrium solution without the wind farm effects is purposely made to match the free stream velocity and thus background errors in the simulations are absent. Besides a more detailed analyses, this method also allows the study of 15 the vertical dependency of the velocity reductions, since the atmospheric background conditions are almost identical for the simulations with the different vertical resolution. Before we validated the results of the schemes in the wake of the wind farm, the a priori unknown initial length scale of the EWP scheme had to be determined. We did this using the turbine power measurements from the most frequently observed wind 20 speed bin. This limitation could not be avoided, since to our knowledge no other longterm measurements from large offshore wind farms are available. This constant was then used for the validation of all wind speed bins. The difference of 0.1 % between the velocity deficit simulated by the WRF-WF and EWP scheme at the end of the wind farm for the 10 m s −1 wind speed bin, facilitated the comparison between the schemes 25 in the wind farm wake, where the parametrisations are not active anymore. One major difference between the EWP and WRF-WF approach is in the treatment of the grid-cell averaged TKE budget equation. The TKE production regulates the vertical profiles of momentum, temperature, and moisture within the PBL. Differences in

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | TKE production would thus affect the local weather (e.g., temperature, humidity, and possibly clouds) response to the presence of large wind farms. For the WRF-WF and EWP scheme, the PBL scheme adjusts the TKE production through the changed vertical shear in horizontal velocity. The term that represents the local turbine-induced turbulence is neglected in the EWP approach, whereas it is parametrised as an addi-5 tional source term in the WRF-WF scheme. The total TKE is more than 3 times larger in the WRF-WF scheme than in the EWP scheme. Future measurement campaigns around large wind farms under suitable atmospheric conditions can help to settle this issue.

10
We introduce a wind farm parametrisation for use in mesoscale models. The EWP approach is based on classical wake theory and parametrises the unresolved expansion of the turbine-induced wake explicitly in the grid-cell that contains turbines, where the largest velocity gradients occur. The associated turbulence shear production is then determined by the PBL scheme in the mesoscale model. The approach has been im- 15 plemented in the WRF mesoscale model and can be used with any PBL scheme. However, we recommend PBL schemes that model the TKE equation.
We analysed the results of simulations from the scheme in the wake of a wind farm. For the validation, we used the averaged wind speeds within a 30 • wind direction sector at two met masts in the wake of the Horns Rev I wind farm. The model was run 20 for several wind direction bins that cover those sampled in the observations. For each wind speed bin, we compared a wind direction averaged wind speed to the similarly averaged measurements. We found that for all 5 velocity bins, the velocity modelled with the EWP scheme agreed well with the met mast measurements 2 and 6 km downstream from the edge of the wind farm. The bias was less than 0.15 m s −1 , except for Introduction To our knowledge, no long-term data-sets are available to validate the details of the vertical structure of the velocity deficit and turbulence in the wake of the wind farm. In a qualitative comparison, we found the vertical structure of the modelled TKE field to agreed with that of actuator-disk simulations by LES models (Wu and Porté-Agel, 2013), with an increased and decreased TKE around the upper and bottom rotor tip, 5 respectively. The TKE field in the EWP scheme led to a symmetric velocity deficit profile around hub-height, similar to velocity deficit profiles in Wu and Porté-Agel (2013);Vermeer et al. (2003). Also, the vertical wind farm wake expansion in the EWP approach was similar to that described in the before mentioned studies. While validation of the WRF-WF parametrisation has before been carried out with measurements within 10 a wind farm (Jiménez et al., 2014), this is the first time that the validation has been done in the wake of a wind farm, at the scales the mesoscale model was designed to simulate.

Appendix
In the EWP scheme, we approximate the velocity deficit profiles on the r.h.s. of Eq. (11), 15 by a Gaussian shaped velocity profile on the l.h.s. of Eq. (11).
To show that these profiles are to a good approximation similar, we compare the difference between the average of 5000 single profile at distances 0 < x < 500 m to the approximated Gaussian velocity deficit. We normalise both profiles by the depth, ∆y, of the wake in the cross-stream direction. For the comparison we used U 0 = 8 m s −1 , 20 R 0 = 40 m, σ 0 = 60 m, C T = 0.8, K = 6 m 2 s −1 , and ∆ y = 1120 m. The wake centre is defined at z = 0 m.
The result in Fig. 10 shows that the differences between the spatial averaged Gaussian profiles and the Gaussian profile with the space averaged spread is in the entire velocity deficit region far less than 0.001 m s −1 . Introduction

Code availability
In this section a short guideline of the usage of the EWP scheme in the WRF model is given. The EWP approach can be run either in serial, shared-memory, or distributed memory options. Currently, it is not possible to run the approach with the mixed shared and distributed memory option. The scheme can be used for idealised and 5 real case simulations. For the real case simulations, the wind farm parametrisation can be activated in any nest of the simulations. The additional namelist.input option bl_turbine should be used to select the wind farm parametrisation. The EWP scheme needs additional input files in ASCII format. In the first file the positions and types of all wind turbines should be listed. The file name has to be specified as a string in the additional namelist.input option windturbines_spec. For every turbine used, the turbine characteristics, i.e., the hub-height and diameter, and the thrust and power coefficients as a function of wind speed, need to be listed in a file. This file name has to start with the turbine type used in the first file, followed by the extension .turbine. 15 Please contact pvol@dtu.dk to obtain the code of the EWP wind farm parametrisation.