Wind turbine wake models developed at the Technical University of Denmark: A review

Wind turbine wakes are one of the most important aspects in wind power meteorology because they decrease the power production and increase the loading of downstream wind turbines. Therefore, there is a continuous need to ﬁ nd a ‘ good ’ wake model to properly plan wind power plant-level control strategies, predict the performance and understand the fatigue loads of turbines. In this paper, six widely used approaches of wake modelling (Jensen, Larsen, Dynamic Wake Meandering, Fuga and, Ellipsys3D LES and RANS together with their interpretations) that were developed at Technical University of Denmark, are described and the model subcomponents are analysed. The models are evaluated using data from the Sexbierum (onshore) and the Lillgrund (offshore) wind farms to understand how to best utilize them. The paper provides a comprehensive conceptual background to wake modelling combined with the overview of the state-of-the-art models including their implementations on operating wind farms.


Introduction
Wind turbine wake modelling concentrates on characterizing the flow behind wind turbines. There are two main physical phenomena of interest in the wake: (1) the momentum (or velocity) deficit, which causes a reduction in the power output of the downstream turbines, and (2) the increased level of turbulence, which gives rise to unsteady loading on downstream turbines.
The wake-induced power losses and blade loadings are studied in two regions within the wake, referred to as near and far wake. The near wake starts right after the turbine and extends to approximately 2-4 rotor diameters (D) downstream [1,2]. In this region, the flow is highly influenced by the rotor geometry, which leads to the formation of the blade tip vortices. In addition, there are steep gradients of pressure and axial velocity, and wake expansion. In the far wake, the effects of the rotor geometry are limited to the reduced wind speeds and increased turbulence intensities. Further, the turbulence is the dominating physical property in the far wake [3]. In addition to the rotor induced turbulence, the region further downstream is influenced by the large scale (or atmospheric) turbulence. The turbulence mixing accelerates the wake recovery in terms of both the velocity deficit and the turbulence intensity. In the far wake, the velocity deficit approaches a Gaussian profile, which is axisymmetric and selfsimilar [4]. Moreover, the meandering of the wake might also contribute to the recovery of the velocity deficit although it significantly increases the unsteady loading on the downstream turbine(s). All these elements lead to different approaches for the development of wind turbine wake models. Out of the many, widely used six models and their interpretations developed at the Technical University of Denmark (DTU) are presented here. First, the components of wake modelling are described in order to demonstrate the differences between the modelling approaches better and then the benchmark study for onshore and offshore cases has been performed. The paper is organized as follows: in Section 2, the Navier-Stokes equations are presented with the incorporated turbulence modelling form of Reynolds Averaged Navier-Stokes (RANS) and Large Eddy Simulation (LES), discussed in Sections 2.1 and 2.2, respectively. In Section 3, the subcomponents of wake modelling, namely the inflow generation, the wake summation, the calculation of the wind speed at the rotor, the modelling of the wind turbine in the simulations, and the assessment of wind direction and speed are listed. Note that all of these concepts are originally much more comprehensive than their descriptions presented here and are only explained in the frame of wake modelling. The wake models Jensen, Larsen, Dynamic Wake Meandering, Fuga and the flow solver Ellipsys3D (both RANS and LES versions) are described in Section 4 and they are implemented on the onshore Sexbierum and offshore Lillgrund wind farms in Section 5. Accordingly, the models are evaluated in terms of their targets of application in Section 6.

Governing equations
It is convenient to say, except for the blade tip region, that the physics of wind turbine wakes can be described by the incompressible Navier-Stokes equations, where the atmospheric flow velocities upstream and downstream of a wind turbine typically range between 4 and 25 m/s. The governing equations in Einstein notation and Cartesian coordinates are: where u is the velocity and x is the position vector, P is the pressure, ρ is the fluid density, ν is the kinematic viscosity, f i ¼ F i ρ are the external body forces, t is the time, i, j are the directional components and S ij is the strain rate tensor defined as Since Eq. (2) includes a non-linear convective term, u j ∂u i ∂x j , especially in complex turbulent flows, some simplifications in both fluid and blade modelling are needed for computational purposes. There are therefore a large number of turbulence models, and from the computational fluid dynamics (CFD) point of view, we will concentrate in the RANS and LES methodologies.

RANS
To account for the turbulence effects, which can be of random chaotic nature, the instantaneous Navier-Stokes equations are time-weighted averaged resulting in the well-known RANS equations. The time-weighted average procedure is based on a statistical approach applied to the main variables of the flow and decomposes velocity into an average, u, and a fluctuation term, u 0 , the so-called Reynolds decomposition [5], When applied to the Navier-Stokes equations, the continuity equation becomes Note that u 0 j ¼ 0 and both u j ðx j ; tÞ and u 0 j ðx j ; tÞ are solenoidal because the flow is assumed to be incompressible which can be written as the left-hand side of Eq. (2) can be rewritten. Using the continuity relation as, The time averaging of Eq. (2), when considering Eq. (7), becomes The non-linear term u i u j À Á can be expanded using Reynolds decomposition.
Substituting Eq. (9) into Eq. (8) results in the final form of RANS equations.
where the term u 0 i u 0 j is the Reynolds stress tensor, which is a result of the non-linearity of the convective terms, and represents the averaged momentum transfer caused by turbulent fluctuations. The Reynolds stress tensor contains 6 new additional unknown variables that must be modeled (only 6 additional variables instead of 9 because of the imposed condition of the RANS angular momentum equation). A first approach to model the Reynolds stresses was first proposed by Boussinesq [6] who introduced the concept of turbulent viscosity or eddy viscosity. His idea was to describe turbulent effects as an increased fluid-fluid and fluidsolid viscosity interaction (surface forces). However, the approach assumes turbulence effects to be isotropic. The hypothesis basically relates the Reynolds stresses to the mean rate of deformation and can be simplified as, where k 1 2 u 0 k u 0 k is the turbulent kinetic energy, δ ij is the Kronecker delta, S ij is the mean strain rate tensor as defined in Eq. (3), and ν T is the turbulent eddy viscosity. The Boussinesq hypothesis states that the transfer of the energy mechanism between turbulent eddies is very much like that between molecular formations [7]. However, in contrast to ν, ν T is not a physical property of the fluid but it represents the turbulent characteristics of the flow. Furthermore, direct numerical simulations (DNS) have indicated that there is no correlation between the terms u 0 i u 0 j and S ij [8]. Therefore there is no physical basis for Eq. (11) and the assumptions are not valid for simple shear flows, anisotropic flows or 3D flows [9]. However, because they improve practicality and maintain robustness, the assumptions are applied to a variety of tools that provide solutions within certain accuracy. Several studies [10][11][12] offer detailed explanations of turbulence modelling and turbulence viscosity models.

LES
LES is a powerful technique to represent the turbulence characteristics of a flow by decomposing it into large and small scales. The small eddies are filtered out, so that the effect of large structures can be solved using the Navier-Stokes equations directly while small scale turbulent mixing is modelled. Eddies smaller than a certain grid size, Δx, are estimated using a subgrid-scale model. To eliminate the subgrid-scale, a filter with a width of Δx is introduced, which corresponds to the convolution of uðx; tÞ by the Therefore, the subgrid-scale field can be described by the difference between the actual and filtered flow, When filtering is applied to the combination of Eqs. (2) and (7) we get where T ij is defined as the subgrid stress tensor, is responsible from the momentum exchange between subgrid and filtered scales and is formulated as By using the Boussinesq hypothesis, the subgrid stress tensor can be rewritten in analogy to the RANS equations as; whereS ij is the filtered strain rate tensor and ν SGS is the subgridscale eddy viscosity which is most widely used as shown in Smagorinsky et al. [13], where jS j ¼ 2S ijS ij 1=2 , C s is the Smagorinsky constant, which varies between 0.1 and 0.2, depending on the properties of the flow [7], and Δ is defined as There are other subgrid-scale models, some of which can be found in the review of Lesieur et al. [14]. Determining Δ, which states the filter width, is crucial in order to represent the physical characteristics of the flow. An example of the accuracy studies performed on LES can be found in Geurts et al. [15].
The major concern regarding LES is the high computational cost for many engineering applications. However, particularly when dealing with flows with a solid wall as a boundary, it is possible to switch to RANS, since Eqs. (16) and (11) have similar characteristics. This is a hybrid approach called detached eddy simulations (DES). Specifically for high Reynolds number flows, where massive detachments may occur, DES has been shown to better represent the flow characteristics than both RANS and LES [16].

Subcomponents of wake modelling
After describing the equations for two of the main CFD approaches, RANS and LES, the "initial" and "boundary" conditions of those equations need to be introduced. Furthermore, concerning the analytical models, the approaches to account for multiple wakes and the way the rotor is characterized for wake simulations are described in this section.

Inflow generation
The characteristics of the atmospheric flow are mainly determined by the orography, the roughness and the roughness changes, together with the atmospheric stability which have been modelled using a variety of approaches. Here, for the description of the vertical velocity profile and turbulence, some fundamental concepts are discussed.

Logarithmic (or linear) law
The logarithmic wind profile can be derived in many ways [17][18][19] and it is formulated in the meteorological context as where U is the mean wind speed, u n is the friction velocity, κ is the von Kármán constant (E 0.4), z is the elevation above ground level and z 0 is the surface roughness length.

Power law
Another way to characterize the vertical wind profile is the power law, which is widely used in many wind engineering applications due to its practicality. It is given as where U ref is the undisturbed mean wind speed at a reference height, z ref , and α is the power law or shear exponent. In general α is a variable quantity ranging from less than 1/7 during daytime and more than 1/2 during the night [20].

Atmospheric stability
The static stability condition of the atmosphere has an effect on the flow characteristics and is normally taken into account by the Obukhov length, L S , which represents the ratio of the mechanical to convective turbulence production, where g is the gravitational acceleration, T is the air layers mean temperature and H; C P , and ρ are the kinematic heat flux, specific heat and density of the air, respectively. As u n Z 0 m=s, an unstable or stable behaviour of the atmosphere is determined by the sign of H. For example when dT dz 4 dT dz À Á adiabatic the atmosphere is considered to be unstable.
In many wind energy applications, the atmosphere is generally assumed neutral. However the atmospheric stability can have a large effect on the atmospheric flow behaviour for the inflow generation for wind turbines [21,22] and therefore it has become a growing research interest [23,24].
Atmospheric stability is taken into account in the log profile by including a correction term Ψ m , Ψ m is the integration of ϕ m which is the dimensionless wind shear ϕ m ¼ κz un ∂u ∂z . The form of ϕ m and Ψ can be found in [25]. Note that, Eq. (21) is most commonly used in wind energy without the term, Ψ m z 0 =L À Á , since turbines are deployed in areas with low z 0 . However, especially for complex terrain problems (forested areas, complex elevations, etc.) logarithmic law with stability correction should be used as described by Eq. (21).

Precursor turbulence box
A more sophisticated way to generate inflow conditions for wind turbine (or wind farm) simulations is by using a precursor turbulence box, in which a separate simulation without wind turbines under specific boundary conditions and assumptions is performed. In the studies of Bechmann and Sørensen [26], a precursor simulation was run over flat terrain with a set of parabolic equations (the Navier-Stokes equations with a boundary layer approximation), in which the pressure gradient is assumed to be constant, and the Coriolis forces are included.

Mann turbulence box
Mann [27] developed a spectral tensor turbulence model which can be used to simulate wind fields with particular turbulence characteristics [28]. Thus, it is now used for inflow turbulence generation for wake modelling, e.g. Dynamic Wake Meandering model and EllipSys3D, as discussed in Sections 4.4.3 and 4.6.2, respectively.
The Mann model looks at the spectral tensor of atmospheric turbulence at neutral stability state. The spectral tensor contains all information on spectra, cross-spectra and coherences that are required for engineering applications in wind energy. In the model, the Rapid Distortion Theory (RDT) [5] is combined with the "eddy lifetime" to describe the amount of shear, which gives the turbulence an anisotropic character. The model involves three adjustable parameters which can roughly be described as (1) a length scale that defines the size of the turbulent eddies, (2) a nondimensional parameter to estimate the eddy lifetime, and (3) a parameter related to the energy dissipation.

Wake summation
One of the subcomponents of wake modelling is the wake superposition concept. In order to include the effects of all the upstream turbines to the total velocity deficit, 4 approaches are mainly used [29]: where n is to the total number of upstream turbines, therefore u n þ 1 refers to the wind speed at the turbine in question. In addition to these approaches, Van Leuven [30] considers the effects of the closest upstream turbine in the WINDPARK model, which shows a good agreement with the measurements obtained in the Zeebrugge wind farm. Although this approach works fine for onshore conditions, for offshore wind farms where the wake effects are shown to be more dominant, Habenicht [31] underlines the importance of wake superposition methods. He has compared the superposition methods for four different offshore wind farms and showed that the linear and quadratic sums give the best results.

Rotor wind speed calculations
In this section, the methodologies to model the physical induction of the atmospheric inflow near the rotor disk are discussed.

Elliptic equations
The vortex system created downstream of the turbine induces a velocity component on the rotor axis in the direction opposite to the incoming atmospheric inflow. Characteristically in the elliptic problems, the disturbance signals, or a sudden change of information inside the domain, travel in all directions and affect the solution everywhere else. Hence, the change that occurred downstream of the turbine will naturally affect the modelled flow around the rotor, including the near upstream flow and the rotor itself. Therefore, elliptic solutions techniques for Navier-Stokes equations around wind turbines do not require to introduce any external induction, as it will appear in the flow automatically. However, the major drawback of the elliptic equations is their complexity, corresponding to higher computational cost.

Vortex equations
In the vortex modelling, the blades and the vortices are considered as lifting lines or surfaces. The vortex strength parameter is defined using the circulation which is highly related to the atmospheric inflow. This inflow is induced using the Biot-Savart law, and for a single vortex element of strength Γ the vorticity is given as where r is the perpendicular distance between the point p and the vortex filament ds, as illustrated in Fig. 1.

Parabolic equations
The parabolic form of the Navier-Stokes equations can be achieved by implementing the boundary-layer or thin shear layer approximation. The latter briefly states that the flow over a surface can be divided into two as the flow inside and outside of the boundary-layer region. Using such approximation, the Navier-Stokes equations can be simplified by omitting the diffusive momentum transport term through the principal direction of the flow. It essentially means neglecting the pressure gradient along the transverse direction, which is assumed to be much smaller than in the principal direction [32]. The solution procedure of the parabolized Navier-Stokes equations is relatively simpler so it is commonly used in many engineering applications [33]. The most common method for induced velocity component in parabolic flows is the actuator disk or 1-D momentum theory.
Actuator disk (1-D momentum theory) approach: This is based on linear momentum theory in which the wind turbine is modelled as an actuator disk, i.e. with an infinite number of blades. The flow before and after the actuator disk is considered to be steady, incompressible, homogeneous, isotropic, asymmetric with constant pressure profile, non-turbulent, inviscid, neutrally stable and non-rotational. Also, the thrust is assumed to be uniformly distributed over the disk area, and the velocity through the disk is considered to be constant. An illustration of the control volume used for the theory is shown in Fig. 2.
The velocity U 1 is induced using an axial induction factor, a, such that [34]; Also, the angular induction is introduced radially by the angular induction factor a 0 as where U 2rot is the induced tangential velocity at the rotor plane, Ω is the rotational velocity of the rotor, and r is the radial distance from the rotational axis. An example of an iterative calculation procedure for a and a 0 can be found in Manwell et al. [35].

Wind turbine model
In this section, some of the methodologies followed throughout the literature to estimate forces applied to the rotor plane are presented.

Inverse 1-D momentum theory approach
The axial force term that appears in Eq. (2) can be modelled using the inverse momentum theory approach, where a general aerodynamic expression for a uniformly loaded actuator disc can be written as where f x is the axial force, c T is the thrust coefficient and A is the rotor swept area. In Eq. (26), the definition of the reference inflow wind speed, U ref , is not always straightforward for a wake-affected downstream turbine. For the upstream turbines U ref ¼ U 1 , whereas for the downstream ones, following the studies of Prospathopoulos et al. [36], U ref is defined using the local velocity field and an induction factor as in the 1-D momentum theory so that U ref ¼ U local =ð1 ÀaÞ.
They proposed an iterative approach by assigning an initial value to U ref and determining the thrust coefficient. The axial induction factor can be approximated accordingly, using the thrust curve of the rotor with c T ¼ 4að1 À aÞ. The process continues until convergence is achieved for U ref . Finally, that value is used to estimate the axial force applied on the rotor using Eq. (26). Note that U ref is often over-predicted (2-3%) when estimated using the inverse 1-D momentum theory, leading to 5% and 8% error in total thrust and power, respectively [37].

Induced thrust curve approach
In this approach, the thrust coefficient is defined as a function of the relative or induced velocity, U ref instead of the free stream wind speed ðU 1 Þ where the conventional thrust curves are calculated accordingly i.e. c T vs. U 1 . In order to create a newly defined thrust curve, a CFD algorithm was created for a relatively simple individual turbine case and run at different wind speeds from which the relative velocity and axial force values can be extracted. As a result, a new thrust curve in terms of induced velocities for that specific wind turbine is constructed and, therefore, can be used in more comprehensive calculations [37].

Blade element momentum (BEM) theory and generalized actuator disc model
BEM theory is one of the first and still the most commonly used methodologies to investigate rotor aerodynamics and it is described in many studies in the literature, e.g. see Hansen [38]. In this section, the application of BEM theory to estimate the forces across the rotor are explained.
The forces on the rotor are calculated using the geometrical components of the aerodynamic sectional lift, L, and drag, D, forces, which strongly depend on the aerodynamic characteristics of the airfoils, where c L and c D are the lift and drag coefficients, respectively, e L and e D are the unit vectors, c is the chord length of the airfoil, B is the number of blades of the turbine and V rel is the relative velocity. The rotational effects are taken into account in the estimation of the V rel as where W z is the induced velocity such that a ¼ W z =U r ef , r l is the local radius of the considered annular section, and W θ is the induced angular velocity. The axial and tangential forces are determined in terms of the flow angle, ϕ l defined between the direction of the V rel and the rotor plane, and the sectional forces become In the BEM method, Eqs. (27)- (30) are numerically solved by iterative algorithms (e.g., the Newton-Raphson algorithm) by estimating the axial induction factor using a similar approach as described in Section 3.4.1, whereas in the generalized actuator disc model, the components V z ¼ U 1 À W z and V θ ¼ ÀW θ are measured on the disc [39].

Sequentially activation method
In this approach, a CFD simulation is performed for only the upstream turbines in the wind farm, eliminating the induction effect of the downstream turbines. As a result, the incoming velocity, U 1 , at the location of the downstream turbine can be defined and a corresponding c T can be determined using the conventional thrust curve provided by the manufacturer of the turbine. The procedure is repeated until the whole wind farm is computed and the axial forces are determined accordingly.

Aero-elastic model approach
The axial force or the thrust is known to be the dominant force in relation to bending moments on the wind turbine, which can be measured using e.g. strain gauges. These measured bending moments are used inversely, together with aero-elastic models, e.g. HAWC2 [40], to estimate the thrust, which is spatially integrated over the rotor. Since the thrust is not the only force causing bending on the turbine structures, but those resulting from the interaction between the turbine and the complex atmospheric flow, aero-elastic models are used to estimate the moments that are dominated by the existence of the axial forcing [41].

Wind direction & speed
Accurate wind direction data is key in wake modelling since the direction defines the path of the wake. Therefore it determines the full and partial wake conditions at the downstream wind turbine positions, which are critical in both wind turbine loading and power production calculations.
Typically in wind-power meteorology, the wind direction is determined either using the measurements from a meteorological mast or the yaw angle extracted from the Supervisory Control And Data Acquisition (SCADA) system of the turbine(s). However, it is well known that the conventional assumptions introduce a considerable amount of uncertainty mainly caused by the physical distance between the meteorological mast and the farm, the sensitivity of yaw measurement techniques also known as "yaw misalignment".
For example, in the interface of Fuga [42], the model is further described in Section 4.5, the turbine site locations and the wind data are input as done in Wind Atlas Analysis and Application Program (WAsP) [43]. The wind direction measurements can be post-processed by either simple averaging, Gaussian averaging [44], or considering the meandering of the wake in which a spatial correlation is activated to account for the direction uncertainty. In Fuga the meandering is taken into account by creating a curve by joining 10-min averaged wind direction values and considering the probability of the difference between this curve and the instantaneous values.

Wake models
In this section, an extensive conceptual review of the wake models developed at the DTU is presented.

Infinite wind farm boundary layer model
An infinite wind farm boundary layer (IWFBL) model was developed by Frandsen [1]. In the model, around the turbine rotors, i.e. the "rotor layer", the velocity profile is reduced compared to that above hub height and both profiles are logarithmic as shown in Fig. 3.
Inside the wind farm the turbines are assumed to be evenly spaced at a distance x and a dimensionless separation between the turbines is defined as; s ¼x/R, the term R being the radius of the turbine. As shown in Fig. 3, the difference in shear stresses around the turbine hub height is where u h is the asymptotic spatial average wind speed at hub height, and t is the simplified thrust term t ¼ ÀρC 0 Therefore the relation between the friction velocities and u h can be calculated using where u n2 is the friction velocity above and u n1 below hub height. Note that, under the rotor layer, the logarithmic wind profile is valid and can be used to relate u h to u n1 using logarithmic law. For the region above that layer, the simplified geostrophic law [45] is applied and the resulting expression is found as; where h is the hub height, f p ¼ f c expðA n Þ with f c being the Coriolis parameter, A n a modified A parameter from the resistance-law constants, and G the geostrophic wind speed. The friction velocities are parametrized as; Substituting them into Eq. (32) and solving for u h yield After solving for u h , the friction velocities, u n1 and u n2 , can also be calculated.
Additionally, Frandsen [1] approximated the wind speed reduction at hub height, R u ¼ u h =u 0 , where u 0 is the undisturbed wind speed at the same height, as where γ ¼ 0:025=ln h=z 0 À Á 1=3 , with z 0 being the surface roughness.

Atmospheric stability correction
Peña and Rathmann [46] added atmospheric stability effects to the IWFBL model extending the logarithmic wind profile to account for atmospheric stability using a correction term depending on the dimensionless wind shear. The wind speed reduction has a similar form as that of Frandsen [1], but both K 1 and K 2 are modified to include atmospheric stability by adding/subtracting the stability function, Ψ m ðz=L s Þ, where A n ¼ ln G fz 0 À κG un and μ 0 ¼ κu n =f c L s , where u n2 for the section above rotor layer is formulated as z 00 is the effective roughness length of the wind farm,

The Jensen wake model
The Jensen wake model is one of the most popular models among engineering applications due to its simplicity, practicality and robustness. The description is based on the studies of Jensen [47] and Katic et al. [48].
Using the control volume presented in Fig. 4, where D ¼ D r is the rotor diameter, and assuming a top-hat inflow profile the mass balance between the rotor plane and the downstream flow yields, Also, the wake is assumed to be expanded linearly as a function of the downstream distance x at a rate α, D w ¼ D r þ 2αx and u r =u 0 ¼ 1 À 2a using the axial induction factor, the fractional decrease in wind speed, a ¼ u 0 À ur u 0 . Putting them into (42), the normalized velocity can be found as Assuming ideal axially symmetric flow, no rotation, no turbulence and conic shape wake profile, the axial induction factor can  also be written as  [43] is based on the Jensen wake model and accounts for the effect of multiple wakes on the velocity. In the original version of Katic et al. [48], the ground interaction of the wake is taken into account by assuming an "underground rotor", which is a reflection of the original one. To derive the efficiency of a wind farm, the combination of the effects of four different overlapping mechanisms is considered: 1. Directly upwind rotor wakes. 2. Reflected upwind "underground rotors". 3. Shading upwind rotors, located left or right of the directly upwind rotor. 4. Reflected shading upwind rotors, located left or right of the wind direction.
The local wakes are superposed to estimate the velocity deficit at the nth turbine δ n ¼ Infinite row of turbines: Jensen [47] already estimated a model for the velocity deficit of an infinite row of turbines based on his wake model. If the velocity at the last partition of the infinite row of turbines is defined as u inf then Infinite Park Wake model: Considering the effects of four overlapping of wakes in the Park wake model the total wake deficit δ T is estimated as the quadratic sum of four types of wakes [22], Rathmann et al. [49] have solved those effects analytically and Peña and Rathmann proposed [46], where δ 0 is the initial wake deficit, δ 0 ¼ 1 À ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À C T p , s r is the dimensionless stream-wise separation between turbines, i.e. s r ¼ x=D and s f ¼ y=D with y being the cross-wind turbine-turbine distance. Wake decay coefficient: When using the Park model in WAsP, the wake decay coefficient term α is by default α¼ 0.075. In the study of Peña and Rathmann [46], the wake decay coefficient was shown to be a function of height roughness, atmospheric stability and turbulence separation. For practical purposes, the below expression is recommended.
which showed very good agreement with data from the Sexbierum [50] and the Horns Rev-I wind farms [51]. Using the similarity theory, α can be related to the turbulence intensity, TI, as α % 0:4TI.

1988 (Early) version
Larsen [52] has introduced a simple wake calculation procedure which was implemented in the commercial software WindPRO [53]. In the model, the axis-symmetric form of the RANS equations with the thin shear layer approximation is used. The pressure term appearing in the parabolic equations was also neglected and the turbulence closure, ν T , was represented using Prandtl's mixinglength theory as where l is the mixing length and S ij is the strain rate tensor. The problem is assumed to be steady, axisymmetric and self-similar along the perpendicular direction to the flow. Larsen considered the solution of the RANS equations using first and second order approximations. In the first order approximation, the expression to be solved together with continuity equation is simplified as: u x is the wake perturbation of the inflow along the axial direction and r is the radial direction and x is the axis of symmetry.
In order to solve Equation (51), two boundary conditions are defined: (1) u x ¼0 on the boundary of the wake, and (2) U 1 ⪢u x , which is obtained by writing the momentum balance assuming the inflow velocity to be much higher than the axial wake perturbations. Using those conditions, the radius of the wake, r w , and the axial (u x ) and radial (u r ) wake perturbations are found as: where C T is the thrust coefficient, A is the rotor swept area, and c 1 is a constant that is defined empirically [52]. The second order system uses the full form of the RANS equations, which were later found to be negligible for most engineering applications [52].

2009 (Later) version
The main improvements in the 2009 version of the Larsen model [54] compared to the 1988 one are the boundary condition (s) and the wind farm approach because the early version was derived considering the single wake case only and provided no solution for multiple wake situations. The later version of the model defines the boundary conditions using the results of the analysis of full scale experiments. The first boundary condition is defined at the rotor plane and the second one is defined at a fixed frame of reference placed at a distance 9.6D downstream.
The second order approximation is neglected in the later version as well, and the wake radius and velocity deficit resulting from the updated boundary conditions are where x 0 ¼ 9:6D with R 9:6D being the wake radius at 9.6D, which is empirically calculated using the analysis performed for the Vindeby offshore wind farm, and expressed using atmospheric turbulence intensity, I a as, where the constants a 1 ; a 2 ; a 3 ; a 4 and b 1 are defined in Larsen [54]. The wind farm approach is considered using two different methodologies to calculate the inflow speed: the geometric (or linear) averaging and momentum balance, which are respectively, where U is the incoming ambient velocity modelled by logarithmic wind profile. The velocity inside the wind farm is calculated using the linear averaging as: where M is the number of upstream rotors that generate wakes affecting the rotor m.
For the non-linear approach, the decomposition of U m cannot be performed linearly, thus, the velocity profile imposed on rotor m may be described as; Eqs. (63) and (64) are solved using a 4-point Gauss integration method, which is explained in detail in Larsen [54]. Additionally, in both of those equations it can be seen that the multiple wake effects are superposed using the linear sum.

Dynamic wake meandering model
The dynamic wake meandering (DWM) model describes the wake as a passive tracer driven by the large-scale turbulence structures in the atmospheric boundary layer. The model may be further investigated using the studies performed by Larsen et al. [55,56] and Madsen et al. [57]. The recent improvements to the model and the validation cases are presented in [58][59][60][61].
The DWM model consists of three elements, which together describe the essential flow characteristics behind a turbine: (1) Velocity or wake deficit; (2) Meandering of the wake; and (3) Rotor added turbulence. Here they will be considered separately.

Velocity deficit
In the DWM model, the velocity deficit is initialized by the pressure gradient and formulated in the meandering frame of reference. The profile behind the turbine is assumed to be axisymmetric and steady. Parabolic Navier-Stokes equations with neglected pressure terms are used and the resulting equations are where V r denote the mean velocity along the radial direction. The eddy viscosity ν T is mainly described by the methodology proposed by Ainslie [62] and manipulated to include ambient turbulence intensity, where k 2 is an empirical constant for the flow field, b is the instantaneous wake half width, U def ;min is the minimum wake wind speed, U H is the wind speed at hub height, I amb is the ambient turbulence intensity at hub height, k amb is a calibration constant, and F 1 and F 2 are filter functions depending on the downstream distance x only.

Meandering of the wake
As mentioned earlier, the DWM model assumes the wake to behave as a passive tracer transported in a large-scale turbulence field, where eddies larger than two rotor diameters. Therefore, the large-scale transversal and vertical velocities, v and w respectively, are important.
The displacement of the wake is defined by the characteristic velocities, where ðx b ; y b ; z b Þ are the inertial coordinate system fixed to the turbulence box introduced, and A f is the averaging area most logically selected as a circle in which the origin is assigned as ðy b ; z b Þ with a diameter D w . The transversal and vertical wake displacements are described as, where T is the time interval considered in the "snapshot" associated with the Pseudo-Lagrangian approach formulated by T ¼L/U with L being the along-width length of the turbulence box considered, and t i is the time when the velocity deficit is released. Additionally, the initial conditions at the time t i are given as, The solution to Eqs. (68) and (69) together with Eqs. (70)-(73) are presented in detail in Larsen [56] including methodologies for their simplification and a numerical algorithm. Finally, it should be noted that the wake deficit is not affected by the meandering progress.

Rotor induced turbulence
The rotor induced turbulence, or wake added turbulence, in the meandering frame of reference corresponds to the small scale turbulence, namely the tip, root and blade bound vortices, as well as the wake shear layer. In the DWM model, the wake added turbulence at a particular downstream position is modelled using an isotropic Mann turbulence box (Section 3.1.5), with cross sections corresponding to one rotor diameter. Additionally, the added wake turbulence intensity is assumed to be rotationally symmetric and does not influence the meandering or velocity deficit processes.
In summary, the resulting turbulence in the DWM model includes components from meandering, added wake turbulence and ambient turbulence. The resulting velocity field may be expressed as where U res , v res and w res are the axial, lateral and vertical velocity components of the resulting velocity field, respectively. The subscript aw represents the added wake component, and the subscript amb is the ambient contribution. Finally, U m denotes the unsteady velocity component obtained from the meandering of the velocity deficit, which is determined as with x m and r m being the downstream axial and radial coordinates, respectively, updated at each time step.

FUGA
Fuga is a fast engineering tool based on the linearized RANS equations. It uses a system of look-up tables to construct the velocity field behind a turbine, and it uses linear summation to consider multiple wake cases. Due to its simplicity in wake modelling, Fuga is one of the most robust computational fluid dynamics (CFD) based models established for wake effects' calculations. The methodology presented in this study is based on the works by Ott et al. [42].
The Cartesian form of the RANS equations are used with a simple closure, where the eddy viscosity is equal to that usually used within the atmospheric surface layer.
Since the equations are not parabolized, there is no need to artificially induce the rotor velocity where the atmospheric inflow is modelled using the logarithmic wind profile including the stability effects. The drag forcing term is modelled using an actuator disk model with a layered control volume as, where δ is the Dirac delta function and Θ is a step function, which is equal to zero for negative and 1 for positive arguments. Due to the fluctuations related to the existence of a step function, the drag calculations are smeared out.
The simplified RANS equations are linearized using Taylor expansion and only the terms with order zero and one are considered. The zeroth order equations correspond to the case without any perturbations to the flow, meaning that there are no turbines. The drag force of order one, f x 1 , is defined by fitting the first order equations to a Chapeau function. The resulting equations are further simplified using Fourier transformation in which two mixed spectral variables are defined along x and y directions. A new numerical scheme is implemented to overcome the difficulties of solving a linearized model for flows over small values of z 0 which is the case for offshore sites with low roughness lengths and where the wakes are more pronounced. The scheme is described in detail in Ott et al. [42] together with the validation of the model for certain test cases.

EllipSys3D
Ellipsys3D [63,64] is a 3D general purpose CFD solver with a block-structured finite volume approach. Both RANS and LES models are available in EllipSys3D and can be further examined in Sanderse et al. [65].

RANS
In the RANS version of Ellipsys3D, the rotor is modelled as an actuator disk, the elliptic form of the Navier-Stokes equations are used thus no external induction is introduced, and the non-linear terms, u j ∂u i ∂x j , are discretized using the QUICK scheme [66]. Ellipsys3D can use a number of turbulence models. One of the latest developments is the k-ε-f P model [67], which is a modified version of the widely used k-ε model from Launder and Spalding [68]. Where the standard k-ε model fails to predict the velocity deficit in the near wind turbine wake [67,[69][70][71][72], the k-ε-f P model has shown good agreement with LES and measurements for single [67], double wake cases [73], and complete wind farms [74]. In the turbulence models, the turbulent eddy-viscosity is defined as: where C μ is a constant, k the turbulent kinetic energy, and ε the turbulent dissipation. In the standard k-ε, f P ¼1 and the effective eddy-viscosity coefficient C μ f P is a constant. In the k-ε-f P model, f P is a scalar function that depends on the local shear parameter . The scalar function f P in the k-ε-f P model is defined as whereσ is the shear parameter in an idealized (logarithmic) neutral atmospheric surface layer and C R is a calibration parameter. In the neutral stability solution, f P ¼1 because σ ¼σ . In regions with a high shear parameter, i.e. σ 4σ , f P o 1 and the turbulent eddy viscosity from Eq. (80) is decreased. The near wind turbine wake is characterized by high velocity gradients, where σ⪢σ . As a result, the k-ε-f P eddy viscosity model delays the wake recovery compared to the standard k-ε. It should be noted that C R controls the magnitude of the delayed wake recovery. The constant C R is calibrated against LES for eight different single wind turbine cases [67]. The same transport equations for k and ε are used in both turbulence models, where P u 0 i u 0 j U i;j is the turbulent production, ν is the kinematic molecular viscosity and C ε;1 ; C ε;2 , σ k , σ ε are constants. The values of the constants are listed in Table 1.
When the standard k-ε is applied to atmospheric flows, it is common to control the ambient turbulence intensity at a reference height I H;1 with C μ using Subsequently, one of the model constants from Table 1 is adjusted to maintain the logarithmic solution [75] ffiffiffiffiffi ffi However, the behaviour of the f P function changes when C μ is altered, which is not desired. Therefore, the ambient turbulence intensity is set with z 0 , and the friction velocity u n is adapted to set the free-stream velocity. As a result, the velocity inflow profile differs from the measured profile, although the difference in the rotor area is only in the order of a few percent.

LES
In the EllipSys3D LES, a combination of Eqs. (14)- (16) is performed. The non-linear terms are discretized using a hybrid scheme formed of QUICK and fourth order central differencing schemes and the empirical constants are chosen based on the studies of Troldborg related with the actuator line [76]. There, the rotor is modelled as an actuator line, and the axial force is defined using BEM, which requires sectional aerodynamic characteristics of the blades but increases the efficiency of the calculations. The inflow velocity profile is defined using the logarithmic wind profile (see Section 3.1.1), and the Mann Turbulence box (see Section 3.1.5) is used as inflow turbulence, where neutral atmospheric conditions are assumed.

Introduction
Sexbierum is an onshore wind farm located in the Northern part of the Netherlands at approximately 4 km from the shore on homogeneous flat terrain, mainly grassland. It consists of 18 turbines with a total installed capacity of 5.4 MW. The layout of the wind farm is presented in Fig. 5.
The turbines in the farm are HOLEC WPS 30-3 [77] with a rated power of 310 kW, a rotor diameter of 30.1 m and 35 m hub height. These turbines are pitch regulated with a cut-in wind speed of 5 m/ s, a rated wind speed of 14 m/s, and a cut-out wind speed of 20 m/s. For Sexbierum case, two benchmarks were defined; the single and double wake cases by Cleijne [78,79].
B1 -single wake: In the single wake test case, the comparison between the simulations and the measurements is performed in the wake of the turbine T18. The met masts are placed 2.5, 5.5 and 8 diameters downstream of T18, and the wind speed measurements during 6 months provided in Cleijne [78] are considered in the benchmark. The observed wind speeds are between 5 and 10 m/s, where for the simulations 8 71 m/s incoming wind speed is considered as this is the most frequent wind speed bin observed. The roughness length and the turbulence intensity at hub height are estimated to be 0.049 m and 9.5%, respectively. Neutral atmospheric stability is assumed for the wake computations but discussed afterwards.
In this benchmark, the results of the Jensen model, Larsen model, Fuga, and Ellipsys3D RANS and LES solvers are presented.
Two RANS turbulence models are tested: the standard k-ε model and the k-ε-f P model. The RANS computations are performed with a domain size of 25D Â 16D Â 8D. The inlet is defined at 5D upstream of T18 and the refinement of the mesh in the wake region is performed in such a way that there are 10 cells per rotor diameter as proposed by van der Laan et al. [67] to obtain good resolution in the near wake region with Ellipsys3D RANS. The high resolution area starts at 3D downstream from the inlet with a width and length of 4D and 14D, respectively. Vertically, it starts from 0.5D below hub height and goes up to 1.5D above the rotor. The mesh has a maximum expansion ratio of 1.2 with an initial height of z 0 on the ground. The computational domain for Ellip-sys3D RANS includes 1.57 million cells, and the boundary conditions for that domain are: (1) rough wall condition at the ground surface, (2) symmetric boundary conditions on the sides, (3) inlet velocity condition at the inlet and top, and (4) far field outlet boundary conditions.
The time required to run Ellipsys3D RANS for the single wake case with a convergence criterion of 10 À 5 , i.e. the iteration is terminated when the difference between two calculation steps falls below 10 À 5 , is 3-min and 3-s with a time step of 0.008-s and 48 CPUs.
The computational domain used for Ellipsys3D LES is the same for both the single and double wake cases and its dimensions are 15D Â 15D Â 23:25D. The inlet boundary is located at 7.35D from the first upstream turbine (T18 for single wake, T38 for double wake case). The grid points are distributed uniformly in such a way that there are 30 points corresponding to each rotor, and two refined regions to resolve the inflow and wake turbulence. The first high resolution region is located at 0.35D upstream of the first turbine and extends to 10.3D downstream with a height and width of 1.8D. The second highly resolved area is located at 1.8D upstream of the first turbine where the inflow turbulence is introduced. It has the same height and width of 1.8D and it extends to 1.9D. The computational domain has 19.7 million grid points and the boundary conditions for that domain are: (1) no slip condition at the ground surface, (2) periodic boundary conditions on the sides, (3) far-field velocity on the top, (4) inlet velocity and turbulence as described in Section 4.6.1, and (5) unsteady convective outlet conditions. The computational time required to run the single wake case using the Ellipsys3D LES for 12-min real-time with a time step of 0.008-s and 150 CPUs is approximately 4 days and 4 h. However, it should be noted that the computational domain used for the single wake case was actually optimized for the double wake case. Therefore it should be expected that the performance of the simulation in terms of the computational costs can be enhanced by simplifying the mesh according to the single wake requirements.
B2 -double wake: In this benchmark, the power measurements of turbine T36 in the wake of T38 and T37 (see Fig. 5), covering a period of 3 months are studied. Similar to the single wake case, the wind speed interval is 5-10 m/s for the dataset. The model simulations are performed at 8 m/s and the roughness length and turbulence intensity are 0.045 m and 9.5%, respectively, as recommended by Cleijne [79].
For the double wake case, the results of the Jensen model, Larsen model, Fuga, and Ellipsys3D RANS (using the same two turbulence models) and LES solvers are presented. Four relative different wind directions [0°, þ 7°, þ 14°, þ21°] are simulated in LES, and 22 min of data with a time step of 0.008-s and 150 CPUs took over a week to run, which corresponds to approximately 27 500 CPU hours.

Results and discussion
B1 -single wake: In the Sexbierum single wake case two different wind direction averaging techniques are applied to both the Jensen (with wake decay coefficient, α¼ 0.04) and the Larsen results; a wind direction sectoral averaging (BinAve) and a Gaussian averaging considering the wind direction uncertainty [44] (GauAve). For these two models, the bin averaging is performed for a 30°wind direction span where 2.5°simulations are run and averaged over 5°bins. The Gaussian averaging method was applied with a standard deviation of 5°in the wind direction, although no information about the wind direction uncertainty is provided in the corresponding report [78]. Note that, the free stream wind direction is measured via the meteorological mast and given relative to the line connecting turbines T18 and T27, denoted as 0°.
The wake model performance in the near wake region is compared in Fig. 6. All the model and solver results, independent of the post-processing method, considerably deviate from the measurements. The data were collected approximately for only 43 h, therefore they are not statistically representative.
The deviation might be due to the atmospheric conditions. However, the atmosphere is very likely to be stable during the measured period [50]. With a low turbulent mixing of the stable atmosphere, the wake takes longer to recover, which explains the depth of the measured wake together with the under-estimation of the models, which are valid for neutral conditions, especially in the near wake regionsee Fig. 6(a).
Looking at the performance of the models, Fuga, the Larsen and the standard k-ε RANS underestimate the velocity deficit, and the k-ε-f P compares well with the LES (note that the LES results are not Gaussian averaged though). The Jensen model, the simplest of all, with Gaussian averaging provides very similar results to those of LES at 2.5D and 5.5D downstream distances. The Jensen model outperforms the others because of the wake decay coefficient used in the simulations. For onshore sites the recommended value is α¼0.075, whereas in our case α¼ 0.04 but as shown in Peña et al. [50] α could be even lower.
B2 -double wake: Fig. 7 is a combination of the results of the Jensen (with a wake decay coefficient of α¼0.04), the Larsen and Fuga models, EllipSys3D RANS k-ε-f P and Ellipsys3D LES solvers where the runs were performed in a similar manner as in the single wake case using a 5°bin and a wind direction step size of 2.5°.
Similar to the single wake case, the models generally underpredict the wake losses which might be due to the stable atmospheric conditions (9.5% turbulence intensity for a 0.045 m roughness length).
Apart from the performance of the Jensen model, especially the bin averaged version, the Larsen model with bin average seems to perform very well compared to the LES results. This might be due to the turbine spacing for the double wake case ð10DÞ which is very close to the distance that the model is calibrated ð9:6DÞ. On the other hand, the Gaussian averaged version of the Larsen, Fuga and RANS k-ε-f P turbulence model results seem to deviate from the measurements. Note that the Gaussian averaging takes into account the wind directions that might highly differ from the mean value during the averaging time period. However, for less turbulent cases, the deviation of the wind direction from the mean value is small. In that case, the wake deficit profile is over-smeared and the details are lost.
Notes and remarks about the Sexbierum wind farm case: The benchmarks included in the Sexbierum case were constructed using the reports of Cleijne [78,79] and all the models except for DWM were used for the simulations. In general, the models seem to deviate from the measurements significantly, which may be a consequence of the probable stable characteristic of the atmosphere. The Jensen model with a low wake decay seems to provide very good results although it is the simplest model used. It is concluded that the post-processing approach for stable cases should be revised and differ from the one developed for the neutral atmospheric conditions.
Note that in all the cases, the measured wake deficit profile is far from being symmetric which may occur due to the onshore effects such as terrain complexity, etc. However, it is not easy to tell since the dataset covers only a short period of time and the observations might be biased in terms of the atmospheric stability by the seasonal variation of the atmospheric stability. The wind farm wake modelling, especially for onshore sites, requires more inputs to model the inflow and the Sexbierum case is another example of data issues encountered in wind farm simulations.

Introduction
The Lillgrund wind farm is located in Øresund, 6-8 km from the Swedish west coast and south of Malmø. It consists of 48 SWT-2.3-93 wind turbines with a total rated capacity of 110 MW. The Lillgrund wind farm has an irregular layout with a gap in between, and the internal spacing of the turbines is 3.3 and 4.3D rotor diameters, as shown in Fig. 8. In the EERA-DTOC report for the Lillgrund wind farm test case [80], four benchmarks were specified as listed below.
The benchmark consists of 4 main cases: B1 -sector variation: The power deficit along complete rows with internal spacing of 3.3D and 4.3D is simulated to test the sensitivity of the models to the flow direction. The roughness length is 0.0001 m, the inflow mean velocity at hub height is 9 m/s and the inflow turbulence intensity at hub height is 6%, which is estimated based on sector wise-long term measurements of the met mast.
Two different rows that do not have a missing turbine are used and for the 3.3D case the wind direction is in the interval 120 715°, whereas for the 4.3D case it is 222 715°. Note that, the runs are performed at every 2.5°step for both arrays.
B2 -speed recovery: The power deficit along a row with missing turbine(s) and internal spacing of 3.3D and 4.3D are observed. In addition, the sensitivity of the models to the flow direction together with the speed recovery due to the missing turbines is tested. The input data and the characteristics of the runs to be performed are the same as in the previous benchmark, B1.
B3 -power deficit as a function of turbulence intensity: The calculations are performed for different inflow turbulence intensity levels at hub height (2-12%) with the same inflow conditions as in the previous benchmarks. Two different runs are performed for both 3.3D and 4.3D spacings using only the first two turbines in the row and the wind direction sectors are 1207 2.5°and 222 72.5°, respectively.
B4 -park efficiency: The wind farm park efficiency is defined as the ratio between the wind farm total output power and the power of the wind farm assuming undisturbed inflow for each turbine. Similar input data as in the previous benchmarks is considered and the inflow sector is taken as 0-360°with a span of 3°.

Results and discussion
We use a wake decay coefficient of 0.04 for the Jensen model with a quadratic sum for the wake summation, whereas for the Larsen model a linear summation is applied. Additionally, the thrust coefficients in both models are those provided by the turbine manufacturer.
B1 -sector variation: In Fig. 9, two different wind direction averaging techniques with 3 different models are run for this case with 3.3D and 4.3D spacings. The simulations are run at 2.5°step wind directions and averaged over 5°bins. The Gaussian averaging is applied for a 5°standard deviation in wind direction. Additionally, the same technique is applied to Fuga where the uncertainties in wind direction are taken into account using a Gaussian distribution of 4.9°. The Larsen model and Fuga under predict the wake losses for the second turbine placed at 3.3D and 4.3D. Both models are however designed to simulate the flow behaviour at much larger downstream distances. On the other hand, as shown in Fig. 10(a), the wake deficit under-prediction is compensated with a good prediction for the following rows, especially for the Larsen model.
The models are shown to perform better for wider wind direction sectors in Lillgrund by Gaumond et al. [82]. In our case, the EllipSys3D RANS k-ε-f P model over-performs to estimate the power deficit at the second wind turbine, because the f P function delays the wake recovery compared to the standard k-ε model.
B2 -speed recovery: In Fig. 11, the recovery point is clearly seen at 16.5D for 1207 2.5°and 17.2D for 222 72.5°. All the models capture the recovery and for this particular case the Larsen and the k-ε-f P model seem to estimate the power production reasonably well, especially after the second turbine. Both the Jensen and the Larsen models produce better results with the post processing of the wind direction uncertainty using a Gaussian distribution, which was also the case in previous benchmark, B1. Fuga seems to over-predict the power production for the first downstream turbine and then under-predicts the power production for the following turbines including the recovery point in the 3.3D spacing case. However, for the 4.3D spacing case, and similar to the previous benchmark, agreement between Fuga and the measurements is improved.
B3 -power deficit as a function of turbulence intensity: The standard uncertainty of the power deficit for the turbulence case is represented by the error bars in Fig. 12 with a confidence level of 68% for the SCADA results [83].
Since the original Jensen model does not consider the variations in turbulence, it remains constant. Both the Larsen model and Fuga significantly deviate from the measurements, especially for high turbulence levels. Increasing turbulence intensity levels show larger lateral wind components. That results in greater wind direction variations. Due to those large variations and narrow sectors, both models might fail to reproduce the observations well. The necessity to use the local turbulence intensity in simulations is addressed later in this study. Furthermore, a pragmatic approach to introduce dynamic effects to the engineering wake model is developed and presented in the following chapter.
B4 -park efficiency: The error bars indicated in Fig. 13 correspond to the uncertainty of power deficit with a 68% confidence level for the SCADA results [83]. The improvement of the model results by post-processing the wind direction uncertainty using a Gaussian distribution is considerable. Those Gaussian averaged versions of the Jensen and the Larsen model show a fair agreement with measurements. However, significant differences around   Remarks about Lillgrund wind farm case: Especially for the first three benchmark cases, the narrow wind sector of 5°(72.5°) is the major source of uncertainty since the data is 10-min averaged and most probably includes wind directions outside of that range.
Overall, the Larsen model and Fuga performed in a similar manner when considering a Gaussian distribution for the direction uncertainty, in agreement with the results obtained for Horns Rev Wind Farm [44]. In such a layout with small turbine spacings, the k-ε-f P closure of the Ellipsys3D RANS is seen to capture well the wind speed at the closest turbines through downstream. It can be said that even though the direction bins are narrow ð 7 2:51Þ for 10-min averaged data, the performances of all the models were considerably good in all benchmarks in general.

Application of the models
In this section, the application of the wake models developed in DTU will be discussed in terms of their typical usage, validity, accuracy, complexity, the uncertainty of the required inputs and computational costs.

Typical usage
The WAsP version of the Park model based on the Jensen model is targeted for wind farm planning and annual energy production (AEP) estimates. Due to its simplicity and practicality, it is often used to perform preliminary studies which are then improved with more sophisticated models.
Similarly, the Larsen model, also implemented in WindPro, is used for both single wind turbine and wind farm design and development stages. Fuga is a relatively new model. However, its robustness, speed, and promising results have already made it popular in the wind energy industry, and it is recently implemented in WAsP. The results showed that especially for the Lillgrund offshore wind farm case, Fuga and the Larsen model provide good results and are comparable to those of the more sophisticated models in offshore.
The DWM model is not only developed to be able to estimate the power production losses due to wake effects but it can also  calculate the loading caused by the wake effects. It is implemented in the aeroelastic code HAWC2 and calibrated accordingly but unfortunately it was not available for the present work.
The LES version of Ellipsys3D, due to its complexity and computational cost, is run by a limited number of users and mainly for academic purposes. Additionally, both RANS and LES simulations are conventionally performed for small number of turbines rather than large scale wind farms for the same practicality reasons.

Accuracy
Most of the models proposed show a fair agreement with the observations especially when they are post-processed to take into account the wind direction uncertainty or the atmospheric stability conditions. Physically, the models with more realistic considerations, or in other words less simplifying assumptions are more successful in simulating wake characteristics in detail. Thus, in general, the more complicated models are more likely to be more accurate. However, the Sexbierum test case in Section 5.1 showed that even the most sophisticated models can fail to reproduce the flow characteristics when the inputs are erroneous or deficient.

Complexity and uncertainty of inputs
The quantity and quality of the modelling inputs are crucial for wakes. In general, all models regarding their complexity require measurements of the turbulence level and the atmospheric stability condition. Additionally the wind speed and direction should contain information about their distribution so that a proper postprocessing can be performed and the results are fairly compared with the observations. Particularly for Ellipsys3D LES, when modelling the wind turbine, the tabulated values of the airfoil aerodynamic properties are required, which are calculated using the airfoil geometry. Such information is hard to obtain from the manufacturers. In addition, the methodologies used to obtain the lift and drag coefficients of a given geometry have their own inaccuracies and limitations.

Computational costs
In general, the computational expenses of the wake models increase with the complexity of the model. Therefore, the Jensen model is the fastest to produce results, followed by the Larsen model and Fuga. The DWM model needs a relatively highly resolved turbulence field to feed back the aeroelastic code in the current version, but yet the computational cost is not implied as a main issue. Ellipsys3D, on the other hand, suffers a lot from high CPU usage especially for the LES version, which eventually limits its application to the super-computers or clusters. There are a lot of studies regarding the hybrid RANS and LES methods which are more accurate and representative than RANS simulations but still more affordable than LES alone. A comprehensive review of various approaches to couple RANS with LES may be found in the study of Fröhlich and von Terzi [84].

Conclusions
Six of the wake models developed at DTU are investigated. The models have different levels of complexity, and overall they represent the wide range of wake models available for the wind energy industry and research community. The models are described and inter-compared using the Sexbierum onshore and Lillgrund offshore wind farms. Both benchmark cases have provided valuable insights in terms of the effects of the turbine spacing (or wind farm layout in general), wind direction averaging sector variations, turbulence intensity and possible atmospheric stability conditions. Finally, the models are briefly evaluated in terms of their application.
The benchmark cases show that the analytical and linearized models of DTU (the Jensen model, the Larsen model and Fuga) are convenient for large wind farm calculations as they are robust and computationally affordable. They provide good results both onshore and offshore implementations as long as the far wake region is considered and the atmospheric conditions are well defined.
The more sophisticated CFD solvers (Ellipsys3D RANS with k-ε and k-ε-f P turbulence closures and LES) are used in the benchmark cases. The k-ε-f P and LES in particular are observed to be in a very good agreement with the measurements. Because of their computational cost however, they are very rarely implemented on large wind farms and their applications are generally limited to the near wake region or highly complex flows.
The benchmarking study also shows that introduction of the wind direction uncertainty significantly improves the accuracy of the power predictions of the Jensen model, Larsen model and Fuga, for the Lillgrund case. For the Sexbierum case, however, even the state-of-the-art model Ellipsys3D LES fails to reproduce the depth of the wake deficit. The limited period of the investigated data and lack of information regarding the characteristics of the inflow are considered to be the reason of the model deficiencies, as they led to erroneous assumptions. Accordingly, the significance of the data set quality, as well as the quantity, for the wind turbine wake model benchmarking has to be underlined.