Atmospheric stability-dependent infinite wind-farm models and the wake-decay coefficient

We extend the inﬁnite wind-farm boundary-layer (IWFBL) model of Frandsen to take into account atmospheric static stability effects. This extended model is compared with the IWFBL model of Emeis and to the Park wake model used in Wind Atlas Analysis and Application Program (WAsP), which is computed for an inﬁnite wind farm. The models show similar behavior for the wind-speed reduction when accounting for a number of surface roughness lengths, turbine to turbine separations and wind speeds under neutral conditions. For a wide range of atmospheric stability and surface roughness length values, the extended IWFBL model of Frandsen shows a much higher wind-speed reduction dependency on atmospheric stability than on roughness length (roughness has been generally thought to have a major effect on the wind-speed reduction). We further adjust the wake-decay coefﬁcient of the Park wake model for an inﬁnite wind farm to match the wind-speed reduction estimated by the extended IWFBL model of Frandsen for different roughness lengths, turbine to turbine separations and atmospheric stability conditions. It is found that the WAsP-recommended values for the wake-decay coefﬁcient of the Park wake model are (i) larger than the adjusted values for a wide range of neutral to stable atmospheric stability conditions, a number of roughness lengths and turbine separations lower than (cid:2) 10 rotor diameters and (ii) too large compared with those obtained by a semiempirical formulation (relating the ratio of the friction to the hub-height free velocity) for all types of roughness and atmospheric stability conditions. Copyright © 2013 John Wiley & Sons, Ltd.


INTRODUCTION
Most wind Park wake (WPW) models are able to estimate wind-speed reductions within the wind farm for a wide range of wind speeds, assuming neutral atmospheric conditions in most cases. These models, such as the Park wake model 1 implemented in the Wind Atlas Analysis and Application Program (WAsP), 2 predict well the energy yield losses due to wakes when analyzing long terms of meteorological (met) and wind-farm data. This is partly because in a long term, most atmospheric static stability conditions at wind turbine sites are generally close to neutral. According to Troen and Petersen, 3 the long-term atmospheric stability is just a little biased to the stable side over land and to the unstable side over water from the analysis of met stations over Northern Europe, which showed a higher amount of positive compared with negative surface heat fluxes at few offshore met masts, whereas for Peña and Hahmann 4 is slightly biased to the stable side over the North Sea on the basis of the analysis of the probability distribution function of stability measures.
When analyzing wind-farm and met data from the Horns Rev I wind farm in the Danish North Sea, Jensen 5 estimated that the annual mean array efficiency reduces from 91.5% under unstable to 85.3% under stable atmospheric conditions. Barthelmie and Jensen 6 also estimated wind-farm efficiency reductions in stable compared with unstable atmospheric conditions up to 9% for the wind-speed range 9-10 m s 1 for the Nysted wind farm in the Danish Baltic Sea. Since wind-farm operators do not want to know the annual energy production of the wind farm only but would also like to forecast the wind-farm energy output for a given set of met conditions, which can rapidly change as shown in Vincent et al., 7 we need to run the WPW models for different met conditions, which include the state of the atmosphere. Inclusion of the atmospheric stability dependency in the WPW models is not straightforward. Therefore, an alternative is to simply adjust the parameters in the models to match the observed/measured data. Barthelmie and Jensen 6 found better agreement when comparing the Park wake model with the Nysted wind-farm data by using a lower wake-decay coefficient (0.03) than that recommended in WAsP for offshore conditions (0.05). Interestingly, the amount of stable atmospheric conditions are relatively large at Nysted 6 compared with those in the North Sea, 8 implying that the more stable the atmosphere, the lower the wake-decay coefficient for the Park wake model.
In this paper, we adjust the wake-decay coefficient of the Park wake model, evaluated for an infinite array of wind turbines, to match the wind-speed reduction estimated by the infinite wind-farm boundary-layer (IWFBL) model of Frandsen. 9 The adjustment can be carried out for different wind speeds, turbine to turbine separations, surface roughness lengths and atmospheric stability conditions. Since Frandsen developed his model for neutral conditions only, we extend it for diabatic atmospheric conditions (assuming horizontal and vertical homogeneity of stability) by using local atmospheric stability corrections to the logarithmic (log) wind profile and the resistance law constants of the geostrophic drag law. The stability corrections to the log wind profile are limited to the range of atmospheric conditions and heights where the surface-layer theory is valid. Therefore, considerations should be made, particularly under stable conditions, where the theory is limited to a few tens of meters above the surface only and where the boundary-layer height (BLH) is about the size of the turbines and therefore influences the shape of the wind profile as shown in Peña et al. 10 and Sathe et al. 8 and when the atmospheric stability changes within the wind farm.
There are other techniques to study the wind-speed reductions due to wind turbine wakes. 11 Computational fluid dynamics (CFD) methods have been extensively applied for multiple wind turbines and in the last couple of years, large eddy simulation (LES) methods have gained popularity compared with the Reynolds averaged Navier-Stokes (RANS) turbulence models. This is partly due to the extension of LESs to account for atmospheric stability conditions other than neutral, which allows LES results to be compared with benchmark cases. 12 The LES technique has also been used to study large arrays of wind farms, [12][13][14] which can also provide the information needed to adjust the wake-decay coefficient of the Park wake model. However, RANS/LES-based wake models are about 6-7 orders of magnitude slower than the WPW or IWFBL models.
In the description of the model by Frandsen and its extension, we compare it with the results for the wind-speed reduction of an infinite array of turbines of the IWFBL model of Emeis and Frandsen (E&F) 15 and those of the Park wake model. We also show the differences between the WAsP-recommended, the semiempirical (described in Section 3) and the IWFBL-adjusted wake-decay coefficients. Results from the Park wake model are not compared with wind-farm data. Analysis of such data is as challenging as the modeling itself, 16 particularly for different atmospheric stability conditions because since most wind farms have no means to estimate stability, it is difficult to separate the effect on wind-farm power production of stability, wind speed and turbulence from real data 6 and because in large wind farms, turbines do not operate concurrently and optimally all the time.

Review
The IWFBL model of Frandsen 9 assumes that within an infinite wind farm with the same turbine type and dimensions, two layers in the atmospheric boundary layer (ABL) are distinguished, one above and one below the turbines' hub height h as shown in Figure 1. At h, the shear of both layers is linked as where u 2 and u 1 are the friction velocities within the above and below layers, respectively, and t is the jump in shear due to the turbines. The latter is given as t D c t u 2 h , being the air density, u h the hub-height spatial average wind speed within the wind farm and c t the areal homogenously distributed thrust coefficient, where s r D x=D r and s f D y=D r , x and y being the along-wind and cross-wind turbine to turbine distances, respectively, C t the turbine's thrust coefficient and D r the rotor diameter.
Since the idea is to derive an expression for u h by using equation (1), Frandsen 9 applied the log wind profile to estimate u 1 and u 2 at hub height from the true (z o ) and the effective (z oo ) wind-farm roughness lengths, respectively. where Ä is the von Kármán constant ( 0:4). Frandsen 9 used a corrected version of the simplified geostrophic drag law of Jensen, 17 where G is the geostrophic wind speed and f p D f c exp.A /, f c being the Coriolis parameter and A a modified A parameter from the resistance-law constants, to derive an expression for the effective roughness as a function of the velocity scales: Equation (5) is then used to eliminate the u 2 dependency on z oo in equation (3, right), which results in For a simple derivation of u h from equation (1), Frandsen 9 defined K 2 D .1=Ä/ lnOEG=.hf p / and K 1 D .1=Ä/ ln.h=z o /; and thus, Substituting both forms in equation (7) into (1) leads to the solution for u h , i.e., the IWFBL model of Frandsen where The wind-speed reduction R u for the IWFBL model of Frandsen is then given as which is the ratio of u h to the asymptotic overall mean wind speed at hub height, i.e., u h .c t ¤ 0/=u h .c t D 0/. Equation (8) guarantees that u h .c t D 0/ is equal to u hfree , i.e., the undisturbed hub-height wind speed, as expected. It should be noted that by using any form of the geostrophic drag law, the ABL is assumed in this model to be barotropic and homogenous in terms of flow, atmospheric stability conditions and roughness (valid for large footprint areas). Emeis and Frandsen 15 further developed another IWFBL model on the basis of the mixing-length theory. Shortly, they assumed the specific turbulent downward momentum flux, = , to be driven by the vertical wind-speed gradient above the wind farm. Using a momentum exchange coefficient K m to parameterize this flux, E&F related it with the mixing length l as where u o is the undisturbed wind speed at a height z D h C z above the wind farm. This flux is assumed to be in balance with the momentum loss due to the turbines.
where c 0 t is the drag coefficient of the turbines. The surface drag coefficient is defined as c s D u 2 =u 2 h with u D u h Ä ln.h=z o / 1 (so it can be shown that c s D K 2 1 ). The effective drag coefficient, c teff , results from the combination of that of the surface and the one from the turbines, i.e., c teff D c 0 t C c s . In E&F, there are two reductions: the first is given by estimating the ratio u h =u o from equation 11; the second is found similarly but by replacing c 0 t by c teff . The wind-speed reduction equivalent to equation 9 is found from the ratio of the aforementioned two reductions.
Since it is straightforward to derive c t from the thrust curve of the wind turbines in the wind farm, we assume that c 0 t D c t in this study. When comparing equations (9) and (12), it can easily be demonstrated that they are identical if z=l D K 2 . So, one of their main differences is that the IWFBL model of Frandsen depends on u hfree through both G (via K 2 ) and C t (which for a wind turbine varies with wind speed), whereas the IWFBL model of E&F, equation (12), is only dependent on u hfree through C t provided that l is wind-speed independent.
For neutral atmospheric conditions, E&F suggested z=l 2=Ä, which in terms of the IWFBL model of Frandsen means lnOEG=.hf p / 2. Such approximation for a place in a rural area (z o D 0:025 m) or at an offshore location (z o D 0:0002 m), both assuming h D 70 m, u hfree D 10 m s 1 , A D 4:53 and latitude of 55.5 ı , is rather low, since lnOEG=.hf p / D 3:54 and 2:73, respectively * . Therefore, R u is generally lower for the IWFBL model of Frandsen than that for E&F. This, as shown later, is valid for a range of wind speeds and atmospheric stability conditions. The main 'drawbacks' of the IWFBL model of E&F are that (i) z is rather difficult to estimate and (ii) the mixinglength concept might be inappropriate when modeling wakes. 18 The approach by Frandsen also has drawbacks related to (i) the assumption that at some level above the hub height, the wind-farm wind speed approaches the geostrophic windspeed value and (ii) the value of A , since as a modified A parameter of the geostrophic drag law, it is rather uncertain and depends on atmospheric stability among others. 19 Figure 2 shows the wind-speed reduction comparison between the IWFBL model of Frandsen and that of E&F for different turbine to turbine separations, s, and roughness lengths † . As mentioned, the values for R u from the IWFBL model of E&F are larger (lower reductions) than those from the approach by Frandsen, following its behavior with turbine to turbine distance. There is a faster change in wind-speed reduction with separation, and the differences (between models) are larger within the range of turbine separations where most wind farms lie, i.e., s < 10.
The previous comparison is based on a fixed undisturbed hub-height wind speed and thus on a fixed C t value. When performing a similar analysis (not shown) but for a fixed turbine to turbine separation and surface roughness and different undisturbed hub-height wind speeds, which translate into a range of C t values, we find similar wind-speed reductions for both IWFBL models, being the R u values of the E&F model generally larger than those of Frandsen (Frandsen's reductions are only lower for very low u hfree values). For the analysis, we chose a C t curve from a 2 MW horizontal-axis wind turbine with a near-constant C t within the low-wind-speed range and decreasing values with increasing wind speed.

Atmospheric stability dependency
Emeis 20 extended the IWFBL model of E&F to account for diabatic atmospheric stability conditions. Since R u is a function of the mixing length l in equation (12), it was rather simple to extend l for diabatic conditions by using Monin-Obukhov similarity theory (MOST), 21 i.e., adding the dimensionless wind shear m , which is a function of the stability parameter z=L, where z is the height above surface and L the Obukhov length. L is a local stability measure, which is assumed in * G is estimated from the simplified geostrophic drag law, equation (4), and related to u hfree through the log wind profile similar to equation (3). † For simplicity, it is hereafter always assumed that s r D s f D s. MOST to be constant with height within the surface layer; in Emeis 20 and in this study, by correcting the wind shear by using a local L value, we extend its use beyond its original intend, since we further have to assume L to be constant with height above the surface layer and within the wind farm (see more details on these issues in the Section 4). Because of the different formulations of l, 22 this extension can be carried out in several ways, and Emeis 20 tried different forms including l D Äz 1 m , which was proven to be valid within the surface layer when compared with turbulence measurements and spectra. 23 Since z=l and K 2 are equivalents, we can also extend equations (8) and (9) by using m .h=L/. Further, as Emeis did it for c s , K 1 has to be extended for stability conditions by means of the stability correction of the log wind profile m , which comes from the integration of m , so m is also a function of z=L. u h , for example, becomes However, there is an inconsistency in equation (13) because u h should equal u hfree when c t D 0. They are indeed equal under neutral conditions because m D 1 and m D 0 for h=L D 0. When h=L ¤ 0, on the other hand, u h .c t D 0/ from equation (13) becomes h-dependent and does not approach the u hfree value.
One way to avoid such 'inconsistency' is by including the stability correction function ( m ) to the log wind profiles within the two layers in equation (3), so that within the wind farm and at hub height, and In the absence of a local stability measure at hub height within the wind farm, it is assumed that L, in equations 14 and 15, is equal to a stability input given upstream the wind farm. One can work out the derivation of u h , in a similar way as that explained in Section 2.1 (see steps between equations 4-6) and finds (by inserting z oo in equation (5) into (15)).
Therefore, to solve u h in equation (1), K 1 and K 2 have to be redefined as .1=Ä/ OEln.h=z o / m .h=L/ and .1=Ä/ lnOEG=.hf p / C m .h=L/ , respectively. The expressions for u h and R u are then formally the same as those in equations (8) and (9), also valid for neutral conditions. However, A (included in the modified Coriolis parameter f p in equation 4) also depends on the state of the atmosphere, as this is close related to the resistance law constants A and B, which have been found to vary with stability. This dependency is classically described by extending A and B to be functions of the dimensionless stability parameter o D Äu =.f c L/. Despite disagreement among researchers on the forms of A. o / and B. o /, due to the variety of data used in the different studies and their high uncertainty, it is possible to get empirical formulations for both. 24 On this basis, the dependency of A on atmospheric stability, A . o /, can be derived (see Appendix A for details).
For a given u hfree value and local stability condition (an L value at a given height upstream the wind farm), G can then iteratively be estimated using the simplified geostrophic drag law, equation (4) (now including the A . o / stability dependence) and the upstream-undisturbed friction velocity u free derived from the diabatic wind profile, Now, it is important to notice that in the non-neutral cases, K 2 depends on u through A : Thus, R u in equation (9) must be modified as where u 2 is found iteratively from the simplified geostrophic drag law but by using the expression for z oo for non-neutral conditions (see Appendix B for its derivation).
This procedure ensures that u h , in equation (8) evaluated for c t D 0, gives u hfree . Figure 3 shows the extended IWFBL model of Frandsen for a range of dimensionless atmospheric stabilities h=L, a number of roughness lengths and a fixed dimensionless turbine to turbine distance. For the range of atmospheric stabilities, R u always increases with high roughness length values (similar to Figure 2). The striking result is that the variation of R u is much higher for this broad range of atmospheric stabilities than for the rather wide range of roughness lengths, with general low-wind-speed and high-wind-speed reductions under unstable and stable atmospheric conditions, respectively. At a given positive value of h=L (dependent on the given z o value), R u starts to increase with increasing stable conditions. This is because u hfree is used as the only wind-speed input parameter to the model and for this case (h D 70 m), and for stable and very stable conditions, such a height is comparable to the BLH that can indeed be below 50-70 m in stable and very stable conditions. MOST (strictly valid for 10% of the BLH only) then predicts a too high and unrealistic stability correction when estimating u free from u hfree , which is used for computing G in a range where the stability corrections to the A and B functions and, thus, A in the simplified geostrophic drag law is less applicable. A way to avoid the increasing R u values for stable conditions is by using, as an input parameter for the model, a wind speed not at hub height but at a height well inside the surface boundary layer, where MOST is well applicable.

Model behavior
We can also compare the results for the model of Emeis, i.e., the atmospheric stability-extended IWFBL model of E&F, and the extended IWFBL model of Frandsen for different atmospheric stability conditions and dimensionless turbine to turbine separations. This is illustrated in Figure 4 for an offshore site, where the models have the same behavior with turbine separation distances and atmospheric stability. As for the previous comparisons, the Frandsen's-type model shows larger wind-speed reductions than those from the Emeis'-type model for neutral conditions. For the computed stable and unstable conditions (L D 100 m and L D 50 m, respectively), the reductions are larger for the Emeis'-type model. This is partly due to the reference height used for the free wind speed; when using a low value for h, R u from Emeis' model becomes higher than that of the extended model of Frandsen under stable conditions.

Wind profiles
We can also compute the spatial average wind profiles within the wind farm (by assuming L to be vertically and horizontally constant), which can be used to analyze the effect of the wind farm in the ABL (not only at hub height) and to compare with the LES results (e.g., those in Lu and Porté-Agel, 12 and Calaf et al., 14 ) by using similar functions as those in equations (14) and (15), first by estimating u free for the different stability conditions, as in equation (17), which is then used for the estimation of G (by using a stability-corrected A value). z oo is directly computed from equation (20). They become and  It is important to notice that both layers (above and below h) are modeled using the approach of MOST. As h increases, e.g., with larger turbines, the above-h layer will be affected by the BLH and the baroclinicity (among others), and the MOST approach is less valid. Also, the stable corrections from MOST predict very large wind shears for the above-h layer, which lead to wind speeds u.s ¤ 0/ above the undisturbed ones u.s D 1/ already at 200 m in Figure 5-bottom.

INFINITE PARK WAKE MODEL
The Park wake model implemented in WAsP is based on the momentum wake model by N. O. Jensen (hereafter known as NOJ), * 25 where the wind speed immediately before the first wake-affected turbine, u 1 , can be estimated as ( Figure 6) where u free is the upstream undisturbed wind speed, a the induction factor (a D 1 p 1 C t ), k w the momentum entrainment or wake-decay coefficient and r o the initial wake radius behind the rotor. In the Park wake model, r o is equal to the turbine's rotor radius r r . Frandsen 9 proposed, by semiempirical means, k w 0:5= ln.h=z o /, which suggests that k w is related to the atmospheric turbulence characteristics because 0:5= ln.h=z o / u free =u hfree under neutral stability conditions. * Although commonly referred to as a momentum conservation model, NOJ model conserves mass and not momentum within the control volume.   25 The wind speed immediately before a wake-affected turbine u 1 is estimated as a function of the distance to the upstream turbine x with rotor radius r r , the free wind speed u free , and the thrust coefficient and the wake-decay coefficient k w .
Therefore, the wind-speed reduction in the Park wake model is also a function of k w . The ground interaction of a wake is modeled by adding the wake reflected at the surface, i.e., a wake seemingly originating from a reflected 'underground rotor'. The surface interaction is then modeled by considering the combined effect of the direct and the reflected wakes. Further, the efficiency of the wind turbine cluster is estimated by combining the effects of four types of overlapping wakes 1 : 1. From directly upwind rotors (NOJs original approach). 2. Reflected 'underground rotors' directly upwind. 3. Shading rotors upwind but left or right to the wind direction. 4. Reflected shading 'underground rotors' upwind but left or right to the wind direction.
The combined effect of two or more overlapping wakes on a downwind rotor was modeled through an empirical quadratic summation rule. 1 Here, the partial wake overlap with a rotor is considered by applying the overlap fraction to the speed deficits of the individual wake. Thus, the rotor diameter is required as an input to the model. Rathmann et al. 26 analytically solved the contribution of the four types of wakes (described earlier) for an asymptotically infinite number of wind turbines. This is described here as the 'infinite Park wake model' (IPW) (Appendix C).

Adjusted wake-decay coefficient for the infinite Park wake model
When using the Park wake model, we are required of the wake-decay coefficient for the particular site. To a first approximation, it is common procedure and recommended to use k w D 0:050 and k w D 0:075 when performing wake analysis in WAsP over sea and land surfaces, respectively. However, k w should be estimated for the particular z o value and atmospheric stability condition of the site, since k w u free =u hfree D Ä=OEln.h=z o / m .h=L/, where the correction due to atmospheric stability m is here already considered from that form proposed by Frandsen. 9 This already results in much lower k w values than those recommended in WAsP for a wide range of roughness values and stabilities; k w is lower for low compared with high roughness values and higher for unstable compared with stable conditions. Figure 7 compares the wind-speed reduction from the IPW model-when using two types of k w -values: from the form k w D u free =u hfree , as shown earlier (so we again assume L to be invariant within the wind farm), and the WAsP-k w values-with that from the IWFBL model by Frandsen for neutral conditions, different dimensionless turbine to turbine separations and a number of roughness lengths. The IPW model shows the same behavior as the IWFBL model, i.e., the higher the wind-speed reductions, the smoother the terrain and the shortest the turbine to turbine separation. The IPW model using the WAsP-k w values shows similar reductions compared with the IWFBL model, slightly higher ones for offshore conditions. The IPW model using k w D u free =u hfree shows much higher reductions compared with the IWFBL model for all roughness lengths and similar ones to the IPW model using the WAsP-k w value for z o D 0:2 m only. This result already illustrates that the recommended WAsP-k w values seem to be too high for wind farms located on flat and homogeneous terrain, except for areas with high roughness values such as forests.
When the IPW model is adjusted so that it gives the same wind-speed reductions as those of the IWFBL model of Frandsen, k w has to be modified/adjusted, and thus, it becomes dependent not only on z o (k w increases with increasing z o values) but also on turbine separation as shown in Figure 8. The 'adjusted' k w values are generally much lower for smaller than for larger turbine separations and lower than the WAsP-k w values for z o D 0:002-0:02 m and s < 10. The lowest k w values are given by the form k w D u free =u hfree , which does not take the turbine separation into account. In a similar fashion, a constant turbine to turbine distance, e.g., s D 7, can be used to analyze the behavior of the windspeed reduction R u for different atmospheric stability conditions and roughness lengths. This is illustrated in Figure 9. The wind-speed reductions for the IPW model using the form k w D u f ree =u hfree are generally higher than those for the extended IWFBL model by Frandsen, except for z o D 0:2 m and a range of neutral to stable conditions. The results of the extended IWFBL model by Frandsen show an increase of R u at a given positive h=L value because, as mentioned earlier, the reference wind speed is too high and MOST might be not longer valid. R u does not vary with atmospheric stability in the IPW model as strongly (stability is accounted for through k w only) as it does in the extended IWFBL model of Frandsen. The IPW model using the WAsP-k w values shows similar reductions compared with those from the extended IWFBL model of Frandsen only for atmospheric conditions close to neutral. The correspondent adjusted wake-decay coefficients, found when matching the wind-speed reduction of the IPW model to that of the extended IWFBL model of Frandsen for the different atmospheric stability conditions in Figure 9, are illustrated in Figure 10. The behavior of k w with stability when assuming k w D u free =u hfree is similar to that found with the adjusted wake-decay coefficients for all roughness lengths, i.e., a decrease of the wake-decay coefficient the more stable the atmosphere, and generally shows the lowest k w values (except when compared with the adjusted values for a narrow range of stable conditions and z o D 0:2 m). The adjusted values tend to be much lower than the WAsP-recommended ones for a wide range of neutral to stable conditions and all roughness lengths. The increase of the adjusted k w is due to the increase in R u at the same stable h=L range, a behavior which disappears when the reference undisturbed wind speed is taken at a height within the surface boundary layer.

DISCUSSION
In this paper, k w in the Park wake model for the infinite wind farm is, among others, adjusted so that the model matches the value of the wind-speed reduction of the extended IWFBL model of Frandsen. In a similar fashion, the adjustment can be performed using the IWFBL model of Emeis, which also shows that the WAsP-recommended values for k w are generally large for land-type roughness and short turbine separations. There are other WPW models, which use the k w parameter, and here, we choose the Park wake model because it considers (among others) the effect of adjacent turbines, is a commonly used model and is the base of wind power calculations in WAsP. However, direct comparison and matching of the models are not totally fair because the wind-speed reduction estimated from the IWFBL models is related to the spatial average wind speed within the wind farm, whereas that from the IPW model is related to the wind speed immediately before the last turbine.
From our findings, we do not suggest to use the results for the wake-decay coefficient when performing wind power calculations directly, since the analysis is performed on the basis of the assumption either that k w D u free =u hfree or that of an infinite wind farm and in reality, only a couple of wind farms might be treated as such. For a finite but large wind farm, the first wake-affected turbines experience a rather different wake regime than those at the final rows perpendicular to the mean undisturbed wind direction and, therefore, k w should be slightly decreased with downstream distance to values similar to ours under close turbine separations and near-stable conditions. An interesting and similar approach was shown by Rathmann et al. 26 , in which after the k w of the IPW model is adjusted to the IWFBL model of Frandsen, the WAsP-recommended k w value in the Park model is gradually modified to asymptotically approach the infinite k w value for 'deep' positions in the wind farm via a relaxation constant. Their approach agreed better than only assuming the WAsP-recommended values when compared with data from the Horns Rev and Nysted wind farms.
As wind turbines become larger, the IWFBL models have to be extended to account for effects that are normally not considered for small-medium turbines, which are well inside the surface boundary layer, such as the BLH. As illustrated in Figure 5, particularly for stable conditions, the wind speed predicted using MOST highly and sometimes unrealistically increases with height, and the addition of the BLH as a wind profile parameter will damp such unrealistic growth, as shown in the wind profile models by Peña et al. 10 over the sea, and Peña et al. 22 and Gryning et al. 27 over land. Unfortunately, from our knowledge, there is still a lack of BLH data for wind energy purposes, which can be used for model comparison.
Frandsen 28 attempted to compare theoretical wind-farm efficiencies by assuming that the change in roughness imposed by the infinite wind farm can be modeled as a smooth-to-rough generated internal boundary layer (IBL) to those derived from 'traditional' IBL models such as that of Miyake, 29 also used in WAsP for natural roughness transitions, finding similar results for an infinite row of turbines. However, as shown in Floors et al., 30 although the IBL models give good results for neutral atmospheric conditions, observations of the IBL height show a strong dependency on atmospheric stability, which is normally neglected in the IBL models because it is generally believed that the mechanical contribution dominates the IBL growth. Since from our results the efficiency of the wind farm strongly depends on atmospheric stability, the question arises on how the IBL models can be adjusted for wind-farm efficiency analysis under different atmospheric conditions. Wind speeds within and downstream the wind farm are therefore needed to validate, for example, how accurate equation (20) describes the wind-farm characteristics by applying it to roughness models such as that in WAsP.
As previously mentioned, our results are difficult to compare with observations, since their analysis is at least as challenging as the modeling itself. Gaumond 16 found, for example, that the degree of accuracy of wake models evaluated for large offshore wind farms is much more dependent on the way data and simulations are post-processed than the physics of the models itself. There are, however, studies intending to reproduce the effects of large arrays of turbines on boundary layers by using LES techniques. Calaf et al. 14 'immersed' medium-size and large-size arrays of turbines in a neutral boundary layer, and for all the different cases (varying the geometrical loading, aspect ratios s r =s f and roughness values), the LES-based average velocity profiles as stated in Calaf et al. 14 'clearly showed the existence of two log laws above and below the turbine region', verifying the fundamental assumption of the IWFBL model of Frandsen. From their LES results, they proposed a new formulation of z oo on the basis of that of Frandsen et al. 31 As their analysis of z oo , our proposed expression in equation (20), is not very sensitive to the background roughness z o but to atmospheric stability, which can be incorporated to those expressions in Calaf et al., 14 and Meyers and Meneveau. 13 Lu and Porté-Agel, 12 on the other hand, studied an infinite large wind farm by horizontally applying periodic boundary conditions in a LES of the stable boundary layer covering a wind turbine. The shapes of their average velocity profiles are in good agreement with those in Calaf et al., 14 and they found that the stable BLH considerably increased with simulation time, attributing this effect to radial and vertical wake expansions. We certainly believe that this is partly the case, but the BLH also increases because for all stability conditions, the above hub-height friction velocity, u 2 , is higher than u free , which ends up lifting the BLH (in our model u 2 is about twice the value of u free for stable conditions).
Here, we assume the local 'undisturbed' hub-height value of L to be representative of the upstream atmospheric stability conditions within the surface layer and higher (in cases where h is above it) and within the infinite wind farm. The local stability conditions inside wind farms actually vary because of the presence of the turbines, which change both the local momentum and buoyancy fluxes (L is related to the ratio between both) as Zhang et al. 32 showed from wind-tunnel experiments. Therefore, the vertical and horizontal stability homogeneity assumption does not probably hold, but this needs to be further investigated as there is a lack of experimental data to verify it. The models here presented predict the change of momentum flux, but there is no formulation for the buoyancy flux. Lu and Porté-Agel 12

CONCLUSIONS
The Park wake model has been used to derive the wind-speed reduction asymptotically reached by an infinite array of wind turbines and compared with the wind-speed reduction of the boundary-layer model of Frandsen for the infinite wind farm. The models show the same behavior for different dimensionless turbine to turbine separations and surface roughness lengths; the higher the roughness length and the longer the turbine separation, the lower the wind-speed reduction. The IPW model generally shows the highest wind-speed reductions, by assuming the form k w D u free =u hfree compared with those of the IWFBL model. The IWFBL model of Frandsen has been extended to account for atmospheric static stability conditions (by assuming horizontal and vertical stability homogeneity) and shows a similar behavior with stability when compared with the IWFBL model of Emeis (the stability extended model of E&F); the more stable the atmosphere, the higher the windspeed reduction, being the relative wind-speed reduction due to stability much higher than that due to surface roughness lengths (a wide range of roughness lengths has been considered). Emeis' model estimates lower wind-speed reductions than Frandsen's model for neutral conditions only and a wide range of wind speeds.
The extended IWFBL model of Frandsen shows that to get an increasing wind-speed reduction when going from neutral to stable conditions, it is required that the inputs to the model (i.e., the reference wind speed and atmospheric stability) are given close to the surface where MOST is fully valid. This results in a good estimation of the stability-corrected undisturbed friction velocity, which in turn provides a reasonable estimation of the geostrophic wind speed by using the stability-corrected A value.
The IWFBL model by Frandsen is the only one of the models here studied that is dependent on wind speed not only through the C t curve of the turbine but also because it is geostrophic wind speed dependent. Thus, for a nearly constant C t value, the wind-speed reduction increases anyway. However, at this stage, it is difficult to judge the quality of the models, since the IWFBL models, for example, either assume, among others, a relation between the mixing length and the height where the wind is undisturbed or a relation between the geostrophic and the hub-height wind-farm wind speed. Also, we do not validate the models against wind-farm data, but we expect to do so in the near future.
Finally, and most importantly, we are able to calibrate the 'more realistic' and more complex Park wake model to match the wind-speed reduction of the extended IWFBL model of Frandsen for different wind speeds, roughness values, turbine separations, turbine characteristics and atmospheric stability conditions by adjusting the value of the wake-decay coefficient, considering an infinite array of wind turbines. Already for neutral conditions, the k w values are generally lower than those normally recommended by WAsP for wind power calculations for a number of surface roughness lengths and turbine separations lower than 10 rotor diameters. They are also strongly dependent on atmospheric stability, decreasing as the atmosphere becomes more stable, higher than the WAsP-recommended values for a number of roughness lengths under very unstable atmospheric conditions and much lower than those WAsP-recommended under neutral and a wide range of stable atmospheric conditions.
In the original formulation of Jensen et al., 24 o D u =.f c L/. However, we choose to include Ä in o , as in Zilitinkevich, 19 and Long and Guffey, 34 because the points used by Jensen et al. 24 to derive equations (25)- (27) were taken from the study of the Wangara data by Clarke and Hess, 35 where o D Äu =.f c L/.
For a range of dimensionless stabilities, G is computed from equation (24) by using the forms in equations (25)- (27). The result can then be used to estimate A from rearranging the simplified geostrophic drag law, in equation (4). Figure 11 illustrates the results of such analysis, where an analytical form for A . o / is also given.

APPENDIX B: EFFECTIVE WIND-FARM ROUGHNESS LENGTH
From equation (15), one can easily derive where the ratio u h =u 2 can be replaced by 1= p .u 1 =u h / 2 C c t from equation (1) 31 deduced a similar form for neutral conditions only). It should be noted that the dependency on atmospheric stability of z oo is through m .h=L/ and K 1 , which also depends on m .h=L/.

APPENDIX C: THE INFINITE PARK WAKE MODEL
The total wake deficit ı T D 1 R u of the Park wake model implemented in WAsP is estimated as the quadratic sum of four types of overlapping wakes: The first term corresponds to the wakes directly upwind: .1 C 2k w s r j / 4 ; where ı o is the initial wake deficit, ı o D 1 p 1 C t , d w is a dimensionless wake diameter, d w D D w =D r , with D w as the wake diameter, D w D D r .1 C 2k w s r /, and j is the number of rows upwind from the considered turbine.
The second term corresponds to the reflected 'underground rotors' directly upwind: where m is the number of minimum rows to the first upwind row from where this type of wake has an effect on the considered turbine; in this case, m II D .2h=D r 0:5/=.k w s r /.
The third term corresponds to the shading rotors upwind but left or right to the wind direction: where N III is the number of turbines to the left and to the right of the wind direction that are able to throw their wake onto the considered turbine, N III .s j / D d w .s j /=.2s f / and m III D .s f 0:5/=.k w s r /, so equation (35) becomes The fourth term corresponds to the reflected shading 'underground rotors' upwind but left or right to the wind direction: where N IV .s j / D q d w .s j / 2 .4h=D r / 2 =.2s f / is the number of turbines to the left and to the right of the wind direction that are able to throw a reflected wake onto the considered turbine and m IV D q s 2 f C .2h=D r / 2 0:5 Á =.k w s r /, so equation (37) becomes When using these last four results, one assumes that C t is constant throughout the array, i.e., it does not strongly vary within the range of reduced wind speeds.