Articles | Volume 7, issue 2
Research article
04 Apr 2022
Research article |  | 04 Apr 2022

RANS modeling of a single wind turbine wake in the unstable surface layer

Mads Baungaard, Maarten Paul van der Laan, and Mark Kelly

Unstable atmospheric conditions are often observed during the daytime over land and for significant periods offshore and are hence relevant for wake studies. A simple kε RANS turbulence model for simulation of wind turbine wakes in the unstable surface layer is presented, which is based on Monin–Obukhov similarity theory (MOST). The turbulence model parametrizes buoyant production of turbulent kinetic energy (TKE) without the use of an active temperature equation, and flow balance is ensured throughout the domain by modifications of the turbulence transport equations. Large eddy simulations and experimental data from the literature are used for validation of the model.

1 Introduction

Wind turbine wakes have been studied for decades using many different methodologies, including wind tunnels, field experiments, analytical engineering models and numerical simulations. A review of these methodologies is given by Porté-Agel et al. (2020), and it is noteworthy that many of the references therein are from the past decade. The motivation for many of these new studies is the large number of new wind farms emerging each year, where wake effects significantly impact the annual energy production (AEP), as well as wind farm lifetime through increased fatigue.

A sub-category of numerical simulations is the Reynolds-averaged Navier–Stokes (RANS) approach, which is a computational fluid dynamics (CFD) method that solves for the mean fields. This means that no time history of the flow is obtained; however, the computational resources required for RANS are very small compared to higher-fidelity CFD methods, making RANS an attractive option for parametric studies or for isolating various physical effects (see van der Laan et al.2021). The wind turbine forces are commonly represented as actuator disks (ADs) in RANS; several types of AD models are reviewed by van der Laan et al. (2015a). Compared to engineering models, an advantage of RANS is that physical features of the flow (e.g., induction zones, wake interaction, shear layers, ground effects and flow over complex geometries) are solved for directly rather than being prescribed through empirical relations. Disadvantages are that fatigue loading cannot be determined due to the steady nature of the method and that the solution relies heavily on the turbulence model.

The part of the atmosphere closest to the ground, i.e., the atmospheric surface layer (ASL), can be parametrized with the similarity theory of Monin and Obukhov (1954) (MOST) and used as inflow for RANS simulations of wind turbine wakes. The kε turbulence model is usually preferred in RANS wake studies, and Crespo et al. (1985) for example utilized this (although in a parabolized RANS setup, which requires less computational resources but is less accurate) to simulate a single wake in stable, neutral and unstable conditions. The wake was found to recover faster (i.e., approach the freestream velocity at a shorter downstream distance) in the unstable ASL and slower (i.e., approach the freestream velocity at a longer downstream distance) in the stable ASL compared to in a neutral ASL, and this was later confirmed in field experiments by Magnusson and Smedman (1994) and full RANS simulations, including the temperature equation, by Alinot and Masson (2002). Rados et al. (2009) added a parametrized buoyancy term to the kε equations based on the MOST expressions, eliminating the need for a temperature equation. The “indirect” method of Rados et al. (2009) was shown by El-Askary et al. (2017) to produce similar wake deficit and turbulence intensity (TI) profiles as the “direct” method that employs a temperature equation.

In all the RANS studies discussed thus far, the simulations suffer from a known imbalance in the k and ε equations; this means that the freestream velocity and turbulence profiles vary horizontally throughout the domain, so that different wake results will be obtained depending on the streamwise position of the simulated turbine. van der Laan et al. (2017) solved this problem via the indirect method by adding analytical terms to the equations to be consistent with the ideal turbulent kinetic energy (TKE) budget under MOST and thus enforce a mean balance at all points. Han et al. (2019) used this approach in the direct method but did not show the extent to which their model is in balance.

Although there seems to be a general consensus that wakes should recover faster in unstable conditions, field measurements by Hansen et al. (2012) and Machefaux et al. (2016) found similar wake deficits for unstable and neutral conditions. This can possibly be attributed to the large uncertainties inherent in such measurement campaigns arising from sensors, post-processing and the unpredictable inflow provided by nature. In contrast, large eddy simulation (LES) offers a controlled environment, where complete statistics of all field variables can be extracted but at a large computational cost compared to RANS. Examples include Churchfield et al. (2012), Abkar and Porté-Agel (2015), Ghaisas et al. (2017) and Xie and Archer (2017), which simulate wakes in both unstable and neutral conditions for a wide variety of cases. All these studies agree with the general consensus and explain it with the increased TI encountered in unstable conditions due to buoyant production of turbulence. Nevertheless, both Alinot and Masson (2002) and Keck et al. (2014) show that a faster wake recovery in unstable conditions is still observed, even when the reference TI is kept fixed (by changing the roughness length); they argue that the enhanced wake recovery must be caused by the increased turbulent length scale associated with the unstable conditions because the turbulent velocity scale is approximately constant for fixed reference TI and wind speed.

The balanced kε MOST model by van der Laan et al. (2017) can be combined with the fP correction, which was originally formulated by Apsley and Castro (1997) and later used by van der Laan (2014) to circumvent the over-diffusiveness of the kε model in the wake region. However, the fP limiter was derived and calibrated for a neutral ASL, where it has been applied in many cases with success, but for non-neutral conditions it has yielded unphysical behavior, especially in unstable cases (van der Laan et al.2021). Modifications to the MOST kεfP equations in the unstable regime are therefore suggested in this paper and validated against various field experiments and LESs.

2 Simulation setup

The wakes are simulated with the incompressible, finite-volume flow solver EllipSys3D (Michelsen1992; Sørensen1995). After continuous development since then, it is now a highly scalable code, which can be run in parallel on large high-performance clusters (HPCs) via a message passing interface (MPI). Thus a typical RANS simulation of a single wake only takes a few minutes to simulate on a contemporary HPC, while a similar case with LES would take several hours to run, even with an order of magnitude more computer resources available. In terms of CPU hours, van der Laan et al. (2015b) estimated LES to be approximately 103 times more expensive than RANS, and this estimate may even be considered conservative because the LES inflow was created using a Mann-model turbulence box and not with the more expensive precursor method. Furthermore, several LES runs are in principle necessary in order to create an ensemble average, which multiplies the cost of LES. This clearly motivates the development of the RANS model as a fast, albeit less accurate, alternative to LES.

The different components of the RANS simulations will be discussed in the following sections.

2.1 Inflow profile for unstable ASL

Numerous articles have been written about MOST, and a historical review is given by Foken (2006). The theory is expressed and applied via the dimensionless stability parameter:

(1) ζ z L ,

where z is the height above ground and L is the Obukhov length. Negative ζ corresponds to unstable conditions, while ζ=0 corresponds to the neutral limit where there is no effect of buoyancy. Neutral conditions are typically defined as |L|-10.002 m−1 (e.g., Gryning et al.2007) and tend to occur most often, with observed distributions of the stability (1/L or ζ) having a peak around zero (Kelly and Gryning2010). The most common unstable Obukhov lengths occur at −0.02 m-1L-1-0.002 m−1 (Kelly and Gryning2010), but offshore there tends to be a bias towards more unstable conditions, i.e., more negative L−1 compared to onshore (Sathe et al.2013). Various parametrizations have been suggested for wind speed, k and ε in terms of ζ; in this paper we use the widely accepted forms of Dyer (1974) for U (namely the Ψm and Φm functions) and those found in Kaimal and Finnigan (1994) for ε and k (see van der Laan et al.2017, for more details). Under unstable conditions these are





(7) Φ ε = 1 - ζ .

The above relations are valid for -2ζ<0, so for a fixed L<0, it means that the equations are in principle only valid up to z-2L, where the free convection regime starts (i.e., buoyant production dominates over shear production of TKE). Although the blade tip of a modern turbine can reach beyond −2L in unstable conditions (e.g., for ztip=200 m this happens when L-1-0.01 m−1, which is not rare), we nevertheless still choose to apply the profiles and in fact use them all the way up to the top boundary. More realistic inflow profiles for RANS covering the whole atmospheric boundary layer (ABL) are indeed a current research topic but will not be discussed further in this paper. Maronga and Reuder (2017) reason that MOST is a pragmatic solution because the parameters needed for more realistic inflow profiles are often not available.

The roughness length z0 and friction velocity u* in Eqs. (2) to (7) can be set using reference values (i.e., defined at z=zref) of wind speed (Uref) and total TI (Iref) along with the stability parameter (ζref):


Note that TI (I) here is not the same as streamwise turbulence intensity (σu/U); in this paper “TI” will refer to total TI (i.e., I23k/U), unless stated otherwise. A typographical (sign) error has been corrected in Eq. (8), compared to the similar expression found in van der Laan et al. (2017). The von Karman constant κ and Cμ parameter are given with the other constants in Table 1.

Table 1Parameters of the kε MOST turbulence model; see Sørensen (1995).

Download Print Version | Download XLSX

Examples of four inflow profiles with identical Uref are shown in Fig. 1. The stability and TI differ among the cases, but they still have approximately the same averaged power density; i.e., 1A12ρU(z)3dA308.9 W m−2 (±0.6 %). One peculiarity of the unstable MOST profiles is that the TI does not go to zero for z→∞; this is connected to the ASL assumptions used by MOST (i.e., the mixed and upper layers above the ASL are lacking, where I becomes constant and then vanishes). Another peculiarity (for both neutral and unstable conditions) is that higher TI leads to larger shear (dU/dz) because the velocity gradient scales with u*, which scales with Iref (Eq. 9); this is a consequence of specifying both TI and hub height velocity.

Figure 1Analytical MOST profiles. Combinations of low/high TI and neutral/unstable stability. The rotor area of a NREL5MW turbine is shown. Dashed lines are used for ζ to make all profiles visible.


The eddy viscosity profile, νt(z)=Cμk2ε, is especially interesting to compare between the cases because νt features as a diffusion coefficient in the Reynolds-averaged momentum equation and therefore is connected with the entrainment of ambient air into the wake. A faster wake recovery is therefore expected for larger νt, which as seen in Fig. 1 can be obtained by increasing either the turbulence strength (Iref) or how unstable the atmosphere is (ζref) or both.

The eddy viscosity is sometimes expressed as a product of turbulent velocity and length scales (see Pope2000): νt=utt, where ut=Cμ1/4k1/2 and t=Cμ3/4k3/2ε-1. These are plotted in Fig. 2 from which one can see that Iref only affects ut, while t is unchanged. Changing ζref mainly alters t, while the hub height ut and rotor-averaged ut are unchanged. The increased νt (which gives faster wake recovery) due to unstable conditions can therefore not only be attributed to a larger Iref, but also to an increased turbulent length scale caused solely by ζref.

Figure 2Turbulence scales in the freestream. Same labels as in Fig. 1.


2.2 Wind turbine representation

A recently developed actuator disk (AD) model by Sørensen et al. (2020) is utilized in this paper. The model can be derived from conservation of energy (Bernoulli's equation), conservation of angular momentum (Euler's turbine equation) and an analytical expression for the near-wake azimuthal velocity distribution. The latter is modeled by a vortex extending from the center of the AD to infinity with constant circulation; hence, it resembles the classical Joukowsky optimum rotor (see Okulov and Sørensen2010), and the AD model is therefore referred to as the “Joukowsky AD”. A summary of the model formulation is given in Appendix A2.

The main advantage of the Joukowsky AD over the widely used “airfoil AD” (e.g., Sørensen and Kock1995; Porté-Agel et al.2011; van der Laan et al.2015a) is that only a few parameters are necessary: the thrust coefficient CT, tip-speed ratio λ, rotor radius R and freestream reference wind speed Uref (in addition to these, the power coefficient CP is also made an input parameter in our implementation as described by van der Laan et al.2020). Nevertheless, it is still able to model non-axisymmetric force distributions and wake rotation, similar to the disk loading of an airfoil AD. Porté-Agel et al. (2011) argued that these are important features to capture the correct wake behavior in the near wake, while van der Laan et al. (2015c) showed that they are only of minor importance for the far wake. Wake deficit and rotor loading of the Joukowsky AD have been found to compare well with several validation cases conducted by Sørensen et al. (2020) and Sørensen and Andersen (2020). This is verified to also be the case for our RANS simulations in Appendix A2.

No nacelle or tower is included in our simulations, which have been shown to be a good approximation for >3D downstream of the turbine according to Kasmi and Masson (2008) and Li et al. (2020).

2.3 RANS

A homogeneous, flat lower surface is assumed for all cases in this paper. The inner part of the mesh surrounding the AD is called the “wake domain” and is shown for a typical case in Fig. 3. In this area, a horizontal resolution of Δx=Δy=D/10 is used (based on the grid study in Appendix A1), while grid stretching is used in the vertical direction with Δz=z0 at the first cell and Δz=D/10 at the cell at z/D=3. The wake domain is however only a small part of the full domain: the full domain extends an additional xin=5 km to the west, north, east and south, respectively, while the top of the grid is at z/D=25. Grid stretching is used in all directions outside of the wake domain to circumvent an excessively large number of cells (the case shown in Fig. 3 has 0.45×106 cells in the wake domain and 2.1×106 cells in total). The choice of having grids with size on the order of ∼10 km is made to avoid tunnel-like blockage effects and to have fully developed inflow profiles at the turbine position, while the drawback is that the majority of the cells are actually outside the region of interest, i.e., the wake domain.

Figure 3Top and side views of the wake domain, the size of which is {lx,ly,lz}/D={16,4,3}. The total grid size is {Lx,Ly,Lz}={lx+2xin,ly+2xin,lz+22D} and is too large to be shown here because xin=5 km. Only every fourth cell is plotted.


The numerical solution strategy of the incompressible RANS equations in EllipSys3D is thoroughly discussed in other publications (Michelsen1992; Sørensen1995; Sørensen et al.2007), so only the main features are discussed here. The SIMPLE method is used with a modified Rhie–Chow algorithm, following Réthoré (2009) and Troldborg et al. (2015), to avoid the numerical wiggles induced by the discrete actuator disk body forces.

As mentioned in the introduction, the flow variables in an empty domain with MOST inflow can be kept in balance by modifying the k and ε equations, as suggested by van der Laan et al. (2017):




The source term, Sk, and the Cε,3 parameter constitute the two modifications compared to the usual kε equations (similar corrections exist for the stable ASL but are not discussed in this paper). Viscous terms have been neglected in the above equations, which is a good approximation in atmospheric flow applications due to the Reynolds number being very large (Wyngaard2010). The Coriolis force is also neglected; hence, no veer is present in the simulations. Definitions of the fP correction (which was in fact set to fP=1 in the work of van der Laan et al.2017) and buoyant production, , are deferred to the Sect. 3. The parameters used in the above equations are summarized in Table 1. Finally, Sk and Cε,3 differ slightly from those printed in van der Laan et al. (2017), with the only difference being that here we choose ΦhΦmσθ1; this modeler's choice for turbulent Prandtl number (σθ) avoids the inconsistency mentioned in that paper and makes the model independent of the temperature similarity function Φh.

3 Modification of the kεfP model in the unstable ASL

The background eddy viscosity shown in Fig. 1 is perturbed in the turbine presence and especially so in the wake region; see Fig. 4 for an example with neutral inflow. As explained in Sect. 2.1, νt is very important for the wake development, and the fP correction effectively attenuates the νt perturbation in the interface between the wake and freestream, known as the shear layer, and in the region around the AD to improve wake predictions. This attenuation can also be viewed as a modification of the turbulence scales, ut and t.

Figure 4Streamwise velocity (a, b) and kinematic eddy viscosity (c, d) are normalized by their freestream values, i.e., Uref and νt,ref, respectively. The neutral Iref=6 % case from Fig. 1 with a single NREL5MW turbine is shown.


Figure 5TKE budgets of the 2017 model (a–c) and cstB model (d–f). The profiles are extracted at y/D=0, i.e., in the center of the wake. ζref=-0.5 and Iref=12 % for both rows.


The cause of the νt perturbation in the first place is the large velocity gradients across the AD and the shear layer, which enhances TKE shear production, but other terms of Eq. (10) are also highly active in these regions, and it is this complex interplay together with the fP formulation that in the end determines the wake recovery. The effect of the buoyancy term in this interplay is discussed first and then afterwards the role of fP in the unstable ASL.

3.1 Buoyant production term

The buoyant production of TKE is Bgθ0wθ, and the heat flux is typically obtained using a temperature equation and a flux–gradient relationship. In this work, we pursue an alternative way and investigate two simple parametrizations:


The “2017 model” is the one utilized by van der Laan et al. (2017), van der Laan et al. (2020), Doubrawa et al. (2020) and van der Laan et al. (2021). This model does not require a temperature equation for closure but instead utilizes the temperature similarity function, Φh, and Prandtl number, σθ. The “cstB” model is as the name suggests simply a constant source term throughout the domain, and again no temperature equation is necessary. The k and ε equations are the same for the two models, except some minor changes to Sk, Eq. (16), and Cε3, Eq. (18), are needed for the 2017 model; see Sect. 2.3. To isolate the effect of the  parametrizations, they are tested here with fP=1. A NREL5MW turbine (Jonkman et al.2009) with Uref=8 m s−1 is used for all plots in this section.

Figure 6Disk-averaged streamwise velocity, UAD, for low/high TI and neutral/unstable. ζref=-0.5 for unstable.


The upstream (x/D=-10) budget in Fig. 5 shows the inconsistency of the 2017 model mentioned by van der Laan et al. (2017): the buoyant production goes to B=-u*3κLΦm, although B-u*3κL is expected in the freestream, which can be derived from the ASL definition Bgθ0θws and the Obukhov length definition. The cstB model on the other hand by definition complies with the freestream ASL limit of . Additionally, it can be noted that the cstB model has B/P1 at zref because ζref=-0.5 was used (see Fig. 5.23 of Stull1988).

A clear distinction between the two parametrizations is seen both in the near-wake (x/D=1) and far-wake (x/D=5) TKE budgets: in the top shear layer of the 2017 model simulation, large buoyant production is produced by the large velocity gradients in this region. This is observed in neither direct RANS simulations by El-Askary et al. (2017) nor in wind tunnel experiments by Zhang et al. (2013) or Hancock and Zhang (2015), so this may be deemed as an unphysical artifact. Indeed the 2017 model is derived for a homogeneous ASL, and applying it to a wind turbine wake violates this assumption. The cstB model on the other hand effectively uncouples the buoyant production and wake dynamics. This assumption can partly be justified by the aforementioned studies, which show that temperature changes very little in the wake from the ambient conditions and that the heat flux actually decreases in the wake.

Another deficiency of the unstable 2017 model is illustrated in Fig. 6: for a given Iref, it unphysically predicts slower wake recovery than in neutral as also noted by van der Laan et al. (2021). This is remedied in the cstB model, where a slightly faster wake recovery is seen. It can be noted that in the near and far wake of the cstB model both  and Sk terms are small compared to the other TKE terms (see Fig. 5), and as such it effectively resembles the neutral model, but with the one difference that it has a larger turbulent length scale (see Fig. 2), which explains the faster wake recovery seen in Fig. 6. The  and Sk terms must nevertheless still be retained to enforce the freestream balance of the k and ε equations throughout the domain.

The faster wake recovery of the unstable cstB model compared to the neutral model (as demonstrated in Fig. 6) is also seen in the top row of Fig. 7. The second row shows the shear parameter, σ, which is a non-dimensional scalar metric describing the normalized velocity gradients; see Eq. (23). Both the magnitude and contours of σ are similar for the considered cases; i.e., it is large in the region around the AD and in the shear layers. The freestream shear parameter, σ̃, is however not equal for the two cases (see inlets in third row), which means that σ/σ̃ differs a lot between the two cases. This observation is important because σ/σ̃ is the main parameter used in the fP correction to be discussed in the next section.

Figure 7Normalized streamwise velocity (a, b), shear parameter (c, d) and normalized shear parameter (e, f). The cstB model is used for the unstable case (b, d, f), which has ζref=-0.5.


3.2 Turbulence closure with fP in non-neutral conditions

As stated in the introduction, kε models tend to predict faster wake recovery compared to experiments and LES. This can be corrected by using fP≠1 in the νt definition, Eq. (15), which clearly affects the velocity deficit as shown in Fig. 4. The form of fP used for wakes in the neutral ASL by van der Laan (2014) can be summarized as


As shown in Fig. 7, σ/σ̃ is large in the region surrounding the rotor and in the shear layer (see Fig. 7), which leads to fP<1 and hence a drop of νt in these regions. This is a desirable feature because it corrects the over-diffusiveness of the standard kε model.

It was recognized by van der Laan et al. (2020) that the freestream shear parameter, σ̃, has to be adjusted for MOST inflow in order to have fP=1 in the freestream (the kεfP model should reduce to the standard kε model in the freestream):

(25) σ ̃ = 1 C μ Φ m Φ ε ( general ) .

This can simply be derived by inserting the freestream profiles of U, k and ε (see Sect. 2.1) into the shear parameter definition, Eq. (23); the form of Eq. (25) has been used in all previous papers utilizing the 2017 model for wake modeling.

A more subtle modification arises recognizing that the f0 parameter is also stability-dependent, i.e.,

(26) f 0 = 1 + C μ σ ̃ 2 C R - 1 ( modification 1 ) .

This is actually the form suggested by Apsley and Leschziner (1998), but they were not considering stability effects, i.e., no variation of σ̃ or CR with stability; it has not been used in previous applications of the 2017 model. We note Eq. (26) is consistent with the neutral limit, since σ̃2Cμ-1 for ζ→0. The Rotta constant was calibrated to CR=4.5 for wind turbine wakes in the neutral ASL in the work of van der Laan (2014), and we therefore require CR→4.5 in the neutral limit (ζ→0). One form that satisfies this is

(27) C R = 4.5 + C B B ̃ ε ̃ ( modification 2 ) ,

where B̃/ε̃ is the freestream buoyant production to dissipation ratio and CB is a new parameter to be calibrated. The effect of the stability-dependent σ̃, modification 1 and modification 2 is shown in Fig. 8. When plotted with σ/σ̃ on the abscissa, the stability-dependent σ̃, Eq. (25), has no effect compared to the neutral model (compare “Neutral” and “Unstable” in Fig. 8); this is clear from the definition of fP in Eq. (21). The two modifications increase fP when σ/σ̃>1, i.e., in the wake. This is necessary to compensate for the larger σ/σ̃ encountered with non-neutral inflow; see Fig. 7.

Figure 8Analytical form of the fP correction as function of normalized shear parameter (fP is independent of TI but depends on ζ; the unstable curves use ζref=-0.5 in this plot). The effect of stability-dependent σ̃ (Eq. 25), modification 1 (Eq. 26) and modification 2 (Eq. 27 with CB=10.0) is shown.


When both modifications are used, faster wake recovery for a given Iref occurs in unstable conditions, as shown in Fig. 9, i.e., similar behavior as when fP was “turned off” (i.e., fP=1); see Fig. 6.

Figure 9Same as Fig. 6, but with fP included. CB={0.0,5.0,10.0} are tested for the unstable cstB model.


Table 2Overview of test cases. SWiFT: Doubrawa et al. (2020). NTK41: Machefaux et al. (2016). V80-Abkar: Abkar and Porté-Agel (2015). V80-Keck: Keck et al. (2014). NREL5MW: Churchfield et al. (2012). The air density used for all cases is ρ=1.225 kg m−3.

Download Print Version | Download XLSX

4 Validation with experiments and LES

The cstB model with the fP modifications described in the previous section is tested for the cases summarized in Table 2. Each case is simulated with a range of CB parameters, CB={0.0,5.0,10.0}, while keeping the first term in Eq. (27) fixed to 4.5. The latter has been calibrated for a suite of neutral EllipSys3D LESs by van der Laan (2014), but in practice if it was calibrated with another LES code, a different, optimal value might have been obtained. In the same way, it cannot be expected that a universally valid CB exists, when we compare with results from a range of different codes and experiments. Hence, no optimal CB will be obtained in this section, but rather the qualitative effect of CB is shown.

The numerical setup for each case follows that described in Sect. 2.3; e.g., the cell size and extent of the wake region are scaled with the rotor diameter.

4.1 SWiFT case

A large wake benchmark study was conducted by Doubrawa et al. (2020) to compare various simulation methodologies and codes against lidar measurements in different atmospheric conditions. The measurements were carried out for a Vestas V27 turbine at the Scaled Wind Farm Technology Facility (SWiFT) in Lubbock, Texas, USA, which is an area of flat terrain.

The inflow parameters of the SWiFT row in Table 2 were obtained from the ensemble average of five 10 min averages from a meteorological mast located 2.5 D upstream of the turbine. Note that the stability parameter was measured to ζ=-0.089 at z=10 m, which at hub height corresponds to ζref=-0.29. Also, the streamwise turbulence intensity was measured at hub height to Iu,ref=12.6 %, which is converted to the total turbulence intensity as Iref0.8Iu,ref (van der Laan et al.2015b). This conversion could actually be slightly different in the unstable ASL because the ratios of velocity variance change with stability (e.g., Chougule et al.2018), but unfortunately only the vertical velocity variance follows MOST, so that no general surface layer formula can be constructed (Panofsky and Dutton1984; Wyngaard2010). The operational state parameters {CT, P and Ω} were taken from the OpenFAST steady-state curves, which were supplied for the benchmark.

For the unstable SWiFT case, the wake profile was only measured at 3 D downstream, and the results are shown in Fig. 10. Three different LES codes were used in the benchmark, and the filled area in Fig. 10 represents the spread of the LES results. It can be seen that all three LESs underpredict velocity deficit compared to the lidar measurements, which highlights the fundamental problem of comparing measurements with numerical models: even highly computational expensive simulations do not always match experimental results. This must be either due to experimental errors in the provided input data or because the idealizations used for the LESs are too simple to capture the wake behavior.

Figure 10The unstable SWiFT case. The fixed frame of reference experimental and LES results were digitized from Doubrawa et al. (2020).


Figure 11The unstable NTK41 case, where LES and experimental results were digitized from Machefaux et al. (2016). Profiles extracted at z=zref.


RANS can generally not be expected to perform better than a well-performed LES, and if it does, it is likely due to fortunate error cancelations. Therefore, from a theoretical point of view one could argue that the performance of RANS should mainly be assessed with regards to how well it matches the LES results. Both RANS and LES use many of the same idealizations (uniform roughness, flat terrain, homogeneous inflow, etc.), and indeed our RANS results in Fig. 10 are also closer to the LES results than to the experimental results.

For the SWiFT case, the 2017 model seemingly performs better than the cstB model, but this is probably due to fortuitous model error and/or unaccounted mesoscale effects. Model error was expected because the neutral EllipSys3D RANS simulation in the SWiFT study (surprisingly) did not compare well with experimental results (Doubrawa et al.2020), despite EllipSys3D having been validated by van der Laan (2014) for numerous neutral cases. The very low wind veer of the stable SWiFT case (see Doubrawa et al.2020) indicates that mesoscale effects were present during the experimental campaign.

4.2 NTK41 case

A Nordtank NTK41 500 kW wind turbine was installed at what is now the Risø campus of the Technical University of Denmark in 1992 and was used for many research studies before its decommissioning in 2021. Among these studies, the NTK41 test case of this paper is based on lidar measurements and LESs conducted by Machefaux et al. (2016). They used two different models for their LESs; the results included in Fig. 11 (along with our RANS results) are from their more advanced model, which they called the “LES-ABL” or “extended model”. The inflow parameters (“NTK41” row in Table 2) and lidar measurements were ensemble-averaged over 20 10 min averages.

Figure 12The unstable V80-Abkar case, where LES results were digitized from Abkar and Porté-Agel (2015). Profiles extracted at z=zref. Note that total TI (d–f) is based on Uref and not the local velocity.


The NTK41 turbine is a stall-regulated wind turbine and is therefore operated at constant rotational speed independent of the inflow wind speed (Hansen2015), in this case at Ω=27.1 rpm. The thrust coefficient for the unstable case of Machefaux et al. (2016) was measured with strain gauges to be CT,meas=0.71, while their LES gave CT,LES=0.83. Looking up the thrust curve of the NTK41 turbine at Uref=6.8 m s−1 also gives CT,curve=0.83, so this will be used in the present RANS simulations. They argued that the lower thrust coefficient of their measurement could be explained by the large uncertainty of the strain gauges. Finally, the measured power was Pmeas=120 kW, while PLES=127 kW and Pcurve=125 kW, where the latter will be used to set CP for the AD model of the RANS simulations.

Figure 11 shows that the cstB model matches the LES and experimental data better than the 2017 model, although a still faster wake recovery is seen for both of these reference data. Compared to more conventional LES setups (e.g., V80-Abkar, V80-Keck and NREL5MW cases), the LES model used by Machefaux et al. (2016) is simplified by using a modified Mann box for its inflow and a slip condition at the bottom wall; together these add some uncertainty to the LES results. The experimental wake data can also be expected to have large uncertainties and/or biases connected to them, but the sources and sizes of these were not discussed in detail by Machefaux et al. (2016).

4.3 V80-Abkar case

Abkar and Porté-Agel (2014) investigated the effect of atmospheric stability using LES and used the results to modify the analytical Bastankhah wake model (Abkar and Porté-Agel2015). They studied the wake of a single Vestas V80 turbine (known from, e.g., the Horns Rev 1, North Hoyle and Princess Amalia wind farms), which has been used in many previous wake studies. The turbine was modeled in their studies by an airfoil AD and operated at Ω=16.1 rpm. Neither CT nor P was mentioned in the two papers, so in the following RANS simulations the values deduced from the power and thrust coefficient curves of Hansen et al. (2012) evaluated at Uref=8 m s−1 were used: CT=0.81 and P=696 kW.

Relative to the LES, the velocity deficits shown in Fig. 12 are overpredicted by the 2017 model for all three downstream distances, while the cstB model corrects this especially well in the far wake. Besides the velocity deficit, the TI based on the freestream velocity was also available for this case and is plotted in the lower row of Fig. 12 with the RANS results. The wake TI is overpredicted by RANS, which is also typically seen for neutral RANS simulations (van der Laan2014). The slower wake recovery of the 2017 model is consistent with less turbulence development in the near-AD region, which might explain its seemingly better TI prediction; however, we prioritize predicting a correct velocity deficit, since this is used for AEP calculations.

Figure 13The unstable V80-Keck case, where LES results were digitized from Keck et al. (2014). Profiles extracted at z=zref. Note that streamwise TI (d–f) is based on Uref and σu23k/0.8 for the RANS simulations.


4.4 V80-Keck case

This case is based on a LES from Keck et al. (2014), where the SOWFA solver was used with a similar setup as in the work of Churchfield et al. (2012), i.e., using a precursor simulation for the inflow and modeling the turbine with actuator lines (ALs). More specifically the unstable North Hoyle row A case is considered here; it features four V80 turbines spaced 11 D apart, using the inflow parameters described in the “V80-Keck” row of Table 2. Wake data are available at x/D={4,5,6} downstream of the first turbine, and the induction effect of the downstream turbines on this data should therefore be minimal; hence, they are omitted from the RANS simulations. The inflow wind speed is Uref=8 m s−1, the same as in the V80-Abkar case; thus, the same operational state of the wind turbine is utilized (see Table 2). The streamwise TI given by Keck et al. (2014) is converted to total TI with Iref0.8Iu,ref, similar to the method used in the SWiFT case.

Figure 13 shows that the cstB model improves the velocity deficit prediction over the old 2017 model when comparing with the LES results, which were digitized from Fig. 11 in Keck et al. (2014) (note that a typo is present in that figure; i.e., label should be “[D]” instead of “[R]”). Streamwise TI was also available at the same downstream distances, and in the RANS simulations it was obtained by converting from the total TI as described above. In both the streamwise TI and velocity deficit LES data, a misalignment of the wake center can be observed, which Keck et al. (2014) explains with the fact that only 10 min of LES data was averaged. This is especially visible in the streamwise TI plots, but nevertheless it seems that the cstB model predicts streamwise TI in the right range except for an overprediction in the wake center, which was also seen in the previous case.

4.5 NREL5MW case

The last case is based on the LES studies by Churchfield et al. (2012), more specifically their “U-L case” (see inflow parameters in the “NREL5MW” row of Table 2). They model two NREL5MW turbines with the actuator line (AL) methodology coupled to the aeroelastic FAST solver, and the turbines are spaced 7 D apart. For the RANS simulations of the present study, we omit the second turbine and only compare with the first wake of the LES study. To avoid biases from the induction zone of the second turbine, we only consider wake results ≥2D upstream of the second turbine.

The steady-state power, thrust coefficient and rotational speed were not given in the paper, so therefore the steady-state curves from the DTU in-house aeroelastic solver, HAWCStab2, were used. These are similar to the curves shown by Jonkman et al. (2009), except that the thrust of Jonkman et al. (2009) also includes gravity and therefore cannot be used to obtain the aerodynamic thrust.

The velocity deficit of the new cstB model compares well with the LES data, especially so for CB=5.0; see Fig. 14.

Figure 14The unstable NREL5MW case, where LES results were digitized from Churchfield et al. (2012). Profiles extracted at z=zref.


5 Conclusions

We have proposed a simple kε RANS model, the cstB model, for simulation of wind turbine wakes in the unstable surface layer. The model does not require an additional temperature equation and instead bases the TKE buoyancy production on MOST and the assumption that it is decoupled from the wake dynamics, which means that the buoyant production of TKE is constant throughout the domain, even in the wake region, hence the name cstB. Wind tunnel studies and simulations have hinted that the latter assumption is reasonable, but a more thorough investigation would be beneficial for developing simple, non-neutral wake models. For example comparisons with more detailed LES data (temperature profiles, shear production, buoyancy production, Reynolds stress tensor, anisotropy tensor, derived kinematic viscosity, etc.) would be useful.

Originally developed to account for the general over-diffusive nature of kε models in wind turbine wakes under neutral conditions, here the fP correction is combined with the new cstB model by making two non-neutral modifications. These introduce a new parameter, CB; it is a free parameter analogous to CR in the original fP formulation. Both modifications are consistent in the sense that the new, non-neutral fP formulation becomes equal to the original neutral fP form for ζ→0. By using this updated fP model with the cstB model, a faster wake recovery is obtained for unstable conditions over neutral conditions, when TI is fixed, as was also the case when no fP model was applied.

The cstB model with the modified fP function was generally found to perform better than the previous model of van der Laan et al. (2017) with the old fP formulation, in terms of velocity deficit profiles from five different reference cases found in studies from the literature. Based on these comparisons, we recommend CB=5.0 to be used but also acknowledge that each reference case was originally conducted with different numerical and experimental setups and that further studies are needed to conclude on a more certain CB value, which could also be slightly code-dependent, as has been seen for CR in the original fP model.

Testing the cstB model behavior in more complicated scenarios, i.e., aligned row cases (see Appendix A3), full wind farms, complex terrain, AEP calculations, etc., is a natural next step to map the applicability and limitations of the model. Also the extension to stable conditions for the cstB model is straightforward, but as the surface layer height in such conditions is small compared to modern turbine dimensions, it is questionable if the cstB model will be viable or if full ABL models are more appropriate. Yet another question to be answered is the effect of freestream turbulence anisotropy, which changes with stability and cannot be modeled with turbulence models based on the standard Boussinesq hypothesis.

Appendix A: Simulation details

A1 Grid study

Earlier studies by van der Laan et al. (2015b) have shown that a domain resolution of eight cells per diameter is sufficient to obtain grid independence for wakes in the neutral ASL. A range of domain and AD resolutions are here tested for the new cstB model and fP modifications with the TI=12 % and CB=5.0 case also used in Fig. 9. The domain size is described in Fig. 3, and the Joukowsky AD (see next section) is used. Disk averages of velocity and TI are evaluated 1 D downstream of the turbine to verify grid independence in Fig. A1.

Figure A1Grid-independence study. “Ref” means reference, i.e., the finest resolution available. Metrics are evaluated at x/D=1. The cores used for the increasing domain resolutions are {45, 54, 63, 54, 60, 57}, respectively (non-constant because the domains are decomposed in different number of blocks).


Based on this small grid study, a domain resolution of 10 cells per diameter and an AD resolution of (Nφ,Nr)=(32,32) are chosen for the current investigation. The difference between this resolution and the reference resolution is less than 2 % for the velocity metric and less than 6 % for the TI metric (the differences decrease with downstream position; e.g., at x/D=5 the differences are only approximately 0.2 % for the velocity metric and 2 % for the TI metric). A simulation with this choice can be executed in about 1 wall clock minute on 63 cores (AMD EPYC 7351 processors are used).

A2 The Joukowsky AD method

In summary, the surface force distributions (unit: N m−2) on the AD are calculated in each iteration as


Figure A2The normal blade loading, Fn [N m−1], and tangential blade loading, Ft [N m−1], are normalized by the density, ρ, rotor radius, R, and freestream hub height velocity, Uref. The blade loadings for the Joukowsky AD and airfoil AD have been obtained by azimuthal averages, while Fn=2πrT/ANb is prescribed a priori for the uniform AD.


Figure A3Disk average of velocity, UAD, and disk average of turbulence intensity, IAD.


Figure A4Aligned row of 10 turbines with inflow similar to the Case 5 of van der Laan et al. (2021).


Here, fn,ij and fθ,ij are the normal and azimuthal surface force distributions at the (i,j)th AD element (i: radial direction, j: azimuthal direction), which are applied in the CFD domain using the methodology described by Réthoré et al. (2014) and Troldborg et al. (2015). Fn0ijfn,ij0Aij is the total normal force of the unscaled distribution, P0Urefλiχijfθ,ij0Aij is the total power of the unscaled distribution, λΩR/Uref is the tip-speed ratio, A is the area of the AD, Aij is the area of the (i,j)th AD element, χiri/R is the local normalized radius, Uij is the normal velocity at the (i,j)th AD element, Nb=3 is the number of blades, q0 is the normalized circulation, g is the root correction and F is the tip correction. The latter two are obtained with Delery's root correction (with parameters a=2.335, b=4.0 and δ=0.25) and Prandtl's tip correction, respectively; see Eqs. (A4) and (A5).

The normal and tangential loadings for the same case as used in Sect. A1 are compared between the uniform AD, airfoil AD and Joukowsky AD, the latter with two different δ's, in Fig. A2. Clearly, the Joukowsky AD with δ=0.25 produces similar loadings as the airfoil AD, e.g., qualitatively correct root behavior, tip behavior and constant tangential loading region, which makes it superior to the uniform AD.

Finally, in Fig. A3 the velocity and TI disk averages follow the same trend for all four AD methods, but with a slightly larger velocity deficit for the airfoil AD, possibly because of its also slightly larger blade loadings; see Fig. A2. The thrust coefficient of the uniform AD and Joukowsky AD is CT=0.77, which from 1D momentum theory should give U/Uref=1-0.51-1-CT0.74 at the rotor plane. This is not exactly observed in Fig. A3, but contrary to ideal 1D momentum theory our CFD simulation also includes atmospheric turbulence and shear.

A3 Aligned row case

Even though the focus of the present paper is on single wakes, wake simulations will in practice most often involve wake interaction. Contrary to engineering models, there is no need for empirical wake superposition methods in RANS, since the wake interaction automatically results from solving the RANS equations.

The Case 5 from the study of van der Laan et al. (2021), which consists of an aligned row of 10 NREL5MW turbines spaced 7 D apart, is simulated with a similar setup, and the normalized, disk-averaged velocity recovery is shown in Fig. A4. For reference, a neutral kεfP simulation with the same inflow speed and TI is simulated. The 2017 model performs poorly in the sense that it actually recovers slower than the neutral reference and that no equilibrium wake–wake interaction seems to occur. On the other hand, the cstB model predicts faster wake recovery for the first three to four turbines, while it goes to the same wake deficit as the neutral model in the fully developed part of the wind farm. More studies and validations with LES/experiments for similar wake interaction cases should be undertaken in the future.

Code availability

EllipSys3D and PyWakeEllipSys are proprietary software of DTU. Information about the latter can, however, be freely accessed at (van der Laan2022).

Data availability

The RANS results were generated with DTU's proprietary software, but the data presented can be made available by contacting the corresponding author. Interested parties are also welcome to hand-digitize the results and use them as reference in other publications.

Author contributions

MB performed the RANS simulations and proposed the modifications to the turbulence model. All authors (MB, MPvdL and MK) contributed to derivation of the new model and article writing.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors are grateful to the creators of the LES and lidar results used in the validation section. We would also like to thank the two reviewers for their feedback and suggestions.

Review statement

This paper was edited by Sandrine Aubrun and reviewed by two anonymous referees.


Abkar, M. and Porté-Agel, F.: The effect of atmospheric stability on wind-turbine wakes: A large-eddy simulation study, J. Phys.: Conf. Ser., 524, 1–9,, 2014. a

Abkar, M. and Porté-Agel, F.: Influence of atmospheric stability on wind-turbine wakes: A large-eddy simulation study, Phys. Fluids, 27, 035104,, 2015. a, b, c, d

Alinot, C. and Masson, C.: Aerodynamic simulations of wind turbines operating in atmospheric boundary layer with various thermal stratifications, ASME 2002 Wind Energy Symposium, WIND2002, 206–215,, 2002. a, b

Apsley, D. D. and Castro, I. P.: A limited-length-scale kϵ model for the neutral and stably-stratified atmospheric boundary layer, Bound.-Lay. Meteorol., 83, 75–98,, 1997. a

Apsley, D. D. and Leschziner, M. A.: A new low-Reynolds-number nonlinear two-equation turbulence model for complex flows, Int. J. Heat Fluid Flow, 19, 209–222,, 1998. a

Chougule, A., Mann, J., Kelly, M., and Larsen, G. C.: Simplification and Validation of a Spectral-Tensor Model for Turbulence Including Atmospheric Stability, Bound.-Lay. Meteorol., 167, 371–397,, 2018. a

Churchfield, M. J., Lee, S., Michalakes, J., and Moriarty, P. J.: A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics, J. Turbulence, 13, 1–32,, 2012. a, b, c, d, e

Crespo, A., Manuel, F., Moreno, D., Fraga, E., and Hernandez, J.: Numerical analysis of wind turbine wakes, in: Proceedings of the Delphi Workshop on Wind Energy Applications, 15–25, (last access: 30 March 2022), 1985. a

Doubrawa, P., Quon, E. W., Martinez‐Tossas, L. A., Shaler, K., Debnath, M., Hamilton, N., Herges, T. G., Maniaci, D., Kelley, C. L., Hsieh, A. S., Blaylock, M. L., Laan, P., Andersen, S. J., Krueger, S., Cathelain, M., Schlez, W., Jonkman, J., Branlard, E., Steinfeld, G., Schmidt, S., Blondel, F., Lukassen, L. J., and Moriarty, P.: Multimodel validation of single wakes in neutral and stratified atmospheric conditions, Wind Energy, 23, we.2543,, 2020. a, b, c, d, e, f

Dyer, A. J.: A review of flux-profile relationships, Bound.-Lay. Meteorol., 7, 363–372, 1974. a

El-Askary, W. A., Sakr, I. M., AbdelSalam, A. M., and Abuhegazy, M. R.: Modeling of wind turbine wakes under thermally-stratified atmospheric boundary layer, J. Wind Eng. Indust. Aerodynam., 160, 1–15,, 2017. a, b

Foken, T.: 50 years of the Monin-Obukhov similarity theory, Bound.-Lay. Meteorol., 119, 431–447,, 2006. a

Ghaisas, N. S., Archer, C. L., Xie, S., Wu, S., and Maguire, E.: Evaluation of layout and atmospheric stability effects in wind farms using large-eddy simulation, Wind Energy, 20, 1227–1240,, 2017. a

Gryning, S. E., Batchvarova, E., Brümmer, B., Jørgensen, H., and Larsen, S.: On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer, Bound.-Lay. Meteorol., 124, 251–268,, 2007. a

Han, X. X., Liu, D. Y., Xu, C., Shen, W. Z., Li, L. M., and Xue, F. F.: A BEM-based actuator disk model for wind turbine wakes considering atmospheric stability, Preprint, 1–24,, 2019. a

Hancock, P. E. and Zhang, S.: A Wind-Tunnel Simulation of the Wake of a Large Wind Turbine in a Weakly Unstable Boundary Layer, Bound.-Lay. Meteorol., 156, 395–413,, 2015. a

Hansen, K. S., Barthelmie, R. J., Jensen, L. E., and Sommer, A.: The impact of turbulence intensity and atmospheric stability on power deficits due to wind turbine wakes at Horns Rev wind farm, Wind Energy, 15, 183–196,, 2012. a, b

Hansen, M. O. L.: Aerodynamics of Wind Turbines, 3rd Edn., Routledge,, 2015. a

Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5-MW Reference Wind Turbine for Offshore System Development, NREL technical report, NREL,, 2009. a, b, c

Kaimal, J. and Finnigan, J. J.: Atmospheric Boundary Layer Flows, Oxford University Press,, 1994. a

Kasmi, A. E. and Masson, C.: An extended k2 model for turbulent flow through horizontal-axis wind turbines, J. Wind Eng. Indust. Aerodynam., 96, 103–122,, 2008. a

Keck, R. E., de Maré, M., Churchfield, M. J., Lee, S., Larsen, G., and Aagaard Madsen, H.: On atmospheric stability in the dynamic wake meandering model, Wind Energy, 17, 1689–1710,, 2014. a, b, c, d, e, f, g

Kelly, M. and Gryning, S. E.: Long-Term Mean Wind Profiles Based on Similarity Theory, Bound.-Lay. Meteorol., 136, 377–390,, 2010. a, b

Li, N., Liu, Y., Li, L., Chang, S., Han, S., Zhao, H., and Meng, H.: Numerical simulation of wind turbine wake based on extended kϵ turbulence model coupling with actuator disc considering nacelle and tower, IET Renewable Power Generation,, 2020. a

Machefaux, E., Larsen, G. C., Koblitz, T., Troldborg, N., Kelly, M. C., Chougule, A., Hansen, K. S., and Rodrigo, J. S.: An experimental and numerical study of the atmospheric stability impact on wind turbine wakes, Wind Energy, 19, 1785–1805,, 2016. a, b, c, d, e, f, g

Magnusson, M. and Smedman, A. S.: Influence of Atmospheric Stability on Wind Turbine Wakes, Wind Eng., 18, 139–152, 1994. a

Maronga, B. and Reuder, J.: On the formulation and universality of Monin–Obukhov similarity functions for mean gradients and standard deviations in the unstable surface layer: Results from surface-layer-resolving large-eddy simulations, J. Atmos. Sci., 74, 989–1010,, 2017. a

Michelsen, J. A.: Basic3D : A platform for development of multiblock PDE solvers, Technical Report AFM 92-05, Technical University of Denmark, Lyngby, 1992. a, b

Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Tr. Akad. Nauk SSSR Geophiz. Inst., 24, 163–187, 1954. a

Okulov, V. L. and Sørensen, J. N.: Maximum efficiency of wind turbine rotors using Joukowsky and Betz approaches, J. Fluid Mech., 649, 497–508,, 2010. a

Panofsky, H. and Dutton, J.: Atmospheric Turbulence, John Wiley & Sons, Ltd, ISBN 0471057142, 1984. a

Pope, S. B.: Turbulent Flows, Cambridge University Press,, 2000. a

Porté-Agel, F., Wu, Y. T., Lu, H., and Conzemius, R. J.: Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms, J. Wind Eng. Indust. Aerodynam., 99, 154–168,, 2011. a, b

Porté-Agel, F., Bastankhah, M., and Shamsoddin, S.: Wind-Turbine and Wind-Farm Flows: A Review, vol. 174, Springer Netherlands,, 2020. a

Rados, K. G., Prospathopoulos, J. M., Stefanatos, N. C., Politis, E. S., Chaviaropoulos, P. K., and Zervos, A.: CFD modeling issues of wind turbine wakes under stable atmospheric conditions, in: European Wind Energy Conference and Exhibition 2009, EWEC 2009, 6, Marseille, 4085–4092, ISBN 9781615677467,, (last access: 30 March 2022), 2009. a, b

Réthoré, P.-E.: Wind Turbine Wake in Atmospheric Turbulence, PhD thesis, Aalborg University, Aalborg, (last access: 30 March 2022), 2009. a

Réthoré, P.-E., van der Laan, P., Troldborg, N., Zahle, F., and Sørensen, N. N.: Verification and validation of an actuator disc model, Wind Energy, 17, 919–937,, 2014. a

Sathe, A., Mann, J., Barlas, T., Bierbooms, W., and van Bussel, G.: Influence of atmospheric stability on wind turbine loads, Wind Energy, 16, 1013–1032,, 2013. a

Sørensen, J. N. and Andersen, S. J.: Validation of analytical body force model for actuator disc computations, J. Phys.: Conf. Ser., 1618, 1–8,, 2020. a

Sørensen, J. N. and Kock, C. W.: A model for unsteady rotor aerodynamics, J. Wind Eng. Indust. Aerodynam., 58, 259–275,, 1995. a

Sørensen, J. N., Nilsson, K., Ivanell, S., Asmuth, H., and Mikkelsen, R. F.: Analytical body forces in numerical actuator disc model of wind turbines, Renew. Energy, 147, 2259–2271,, 2020. a, b

Sørensen, N. N.: General purpose flow solver applied to flow over hills, PhD thesis, Technical university of Denmark, Risø National Laboratory, (last access: 30 March 2022), 1995. a, b, c

Sørensen, N. N., Bechmann, A., Johansen, J., Myllerup, L., Botha, P., Vinther, S., and Nielsen, B. S.: Identification of severe wind conditions using a Reynolds Averaged Navier-Stokes solver, J. Phys.: Conf. Ser., 75, 012053,, 2007. a

Stull, R. B.: An introduction to boundary layer meteorology, An introduction to boundary layer meteorology, Kluwer Academic Publishers,, 1988. a

Troldborg, N., Sørensen, N. N., Réthoré, P.-E., and van der Laan, M. P.: A consistent method for finite volume discretization of body forces on collocated grids applied to flow through an actuator disk, Comput. Fluids, 119, 197–203,, 2015. a, b

van der Laan, M. P.: Efficient Turbulence Modeling for CFD Wake Simulations, PhD thesis, DTU Wind Energy, Technical University of Denmark, (last access: 30 March 2022), 2014. a, b, c, d, e, f

van der Laan, M. P.: PyWakeEllipSys documentation,, last access: 30 March 2022. a

van der Laan, M. P., Sørensen, N. N., Réthoré, P.-E., Mann, J., Kelly, M. C., and Troldborg, N.: The k-ϵ-fp model applied to double wind turbine wakes using different actuator disk force methods, Wind Energy, 18, 2223–2240,, 2015a. a, b

van der Laan, M. P., Sørensen, N. N., Réthoré, P. E., Mann, J., Kelly, M. C., Troldborg, N., Schepers, J. G., and Machefaux, E.: An improved k-ϵ model applied to a wind turbine wake in atmospheric turbulence, Wind Energy, 18, 889–907,, 2015b. a, b, c

van der Laan, P., Sørensen, N. N., Réthoré, P.-E., Mann, J., Kelly, M. C., Troldborg, N., Hansen, K. S., and Murcia, J. P.: The k-ϵ-fp model applied to wind farms, Wind Energy, 18, 2065–2084,, 2015c.  a

van der Laan, M. P., Kelly, M. C., and Sørensen, N. N.: A new k-ϵ model consistent with Monin-Obukhov similarity theory, Wind Energy, 20, 479–489,, 2017. a, b, c, d, e, f, g, h, i, j

van der Laan, M. P., Andersen, S. J., Kelly, M., and Baungaard, M. C.: Fluid scaling laws of idealized wind farm simulations, J. Phys.: Conf. Ser., 1618, 1–10,, 2020. a, b, c

van der Laan, M. P., Baungaard, M., and Kelly, M.: Inflow modeling for wind farm flows in RANS, J. Phys.: Conf. Seri., 1934, 012012,, 2021. a, b, c, d, e, f

Wyngaard, J. C.: Turbulence in the atmosphere, Cambridge University Press,, 2010. a, b

Xie, S. and Archer, C. L.: A Numerical Study of Wind-Turbine Wakes for Three Atmospheric Stability Conditions, Bound.-Lay. Meteorol., 165, 87–112,, 2017. a

Zhang, W., Markfort, C. D., and Porté-Agel, F.: Wind-Turbine Wakes in a Convective Boundary Layer: A Wind-Tunnel Study, Bound.-Lay. Meteorol., 146, 161–179,, 2013. a

Short summary
Wind turbine wakes are dependent on the atmospheric conditions, and it is therefore important to be able to simulate in various different atmospheric conditions. This paper concerns the specific case of an unstable atmospheric surface layer, which is the lower part of the typical daytime atmospheric boundary layer. A simple flow model is suggested and tested for a range of single-wake scenarios, and it shows promising results for velocity deficit predictions.
Final-revised paper