Articles | Volume 9, issue 5
Research article
08 May 2024
Research article |  | 08 May 2024

The multi-scale coupled model: a new framework capturing wind farm–atmosphere interaction and global blockage effects

Sebastiano Stipa, Arjun Ajay, Dries Allaerts, and Joshua Brinkerhoff

The growth in the number and size of wind energy projects in the last decade has revealed structural limitations in the current approach adopted by the wind industry to assess potential wind farm sites. These limitations are the result of neglecting the mutual interaction of large wind farms and the thermally stratified atmospheric boundary layer. While currently available analytical models are sufficiently accurate to conduct site assessments for isolated rotors or small wind turbine clusters, the wind farm's interaction with the atmosphere cannot be neglected for large-size arrays. Specifically, the wind farm displaces the boundary layer vertically, triggering atmospheric gravity waves that induce large-scale horizontal pressure gradients. These perturbations in pressure alter the velocity field at the turbine locations, ultimately affecting global wind farm power production. The implication of such dynamics can also produce an extended blockage region upstream of the first turbines and a favorable pressure gradient inside the wind farm. In this paper, we present the multi-scale coupled (MSC) model, a novel approach that allows the simultaneous prediction of micro-scale effects occurring at the wind turbine scale, such as individual wake interactions and rotor induction, and meso-scale phenomena occurring at the wind farm scale and larger, such as atmospheric gravity waves. This is achieved by evaluating wake models on a spatially heterogeneous background velocity field obtained from a reduced-order meso-scale model. Verification of the MSC model is performed against two large-eddy simulations (LESs) with similar average inflow velocity profiles and a different capping inversion strength, so that two distinct interfacial gravity wave regimes are produced, i.e. subcritical and supercritical. Interfacial waves can produce high blockage in the first case, as they are allowed to propagate upstream. On the other hand, in the supercritical regime their propagation speed is less than their advection velocity, and upstream blockage is only operated by internal waves. The MSC model not only proves to successfully capture both local induction and global blockage effects in the two analyzed regimes, but also captures the interaction between the wind farm and gravity waves, underestimating wind farm power by about only 2 % compared with the LES results. Conversely, wake models alone cannot distinguish between differences in thermal stratification, even if combined with a local induction model. Specifically, they are affected by a first-row overprediction bias that leads to an overestimation of the wind farm power by 13 % to 20 % for the modeled regimes.

1 Introduction

The growth in the number and size of wind energy projects in the last decade has revealed structural limitations in the current approach adopted by the wind industry to assess potential wind farm sites. According to updates published by Ørsted (2019) on its long-term financial targets, the existing engineering models used for wind farm site assessments are unable to accurately represent two important physical aspects that affect wind farm power production forecasts, namely, wake and blockage effects. Wake effects are the interaction of individual turbine wakes as they merge to form an extended region of reduced momentum downstream of the wind farm. Despite being characterized by complex flow dynamics, such a physical process is frequently described using a simple combination of individual wake deficits. This leads to some structural limitations in capturing additional phenomena such as the effect of increased shear stress on vertical momentum fluxes (Calaf et al.2010), their effect on the wake, and the formation of an internal boundary layer (IBL) that grows from the wind farm start (Meneveau2011). In addition, since more plants are being constructed in the proximity of pre-existing ones (Nygaard et al.2020), cluster wake evolution is of great interest in the context of farm–farm wake interactions. Regarding blockage (Bleeg et al.2018), this is defined as the flow deceleration upstream of the wind farm produced by the combination of individual (local) blockage created by individual turbines with a bulk (global) deceleration caused by the wind farm interacting with the thermally stratified atmosphere.

Many years of research on turbine-level induction resulted in different local-blockage models, such as the empirical formulation of Troldborg and Meyer Forsting (2017), the Rankine half-body of Gribben and Hawkes (2019), the vortex cylinder model (Branlard and Gaunaa2014; Branlard and Meyer Forsting2020), the self-similar model of Branlard et al. (2020), and the linearized model by Segalini (2021). Initially, global blockage effects produced by the entire cluster have been assessed by linearly superimposing local blockage models. As long as thermal stratification of the Earth's atmosphere was not considered, this superposition compared well against numerical simulations (see for example Branlard and Meyer Forsting2020). However, once atmospheric stability is considered, Centurelli et al. (2021) showed that linearly combining individual effects underpredicts the global blockage, and the inaccuracy worsens for larger wind farms. On the other hand, several numerical studies provided evidence that global blockage is a complex physical phenomenon that mainly depends on the structure of the potential temperature profile (Wu and Porté-Agel2017; Allaerts and Meyers2017, 2018). Specifically, upstream wind slowdown is produced by large-scale pressure perturbations that arise as a consequence of interfacial waves in the inversion layer and internal gravity waves aloft. To model these effects, Allaerts and Meyers (2019) developed the three-layer model (3LM), enhancing the earlier two-layer model (2LM) by Smith (2010), which enabled them to capture wind farm–gravity wave interactions in a fast engineering model for the first time. The idea behind the 3LM is to provide a meso-scale velocity correction for the free-stream wind speed used in conventional wake models (Katić et al.1987; Ainslie1988; Larsen1988; OpenWind2010; Bastankhah and Porté-Agel2014; Niayifar and Porté-Agel2016). Such a correction has the effect of reducing the free-stream velocity, reflecting the influence of atmospheric gravity waves in decelerating the wind incident to the wind farm. However, despite the important physical insights at the meso-scale level offered by the 3LM, the coupling between the turbine-scale wake effects and the meso-scale global effects is weak. In fact, the model is unable to adequately transfer the entirety of the large-scale physical processes to the flow at the turbine scale, failing to capture, for example, the beneficial effects sometimes produced inside the wind farm.

Although it has been used to enhance model performance in capturing wake effects rather than blockage, the idea of coupling two models operating at different scales was initially introduced by Frandsen et al. (2006). Later, Stevens et al. (2015) introduced the coupled wake boundary layer (CWBL) model, where the top-down model of Calaf et al. (2010) (later extended by Meneveau2011, and Stevens2016) was coupled with a wake model to enhance wind farm power predictions, assuming that the wind farm internal boundary layer (IBL) is fully developed. Further improvements of the CWBL model were made by Shapiro et al. (2019), and later by Starke et al. (2021), in what is currently referred to as the area-localized coupled (ALC) model, extending its applicability to those regions where the IBL is not fully developed. The coupling idea has also been used by Nishino and Dunstan (2020) to develop a two-scale momentum theory, which couples conventional momentum theory with the large-scale momentum balance above the wind farm layer.

Except for the 2LM and 3LM, all coupled models mentioned above do not aim to capture the meso-scale wind farm blockage effects. In particular, they do not account for thermal stratification above the atmospheric boundary layer (ABL), with the result that meso-scale stability and its effects on blockage remain a major source of uncertainty in wind farm power predictions. To fill this gap, we introduce the multi-scale coupled (MSC) model, a framework that allows one to couple the 3LM with any wake model while also including local-induction effects. The MSC overcomes the above-mentioned deficiency characterizing the original 3LM and its coupling technique, establishing a new modular framework in which additional meso-scale effects can be easily added in the future. In the present study, verification of the MSC model is performed against large-eddy simulations (LESs) of a wind farm consisting of 100 NREL 5 MW turbines in an aligned configuration, under subcritical and supercritical interface wave regimes. These conditions are obtained by adjusting the so-called inversion Froude number Fr=Ub/gH, where Ub is the bulk velocity inside the boundary layer, g is the reduced gravity, and H is the inversion height (Smith2010). In the subcritical regime, characterized by Fr<1, interface waves exhibit a higher group velocity than the advection speed, allowing them to propagate upstream. Conversely, in supercritical conditions (Fr>1) upstream perturbations are mainly due to internal waves, as interface waves can only propagate downstream. To highlight the improved accuracy of the MSC model, the LES results are also compared against the wake-only approach, with and without local induction, and to the original 3LM formulation by Allaerts and Meyers (2019).

The present paper is organized as follows. Section 2 briefly introduces the original 3LM model, followed by a detailed description of the proposed MSC model framework. Then, Sect. 3 presents the LES that will be used for model verification. Section 4 focuses on evaluating the predictive accuracy of the proposed engineering model, from the viewpoint of both the turbine data and the mean velocity field surrounding the wind farm. Finally, Sect. 5 highlights the conclusions of the present study.

2 Methodology

The multi-scale coupled (MSC) model describes the interaction between physical processes that take place inside the atmospheric boundary layer, such as turbine wakes or local rotor induction, and phenomena that exist aloft, characterized by larger spatial scales, such as wind farm–gravity wave interaction. The micro-scale to meso-scale coupling is deemed possible by introducing the notion of a background velocity field, which represents the effect of gravity wave-induced pressure gradients on an initially uniform velocity field. This background wind is then combined with wake and local induction effects, allowing the MSC model to capture not only the latter physical processes, but also additional phenomena such as global blockage effects, the favorable pressure gradient inside the wind farm, and large-scale velocity oscillations in the wind farm wake for those conditions where gravity wave trains induce alternating favorable–unfavorable pressure gradients past the wind farm. The MSC model is formulated as a modular combination of different sub-models, which are suitably integrated and made interdependent. This approach does not require additional tuning parameters or individual sub-model re-tuning.

Although the MSC framework generally differs from the original 3LM developed by Allaerts and Meyers (2019), especially in the solution strategy and sub-model coupling, the two models share some common aspects, such as the chosen wake model or the meso-scale modeling approach. For this reason, we first introduce the original 3LM formulation in Sect. 2.1, together with all common aspects between the two models. Later, in Sect. 2.2, we present the new methodology, focusing on its differences and improvements with respect to the original formulation.

2.1 The three-layer model

Developed by Allaerts and Meyers (2019), the three-layer model (3LM) aims at capturing wind farm–gravity wave interaction by combining a micro- and a meso-scale model. The well-known Bastankhah and Porté-Agel (2014) wake model is used at the micro-scale level, while a substantially improved evolution of the two-layer model (2LM) of Smith (2010) is adopted to model meso-scale effects produced by gravity waves. Although very promising due to the number of physical processes that it can capture, the 3LM suffers from two limitations. Firstly, the fact that it solves for depth-averaged velocities complicates its coupling with wake models, which are usually formulated in terms of the actual wind speed. Secondly, the coupling method only transfers the gross wind farm blockage at the micro-scale level, omitting the effect of gravity waves inside the wind farm and in its wake. The present section is organized as follows. Section 2.1.1 explains the meso-scale model used within the 3LM, Sect. 2.1.2 outlines the adopted wake model, Sect. 2.1.3 describes the turbulence intensity (TI) model that controls wake expansion within the wake model, and Sect. 2.1.4 explains the wake superposition strategy adopted in the original 3LM formulation. Finally, Sect. 2.1.5 outlines the technique used by Allaerts and Meyers (2019) to couple each sub-model, highlighting its issues and limitations.

2.1.1 Meso-scale model

The 3LM exploits the theory for interacting gravity waves and boundary layers developed by Smith (1980, 2007, 2010), with extra features such as the Coriolis force, the additional wind farm layer that relaxes Smith's homogeneous vertical mixing assumption, and the wind farm–gravity wave coupling mechanism. Specifically, the model is formulated by dividing the vertical ABL structure into three layers, i.e. the wind farm layer, the upper layer, and the free atmosphere or geostrophic layer aloft (a conceptual sketch is given in Fig. 1).

Figure 1Conceptual sketch of the three-layer model, reprinted with permission from Allaerts and Meyers (2019).

The depth-averaged Navier–Stokes equations used in the first two layers, linearized around a reference ABL state, are

(1) U 1 u 1 x + V 1 u 1 y + 1 ρ p x = - f c v 1 + ν t , 1 2 u 1 x 2 + 2 u 1 y 2 + D 11 H 1 u 2 - u 1 + D 12 H 1 v 2 - v 1 - C 11 H 1 u 1 - C 12 H 1 v 1 - f x H 1 U 1 v 1 x + V 1 v 1 y + 1 ρ p y = f c u 1 + ν t , 1 2 v 1 x 2 + 2 v 1 y 2 + D 21 H 1 u 2 - u 1 + D 22 H 1 v 2 - v 1 - C 21 H 1 u 1 - C 22 H 1 v 1 - f y H 1 U 1 η 1 x + V 1 η 1 y + H 1 u 1 x + v 1 y = 0


(2) U 2 u 2 x + V 2 u 2 y + 1 ρ p x = - f c v 2 + ν t , 2 2 u 2 x 2 + 2 u 2 y 2 - D 11 H 2 u 2 - u 1 - D 12 H 2 v 2 - v 1 U 2 v 2 x + V 2 v 2 y + 1 ρ p y = f c u 2 + ν t , 2 2 v 2 x 2 + 2 v 2 y 2 - D 21 H 2 u 2 - u 1 - D 22 H 2 v 2 - v 1 U 2 η 2 x + V 2 η 2 y + H 2 u 2 x + v 2 y = 0 ,

where subscripts 1 and 2 indicate the wind farm and upper layer respectively, η is the layer displacement, H is the layer height, U=(U,V) and u=(u,v) are the depth-averaged background and perturbation velocities respectively, νt is the eddy viscosity (see Sect. 2.3), and fc=2ωz is the Coriolis parameter, with ωz the component of the Earth's rotation rate vector normal to the sea-level geopotential surface approximated as a sphere (Gill1982). The term f/H1 in Eq. (1) is the wind farm forcing term and enables the wind farm–gravity wave feedback mechanism (see Sect. 2.1.5 for the definition of this term). Tensors Cij and Dij describe the perturbation of friction at the ground and at the top of the wind farm layer respectively, defined as

(3) C i j = τ | z = 0 U 1 3 2 U 1 2 + V 1 2 U 1 V 1 V 1 U 1 U 1 2 + 2 V 1 2


(4) D i j = τ | z = H 1 Δ U 3 2 Δ U 2 + Δ V 2 Δ U Δ V Δ V Δ U Δ U 2 + 2 Δ V 2 ,

where ΔU=(ΔU,ΔV)=(U2-U1,V2-V1) is the velocity difference, component by component, between the upper and wind farm layer, while ∥∥ denotes L2 vector norm.

The third layer, which is effectively a boundary condition, relates the total ABL vertical displacement to the perturbation pressure felt inside the ABL by

(5) 1 ρ 0 p ^ = i N 2 - Ω 2 m + g η ^

based on the linear, three-dimensional, non-hydrostatic, non-rotating gravity wave theory of Smith (1980) (see also Lin2007, and Nappo2012). The top-hat symbol is used to indicate that Eq. (5) is expressed in Fourier space, η=η1+η2 is the total vertical ABL displacement, N=γg/θ0 is the Brunt–Väisälä frequency, where the free-atmosphere lapse rate γ=dθ/dz is assumed constant in this study, g=gΔθ/θ0 is referred to as the reduced gravity, where Δθ is the potential temperature jump across the inversion layer, and θ0 is the reference potential temperature. The quantity m identifies the vertical wavenumber, evaluated using the dispersion relation

(6) m 2 = k 2 + l 2 N 2 Ω 2 - 1 ,

in which κ=(k,l) is the horizontal wavenumber and Ω=-κU3 is the intrinsic wave frequency.

After determining the background state (see Sect. 2.3), Eqs. (1), (2), and (5) are numerically solved on a 2D Cartesian grid using a Fourier–Galerkin discretization to yield the perturbation velocity u and displacement η in the wind farm and upper layers, as well as the perturbation pressure p. It is important to highlight the significance of the pressure variable, as this concept will be essential in the formulation of the MSC framework later on. First of all, Allaerts and Meyers (2019) make the assumption of zero vertical pressure gradient inside the ABL; hence p only varies in the horizontal directions. Secondly, looking at Eq. (5), it can be stated that p can only change in response to a vertical displacement of the ABL or if the latter experiences a capping inversion jump Δθ and/or a stable lapse rate γ. For this reason, the pressure variable considered by Allaerts and Meyers (2019) contains the effect of internal and interfacial gravity waves through the N and g terms, respectively, as well as the hydrodynamic pressure perturbation produced by the Ω term. In a truly neutral ABL, where N and g are zero, the hydrodynamic pressure perturbation is produced by the combination of individual turbine-level induction effects. However, as shown by Devesse et al. (2023), Eqs. (1), (2), and (5) are horizontally filtered to yield a horizontal resolution that is larger than the ABL height (they used a filter size of 1 km). As a result, the meso-scale model lacks the capability to resolve both the spatial extent and intensity of the pressure increase in front of each individual wind turbine, along with its impact on velocity, referred to as local blockage.

2.1.2 Wake model

The 3LM has to be complemented with a wake model in order to predict wind turbine thrust and power. In this regard, for each wind turbine, we start by considering a wake coordinate system ek having its origin at the turbine rotor center, the x axis aligned with the local free-stream velocity, the z axis directed as the turbine tower from base to top, and the y axis to form a right-handed frame of reference. Bastankhah and Porté-Agel (2014) have shown that the non-dimensional velocity deficit at a point P(x) in the wind turbine wake assumes a self-similar profile and can be expressed as

(7) U - U ( x ) U = W ( x ) = C ( x ) f ( r ( x ) , σ ( x ) ) H ( x ) ,

where U is the free-stream velocity magnitude, C(x) is the maximum normalized velocity deficit, r(x) is the radial distance between P and the wake centerline, σ(x) is a length scale related to the wake width, which increases with increasing downstream distance from the rotor, f(r,σ) is a similarity function describing the shape of the velocity deficit profile, and  is the Heaviside function. Different expressions exist for C(x) and f(r,σ) (see Lanzilao and Meyers2022, for a review), which identify the specific wake model.

In this study, following Allaerts and Meyers (2019), we use the formulas developed by Bastankhah and Porté-Agel (2014), corresponding to the Gaussian wake model, which evaluate f(r,σ) and C(x) as


where CT and D are the turbine thrust coefficient and diameter respectively. A known problem of the Gaussian wake model is that it does not conserve momentum when CT>8σ2/D2. Specifically, depending on the value of the CT coefficient, Eq. (9) can be undefined up to a certain downstream distance from the wind turbine. Such a limit in the model formulation can be disregarded if Eq. (7) is evaluated at wind turbine locations only, as real-life streamwise spacing generally satisfies CT<8σ2/D2. Conversely, if one aims at predicting the entire flow field using the wake model, a correction in the near-wake region is needed to avoid deficit saturation at high CT values. In the present study, we deal with this issue by using a damping function that gradually shifts from the Gaussian to the super-Gaussian wake model in the near wake (see Appendix A for details). The resulting wake model is fully equivalent to the Gaussian model after 4 downstream diameters, but it allows one to compute the velocity in the near wake by removing the singularity in Eq. (9). Similarly to Jensen (1983), Bastankhah and Porté-Agel (2014) express the wake width σ as a linear function of x, namely


where, according to Niayifar and Porté-Agel (2016), as=0.3837, bs=0.003678, and cs=0.2. The calculation of TI, namely the turbulence intensity at the wind turbine which is shedding the wake, is detailed in Sect. 2.1.3. Regarding CT, this is calculated with U using curves provided in Appendix B. Note that the concept of free-stream wind becomes ambiguous for waked wind turbines and will be discussed in Sect. 2.2.4.

2.1.3 Turbulence intensity model

A key parameter of any wake model is the wake expansion coefficient k*. Equation (12) shows that this parameter is related to the TI level at the wind turbine that is shedding the wake. This normally varies throughout the wind farm, with the strongest variation happening between the first row, where the background ambient TI is experienced, and the downstream rows, where the TI level is increased due to enhanced mixing operated by the wakes. In the present study, we adopted the TI model proposed by Niayifar and Porté-Agel (2016), where the actual turbulence intensity at each turbine location is given by


where TI is the background turbulence intensity provided as a user-defined input, while TI+,k is the added turbulence intensity, evaluated as

(13) TI + , k = max j H x k | j I j x k | j A w π R k 2 ,

where  is again the Heaviside function, xk|j is the distance between turbine k and turbine j rotor center, in the wake reference frame of turbine j, Rk is rotor radius of turbine k, and Aw is an area whose value goes from 0 to πRk2, depending on whether turbine j is waking k or not, with 0 corresponding to the non-waking situation (see Niayifar and Porté-Agel2016). The rationale behind Eq. (13) is that the value of TI+ at turbine k is the maximum among all TI values produced by other turbines, as if they were waking turbine k individually. The parameter Ij is an empirical function which, according to Crespo and Hernandez (1996), is given by

(14) I j ( x ) = d s 1 2 1 - 1 - C T , j e s TI f s x D j g s .

In the present study, we use the values proposed by Niayifar and Porté-Agel (2016) for es, fs, and gs, equal to 0.8325, 0.0325, and −0.32 respectively, while we vary ds in order to minimize TIk,LES-TIk,Mod, i.e. the L2 norm of the TI difference at each turbine location between our LESs (Sect. 3.3) and the model. This leads to ds=0.8798 instead of the original value of 0.73 proposed by Niayifar and Porté-Agel (2016). The reason for re-tuning the TI model comes from the fact that Eq. (14) starts to be valid after 5 downstream diameters from turbine k, which is exactly equal to the streamwise turbine spacing chosen for our LESs. This, combined with the choice to simulate an aligned wind farm, determined that the added turbulence intensity at each waked turbine is always produced by the turbine directly upstream, due to Eq. (13). As a consequence, the TI model operates at the edge of its validity bounds inside the wind farm, which we observed to produce poor results (in particular, underestimating TI+,k). In general, we argue that TI model tuning is always beneficial if data are available, as they allow running the wake model with representative TI values. This yields a more accurate representation throughout the wind farm of the power distribution relative to the first row. In fact, regardless of the wake model being used alone or coupled with meso-scale models, TI model tuning allows to the elimination of the a source of error that is not intrinsic to the wake model but instead arises due to erroneous TI input values. Finally, we emphasize that the same set of coefficients for the turbulence intensity model has been used throughout the present paper.

2.1.4 Wake superposition

When dealing with large wind farm arrays, turbine wakes can overlap, resulting in full or partial waking of downstream turbines. For this reason, individual velocity deficits have to be combined by means of wake superposition methods. Two main techniques have been extensively used:


where k denotes the kth wind turbine, Nt is the number of turbines in the array, uk is a velocity to be defined, associated with turbine k, and U is the free-stream velocity, assumed uniform for now, meaning that wake coordinate systems ek are equal for all k. While Eq. (15), first proposed by Lissaman (1979), consists in summing velocity deficits, Eq. (16) sums the squares of velocity deficits and was introduced by Katić et al. (1987). These studies made the additional assumption of uk=U, leading to a simple summation of W(x) or W(x)2. This produces deficits that are too high and has been addressed by Niayifar and Porté-Agel (2016) and Voutsinas et al. (1990), who used uk=U(xk) with xk, the locations of the kth wind turbine for Eqs. (15) and (16), respectively. This introduces some sort of weighting based on the velocity experienced at the turbine that is shedding the wake.

To be able to account for partial waking of downstream turbines, we follow the approach of Allaerts and Meyers (2019) and model each turbine using Nq=16 quadrature points, evaluating uk as

(17) u k = 1 N q q = 1 N q U x k , q ,

where xk,q indicates the coordinates of the qth quadrature point belonging to the kth turbine. This greatly enhances the wake model performance in partial waking or when a vertical velocity gradient is considered. In the current study, we use a cross-type quadrature points distribution rather than the star-type arrangement adopted by Allaerts and Meyers (2019) (see Fig. 2 for a comparison between the two distributions). Despite the cross-type arrangement only featuring two azimuthal directions (horizontal and vertical), it allows for a finer resolution along the rotor radius with the same number of quadrature points. This is, in our opinion, more suited to capture the effects produced by vertical velocity gradients and partial waking conditions.

Figure 2Comparison between the quadrature point distribution adopted in the present study and in Allaerts and Meyers (2019).


When presenting results from the original 3LM formulation of Allaerts and Meyers (2019), we use Eq. (15) to combine velocity deficits from the Bastankhah and Porté-Agel (2014) model. Conversely, the coupling methodology of the MSC model requires a different superposition method, which is described in Sect. 2.2.

2.1.5 Original coupling technique

In order to capture wind farm–gravity wave feedback, the meso-scale model and the wake model need to be suitably coupled. For conventional wake models, the thrust force associated with the kth turbine ultimately depends on the unperturbed free-stream velocity, i.e. fk=fk(U). In the coupling method proposed by Allaerts and Meyers (2019), fk is corrected for blockage effects based on a perturbation velocity uup so that fk=fk(U), where U=U+uup. The point-wise perturbation velocity uup is calculated by evaluating the meso-scale perturbation velocity u1(x) at xup, a point located upstream the first wind farm row, which is deemed representative of the global blockage effect for the entire farm. For a free-stream wind speed aligned with the x direction, Allaerts and Meyers (2019) evaluated the coordinates of the upstream point as xup=min(xk)-10D and yup=1/Ntk=1Ntyk, where xk and yk are the rotor center coordinates of turbine k. Finally, fk is filtered on the 3LM numerical mesh using a Gaussian kernel G. Such an operation returns a smoothed wind-farm body-force field and can be expressed as

(18) f ( x , y ) = 0 L x 0 L y G ( x - x , y - y ) k = 1 N t f k δ x - x k , y - y k d x d y ,


(19) G ( x , y ) = 1 π L 2 exp - x 2 + y 2 L 2 ,

where δ(x,y) is the Dirac function and L is the projection width, set to the greatest between the streamwise and spanwise meso-scale grid spacing.

Regarding the numerical algorithm used to solve the coupled problem, Allaerts and Meyers (2019) expressed fk's dependency on uup by expanding fk(U) in the Taylor series around U, truncating the polynomial at the first order. Such an operation required the Jacobian of the wake model – specific to the chosen model – with respect to uup to be derived. Moreover, it produced a convolution term in the Fourier transform of Eq. (1), yielding a system in Fourier space where wavenumbers are not mutually independent. In the present paper, we employ a slightly different but more flexible algorithm to solve the coupled problem that eliminates the need for the Jacobian computation. Specifically, we iterate between the meso-scale and the wake model, initializing the Gaussian-filtered thrust distribution by running the wake model with the uniform background wind U. This returns uk, which is calculated at each turbine location using Eqs. (7) and (15). At this point, turbine thrust is evaluated as

(20) f k = 1 2 u k 2 π D 2 4 C T , k e k ,

where ek is defined in Sect. 2.1.2, and CT,k is the thrust coefficient of the kth turbine. After filtering fk using Eq. (18), the 3LM is solved using Eqs. (1), (2), and (5), and the free-stream wind for the wake model is corrected as Ui+1=U+uup. The wake model is run again with the updated free-stream wind by setting U=Ui+1 in Eqs. (7) and (15), and the whole process is repeated until the relative residual between two consecutive iterations falls below a specified tolerance (we use the 3LM pressure to calculate the residual, i.e. r=pi+1-pi2/pi+12, where subscript 2 is the L2 norm). It is interesting to note that, in our approach, the Gaussian-filtered turbine thrust distribution resulting from the initialization step, namely Eq. (18), is equivalent to the Gaussian-filtered zeroth-order term of the Taylor series expansion proposed by Allaerts and Meyers (2019), while the iterative procedure removes the error initially represented by higher-order terms, correcting the entire thrust distribution at each iteration. Moreover, in the system of equations resulting from our procedure each wavenumber is independent of the remaining ones, leading to a block-diagonal system matrix (see Appendix C).

Although the employed numerical algorithm is slightly different, we emphasize that the physical coupling between the meso-scale and the micro-scale is identical to the one proposed by Allaerts and Meyers (2019). Such a coupling method is characterized by a non-local point-wise velocity correction that leads to the following structural limitations. First, the new free-stream velocity – although corrected for global blockage effects – is always the same for all turbines. Secondly, the upstream point location becomes ambiguous if an arbitrary wind farm geometry is considered or if the wind farm contains multiple sizes of turbines. Finally, such a coupling method completely disregards the beneficial effect of gravity wave-induced pressure gradients inside the wind farm.

2.2 The multi-scale coupled model

In the present section, we introduce the multi-scale coupled (MSC) model, a framework and coupling technique aimed at modeling turbine–wake interactions, local turbine blockage, wind farm–gravity wave interactions, and global blockage effects. Specifically, this is achieved by combining different sub-models and developing a suitable coupling technique that allows their inter-dependence. Strengths of the MSC model are its modularity and the fact that it builds upon already-existing sub-models, so that additional tuning parameters or individual sub-model re-tuning are not required. Ultimately, the MSC model aims to resolve the structural deficiencies of the 3LM model, presenting itself as a novel modular framework to capture the interaction between micro- and meso-scale phenomena in wind farms in a stratified atmosphere. The present section is organized as follows. In the remainder we will introduce the MSC model, detailing the new coupling method in Sect. 2.2.1. Then, Sect. 2.2.2 explains the reconstruction technique used to transform the depth-averaged perturbation velocity produced by the meso-scale model into the actual background velocity field required by the wake model. Since the obtained background wind is not uniform in space, we use a new wake superposition strategy, described in Sect. 2.2.3, that allows us to combine individual wake deficits with any heterogeneous background wind. Section 2.2.4 provides details on how local-induction effects are introduced, while Sect. 2.2.5 finally summarizes the MSC model solution procedure.

We start by defining a meso-scale and a micro-scale level, each associated with one or more respective sub-models. The meso-scale model is forced by the turbine thrust distribution, which is provided by the micro-scale model, while the latter is forced by a heterogeneous change (both in magnitude and direction) of the background velocity experienced at each turbine location. A conceptual sketch of this formulation is depicted in Fig. 3. Two assumptions underlie this approach. Firstly, we assume that turbine and atmospheric processes occur at separate spatial scales. As a result, the micro-scale sub-model only accounts for effects at the turbine scale (such as wakes or local blockage), while the meso-scale model is responsible for modeling the gross background velocity change produced by physical processes happening at the wind farm and/or ABL scale (such as wind farm–gravity wave interaction). Secondly, and most importantly, we assume that these effects can be linearly superimposed; i.e. the effective velocity field is given by the linear superposition of the meso-scale background wind speed and the wind speed changes predicted at the micro-scale level.

Figure 3Conceptual sketch of the MSC model. After step 0, used to initialize turbine thrust distribution, steps 1, 2, and 3 are solved iteratively, updating turbine thrust and background velocity, until the residual on perturbation pressure has converged to a specified tolerance.


In the present study, we use a modified version of the 3LM to model meso-scale processes, while the Bastankhah and Porté-Agel (2014) wake model and the Branlard and Gaunaa (2014) vortex cylinder models are used to capture micro-scale processes. To couple the micro- and meso-scales without double-counting their effects, both the superposition and the separation-of-scales assumptions must hold. At the micro-scale level, wake models are usually defined on the free-stream velocity for an isolated wind turbine wake, namely Eq. (7), to be consistent with the same definition used in theoretical power/CT curves. For this reason, except when ad hoc re-calibration is performed, they do not directly account for turbine–ABL interaction. Also the vortex cylinder model, derived using a vorticity formulation, was developed to specifically account for turbine induction, and its results only depend on turbine characteristics, such as thrust coefficient, diameter, and free-stream velocity at the turbine location. At the meso-scale level, on the other hand, the perturbation velocity predicted by the original 3LM contains both gravity wave and wake effects, as the equations are forced with the turbine thrust distribution. As a result, the original 3LM does not satisfy the superposition and separation-of-scales assumptions and needs to be suitably extended with an additional solution step. The objective is to exclude wake effects from the 3LM solution. Specifically, we recall that the 3LM pressure only accounts for atmospheric stratification and hydrodynamic effects at the large scale. Hence, after solving the original 3LM, i.e. Eqs. (1), (2, and (5) and obtaining p, we solve an additional modified set of equations, where perturbation pressure p is the forcing term and turbine thrust is removed. We denote such a procedure as the “three-layer model reconstruction” (3LMR) step, performed inside the wind farm and upper layer. The result is a so-called background perturbation velocity that only contains the effect of large-scale gravity-wave-induced and hydrodynamic pressure gradients.

2.2.1 Three-layer model reconstruction

Within the three-layer model reconstruction (3LMR) step, we use the perturbation pressure field obtained from the original 3LM to force a modified set of three-layer equations. It is important to note that, by fixing p, we are implicitly fixing the boundary layer displacement η due to Eq. (5). Hence, Eq. (5) is automatically satisfied for the third layer, and only the momentum equations in the wind farm and upper layers are retained, reading


where subscripts 1 and 2 again refer to the wind farm and upper layers, respectively, and ubk=(ubk,vbk) is the background perturbation velocity.

We highlight that perturbation velocities resulting from Eqs. (21) and (22) only contain the effect of large-scale gravity wave and hydrodynamic pressure gradients produced by the wind farm, while the direct influence of the latter is removed. However, since it is a depth-averaged perturbation quantity, ubk cannot be directly used within the wake model. In the next section, we deal with its conversion from depth-averaged to an actual, height-dependent velocity field.

2.2.2 Background velocity reconstruction

In order to transform the depth-averaged perturbation field resulting from the 3LMR step into an actual velocity field – also introducing height dependency – we developed a method that consistently matches the local background velocity profile Ub(x) with the layer average U1+u1bk in the wind farm layer. Following Panofsky (1963), who used the similarity hypothesis of Monin and Obukhov (1954), the velocity profile in the surface layer can be expressed as

(23) U b ( z ) = u * κ ln z z 0 - Ψ z L ,

where u* is the friction velocity, z0 is the equivalent roughness length, κ is the von Karman constant (which we set to 0.4 in the present study), and Ψ is a function of the Obukhov length scale L. For neutral stratification below the ABL height, Ψ→0, while for stable or unstable stratification different methods have been proposed to express Ψ, such as the method of Businger (1966) and Dyer (unpublished) for unstable conditions or the method of Etling (1996) for the stable ABL. In the context of this study, we only consider conventionally neutral ABLs, so we set Ψ=0. If the background flow exhibits small variations due to the effect of large-scale pressure gradients, it is reasonable to assume that the friction velocity u* varies horizontally, while the profile shape is preserved. In the present study, we evaluate u*(x,y) by matching the average velocity magnitude calculated exploiting Eq. (23) with results obtained from the 3LMR, namely

(24) U 1 + u 1 b k 2 + V 1 + v 1 b k 2 = u * ( x , y ) κ H 1 - z 0 z 0 H 1 ln z z 0 d z ,

which yields

(25) u * ( x , y ) = κ H 1 - z 0 U 1 + u 1 b k 2 + V 1 + v 1 b k 2 H 1 ln H 1 z 0 - 1 + z 0 .

To include the directional information, we assume that the background velocity profile used to initialize the 3LM, characterized in principle by an angle profile Φ(z), is rotated by a perturbation angle that only depends on the horizontal location, which we evaluate as

(26) ϕ ( x , y ) = arctan V 1 + v 1 b k U 1 + u 1 b k - arctan V 1 U 1

so that the overall flow angle at a given location is ϕ(x)=Φ(z)+ϕ(x,y). The background velocity profile inside the wind farm layer can finally be expressed as

(27) U b ( x ) = u * ( x , y ) κ ( ln z z 0 cos ( ϕ ( x ) ) i + ln z z 0 sin ( ϕ ( x ) ) j ) .

Note that this procedure can also be used in those cases where the temperature profile is stable or unstable below H by including expressions for Ψ and matching the integration numerically. Moreover, we note that such a reconstruction from depth-averaged to actual velocity only matches the vertical average of Eq. (23) in the wind farm layer. Adding the extra constraint of matching also the upper layer would result in an overdetermined system, which could still be solved in a least-squares sense. As this is out of the scope of the present study, we use Eqs. (25), (26), and (27) to reconstruct Ub(x) from u1bk.

2.2.3 New wake superposition

The framework described in Sect.2.2 requires a wake superposition method that accommodates the heterogeneous background velocity field Ub(x). Brogna et al. (2020) developed a variation of Eq. (16) that accounts for a variable background wind speed and direction, while Farrell et al. (2021) combined wake deficits using Eq. (15) with U=uk=Ub(x). Lanzilao and Meyers (2022) also developed new superposition techniques for combining wake deficits with a background flow with variable magnitude or both variable magnitude and direction. Although these methods are all compatible with the MSC framework, we introduce a very simple superposition strategy, which has the advantage of being equivalent to Eqs. (15) or (16) when the background flow is uniform in both speed and direction. This ensures that wake deficits are combined consistently with previous literature, avoiding the uncertainty on whether or not the wake model must be re-tuned.

First, we assume that the background flow angle varies on a much larger spatial scale than the turbine spacing. As a result, for a given turbine k, the wake is assumed to be aligned with the local background flow at the rotor, which is rotated accordingly. As a consequence, W=W(x|k), where |k means that x is expressed based on the local wake coordinate system of turbine k (see Sect. 2.1.2). In addition, we set U=∥Ub(x)∥. With the above variations, Eqs. (15 and (16) become


where Uw(x) now contains the combined effect of turbine wakes and large-scale pressure gradients. Equations (28) and (29) are different to the formulation of Farrell et al. (2021) in that uk is equal to the average of Uw(x) among quadrature points of turbine k, rather then the free-stream velocity at location x. Moreover, it is worth noticing that Uw(x) is, in reality, a velocity magnitude, as the directional information contained in Ub(x) only enters in the above equations by varying the direction of the wake shed by a given wind turbine. In the present study, we follow the combination of wake and turbulence intensity models used by Niayifar and Porté-Agel (2016) and Allaerts and Meyers (2019), who both used Eq, 15. In the context of the MSC model, to remain consistent with such a combination, we use Eq. (28) to combine wake deficits with the varying background wind Ub(x). Finally, we include ground effects by mirroring each wind turbine with respect to the ground plane (Nygaard et al.2022; Stevens et al.2015).

2.2.4 Local blockage model

As explained at the end of Sect. 2.1.1, the atmospheric perturbation model is incapable of resolving local blockage effects at the turbine level. As a consequence, these have to be suitably modeled at the micro-scale using a local blockage model that is compatible with the MSC framework. Specifically, the local blockage model should satisfy the superposition and separation-of-scales principles; i.e. it should be defined for an isolated wind turbine and should not satisfy mass conservation at the wind farm scale. Different models have been developed by the research community that satisfy these constraints (see Branlard et al.2020, for a review). In the present study, we choose to model local blockage using the vortex cylinder model (Branlard and Gaunaa2014; Branlard and Meyer Forsting2020). To this end, referring to the local wake coordinate system of turbine k, the axial perturbation velocity produced by turbine k at given point is given by

(30) u k v c ( r , x ) = γ t 2 π [ R - r + | R - r | 2 | R - r | + x k ( r , x ) 2 π r R ( K k 2 ( r , x ) + R - r R + r Π k 2 ( r , 0 ) , k 2 ( r , x ) ) ] ,

where x=x-xk is the axial distance from the rotor and r=(y-yk)2+(z-zk)2 is the radial distance from the wake axis (with xk, yk, and zk the coordinates of the rotor center of turbine k). R is the turbine radius, and K and Π are the complete elliptic integrals of the first and third kind respectively, with elliptic parameter

(31) k ( r , x ) = 4 r R ( R + r ) 2 + x 2 .

The momentum theory calibration relation (Branlard and Meyer Forsting2020) is used to express γt as

(32) γ t = - U b x k 1 - 1 - C T .

The total velocity at a point, containing wake and local induction effects, is then

(33) U ( x ) = U w ( x ) + k = 1 N t u k v c ( r , x ) .

Ground effects are included by mirroring the vortex cylinder with respect to the ground plane. Local induction at a given point is then computed using Eq. (33), i.e. by summing each turbine's individual effect, including images. Regarding such an operation it is important to notice that, at each coupling iteration, the final velocity U(x), containing meso- and micro-scale effects and averaged among each turbine quadrature points, is used to update turbine thrust, thrust coefficient CT, and power coefficient CP. These operations require that U(x) is the free-stream wind speed for the wind turbine of interest, which becomes a concept somewhat ambiguous in waked conditions. In this paper, we follow the approach of Branlard and Meyer Forsting (2020), where thrust, CT, and CP are evaluated using a free-stream velocity that contains induction effects from all cylinders except the ones associated with the current and mirrored turbines. We choose this method as it is consistent with what has been commonly done so far in the context of engineering models. We also tried another approach exploiting 1D momentum theory (in line with the methodology followed in most LES studies), but the results are similar (not shown here). In this approach, self-induction is retained in the model, which returns at this point the averaged disk velocity Udisk for each turbine. Then, the hypothetical free-stream wind speed is calculated as U=Udisk/(1-a), where a=(1-1-CT)/2 is the axial induction factor. This free-stream wind speed U is then used to update CT and CP coefficients from the turbine manufacturer curves. Finally, thrust and power are calculated using the disk-averaged velocity and the disk-based thrust and power coefficients proposed by Meyers and Meneveau (2010).

2.2.5 Solution procedure

This section details the overall MSC model solution procedure, which is also sketched in Fig. 4. The evaluation of the background state is described in Sect. 2.3 as it is not strictly part of the model and can be attained in different ways. As an initialization step, the wake model is first run using a uniform inflow velocity (optionally adding local induction effects), and initial thrust and CT are calculated for each turbine. At this point, the coupling loop is entered, and all models are solved in sequence at each iteration. Specifically, turbine thrust is filtered on the 3LM grid using Eq. (18); then Eqs. (1), (2), and (5) are numerically solved, yielding perturbation fields u1, u2, η, and p. In Fourier space, these equations represent a linear system with a coefficient matrix that is non-symmetric and block-diagonal (see Appendix C for details). Equations (21) and (22) are then used to solve for ubk, and Eqs. (25), (26), and (27) are used to reconstruct the background velocity field Ub(x) at turbine quadrature points from the 3LM grid by means of bilinear interpolation formulas. Successively, the wake model is run to yield the average of Uw(x) at the quadrature points. In particular, combining Eqs. (17) and (28) yields

(34) 1 N q q = 1 N q U w x i , q + j j i N t 1 N q q = 1 N q U w x j , q W x i , q | j = 1 N q q = 1 N q U b x i , q | j ,

which produces the linear system AuW=b, where


In Eqs. (34), (35), (36), and  (37), the notation xi,q|j denotes quadrature point q of turbine i, seen in the local wind coordinate system of turbine  j.

Figure 4Sketch of the MSC model solution procedure. We adopted the same colors of Fig. 3 to denote the micro-scale and meso-scale levels of the model.


At this point, local induction effects are added to uiW without considering turbine self-induction. Using Eq. (30) and averaging on quadrature points yields

(38) u i W + I = u i W + j j i N t 1 N q q = 1 N q u j v c r i , q , x i , q | j ,

where ri,q and xi,q|k are the radial and axial distances between the quadrature point q of turbine i and the rotor axis and center of turbine j respectively. Finally, the CT coefficient of each turbine is updated using uiW+I and the input CT curves, while turbine thrust is evaluated with Eq. (20), where we set uk=ukW+I. The coupling iteration is then repeated until the L2 norm of the relative pressure residual is lower than a specified tolerance. For instance, we noticed that around four or five iterations are needed to satisfy a specified tolerance of 10−4.

2.3 Background state

The MSC model requires information about the background ABL state, as well as the wind farm layout and turbine data. Regarding the first set of information, the required input parameters are reported in Table 1.

Table 1Input parameters for the MSC model.

Download Print Version | Download XLSX

In the context of the present study, we chose H1 as twice the average turbine hub height, while H2=H-H1. Reference potential temperature, inversion jump, and lapse rate may be prescribed based on reanalysis data such as ERA5, meso-scale simulations, or, if they are gathered in future wind resource campaigns, from observations. Latitude and roughness length can be estimated depending on the specific wind farm site. Average background velocities, wind angle, and shear stress magnitudes at the ground and at z=H1 require knowledge about the veered velocity and shear stress magnitude profiles. These can again be evaluated from observations (Allaerts et al.2018), from analytical models (see for example Nieuwstadt1983), or from numerical simulations. The same holds for background friction velocity and ambient turbulence intensity at the hub height, where the condition u*2=τ|z=0 has to be satisfied, as friction velocity and the shear stress at the wall are not independent. The depth-averaged eddy viscosity profiles can be calculated according to Nieuwstadt (1983) as

(39) ν t ( z ) = κ u * z 1 - z H 2 z < H 0 z >= H .

In the present study, as our aim is to perform verification of the MSC model against wind farm LES, we evaluate the background state using results from our conventionally neutral boundary layer (CNBL) precursor simulations. This avoids any mismatch in the initial background state, making sure that any discrepancy in turbine performance is due to the model only.

3 High-fidelity verification data

The present section describes methodology and results of the LESs used for the MSC model verification. First, the adopted numerical setup is described in Sect. 3.1. Then, Sect. 3.2 reports the choice of input parameters used for the CNBL simulations, characterizing the resulting background atmospheric state. The latter is used both to define the background ABL states required by the 3LM and MSC models and to prescribe the inflow condition for the wind farm analyses discussed in Sect. 3.3.

3.1 LES methodology

For the LESs presented in this paper, we use the open-source finite-volume code TOSCA (Toolbox fOr Stratified Convective Atmospheres) developed at the University of British Columbia. Details about the code, governing equations, and LES methodology are thoroughly explained in Stipa et al. (2024). For the sake of brevity, here we only report the adopted setup for the wind farm simulations.

To avoid gravity wave reflections, we use both periodic boundary conditions and a fringe region (Lanzilao and Meyers2023) located at the domain inlet. This also provides a suitable turbulent inflow, eliminating the wind farm wake re-advected at the inlet by periodic boundaries. At the upper boundary, we use a Rayleigh damping layer with a thickness of 12 km, i.e. slightly more than one expected vertical wavelength (this parameter can be calculated as λz=2πUg/N, where N is the Brunt—Väisälä frequency and Ug is the geostrophic wind). Lateral boundaries are periodic, implying that gravity waves induced by the wind farm will interact with their periodic images. This requires the domain to be sufficiently large for these interactions to happen far from and downstream of the wind turbines. Moreover, we use the advection damping technique developed by Lanzilao and Meyers (2023) to ensure that interactions between fringe-generated and physical waves are not advected downstream but remain trapped inside the advection damping region. After a reflectivity study that employed a computationally cheap canopy model (not shown here), we found that a Rayleigh damping coefficient of νRDL=0.05 s−1 and a fringe damping coefficient of νFR=0.03 s−1 yielded minimal gravity wave reflection. We used the same damping functions as Lanzilao and Meyers (2023), and in Table 2 we report their parameters for our wind farm simulations.

Table 2Fringe and advection damping region information.

Download Print Version | Download XLSX

The size of the successor domain is 40 km × 21 km × 28 km in the streamwise, spanwise, and vertical direction respectively, discretized with 1554×1194×345 cells. All directions are graded to reach a mesh resolution of 30 m × 12.5 m × 10 m around the wind farm. The concurrent precursor mesh coincides with the portion of the successor domain located inside the fringe region. As a consequence, it extends for 5 km × 21 km × 28 km. The mesh resolution in the streamwise direction is 15 m, while in the spanwise and vertical directions, it is the same as the successor. Detailed mesh information is reported in Table 3.

Table 3Mesh information for the wind farm cases.

Download Print Version | Download XLSX

To obtain a fully developed turbulence state in the concurrent precursor, cutting down computational cost at the same time, we use the hybrid offline/concurrent precursor technique described in Stipa et al. (2024). Specifically, a separate offline precursor is first conducted on a domain smaller than the fringe region for 1.2×105 s, and inflow slices are saved at each time step during the last 20 000 s of simulations. This first phase defines the background inflow conditions, and it is described in more detail in Sect. 3.2. Inflow slices are then periodized along the spanwise direction and mapped at the concurrent precursor inlet, as it uses inflow–outflow boundary conditions in this initial phase. After one concurrent precursor flow turnover time, inlet and outlet boundaries are switched to periodic, and the solution is then progressed for 30 000 s. Note that successor and concurrent precursor always proceed in sync with each other, and the former always employs periodic boundary conditions at the streamwise boundaries. After one successor flow turnover time (≈5000 s) flow statistics are gathered for 25 000 s.

The wind farm has a rectangular planform, with 20 rows and 5 columns. The first row is located at x=0 and extends from 300 to 2700 m in the spanwise direction. This determines a lateral spacing of 600 m (4.76 D), while streamwise spacing is set to 630 m (5 D). Wind turbines correspond to the NREL 5 MW reference turbine and are equipped with angular velocity and pitch controllers described in Jonkman et al. (2009). A very simple yaw controller is also added, which rotates wind turbines independently using a uniform speed of 0.5° per second when flow misalignment exceeds 1°. Flow angle is calculated by filtering the wind velocity at a sampling point located 1 D upstream of the rotor center, using a time constant of 600 s. Turbines are modeled using the actuator disk model (ADM) described in Stipa et al. (2024), while tower and nacelle are not included. The ADM force projection width is set to 18.75 m.

In both the concurrent precursor and successor simulations, velocity is controlled using a constant source term, obtained by averaging the offline precursor sources for the last 20 000 s. Temperature controller is retained in the concurrent precursor simulation, but it is switched off in the successor to allow the inversion height to be freely perturbed by the wind farm.

Table 4ABL parameters used for the finite wind farm simulation presented in this section.

Download Print Version | Download XLSX

3.2 CNBL precursors

In order to highlight the differences between subcritical and supercritical regimes, at the same time isolating the impact that atmospheric stability has on wind farm efficiency, we simulate two CNBL cases that only differ in the capping inversion strength but that have the same wind profile and free-atmospheric lapse rate. The choice of input parameters, summarized in Table 4 for both simulations, is based on the sensitivity study performed by Allaerts and Meyers (2019). We set the inversion strength of the subcritical case (N1) to 7.312 K, corresponding to Fr≈0.9, while we chose a value of 4.895 K, corresponding to Fr≈1.1, for the supercritical case (N2). The Coriolis parameter fc corresponds to a latitude of 41.33°.

The Rampanelli and Zardi (2004) model is used to initialize the potential temperature profile, where H is taken as the center of the capping inversion layer and Δh is the capping inversion width. The two CNBL precursors are advanced in time for 105 s; then data are averaged for 2×104 s. The domain size of the two CNBL simulations is 6 km × 3 km × 1 km in the streamwise, spanwise, and vertical directions respectively. The mesh has a horizontal resolution of 15 m, while in the vertical direction it is graded equally as the concurrent precursor and successor simulations. A driving pressure controller with geostrophic damping is used to fix the average velocity at href, chosen as the hub height, while a potential temperature controller is used to fix the average potential temperature profile throughout the simulation (both controllers and geostrophic damping use the same settings reported in Stipa et al.2024).

Figure 5 shows the resulting profiles of velocity magnitude, wind angle, potential temperature, and shear stress from simulations N1 and N2. In Table 5, quantitative data extracted from the two simulations are also reported. The capping inversion center H, ground temperature θ0, inversion strength Δθ, and inversion width Δh are calculated by fitting the Rampanelli and Zardi (2004) model in a least-squares sense. It is clear both from Table 5 and Fig. 5 that cases N1 and N2 have almost identical wind and shear stress profiles, while the background temperature stratification differs in the inversion strength, as explained earlier. Moreover, the pressure and potential temperature controllers allow us to obtain both an average hub-height wind speed and a ground potential temperature that match the input parameters listed in Table 4. Both cases result in the same average hub-height turbulent intensity of ≈9 %. According to the obtained CNBL profiles, it is clear that engineering wake models alone (optionally combined with a local induction model) would predict very similar power production for each individual turbine in the two cases. In reality, as will be shown later, wind farm interaction with the thermally stratified boundary layer leads to very different trends in power production within the farm.

Table 5ABL parameters obtained by fitting the Rampanelli and Zardi (2004) model for the subcritical (N1, Fr≈0.9) and supercritical (N2, Fr≈1.1) CNBL cases presented in this section, together with resulting friction velocity, minimum heat flux qmin within the boundary layer, and geostrophic wind angle.

Download Print Version | Download XLSX

Figure 5Velocity magnitude, wind direction, potential temperature, and shear stress profiles from precursor simulations.


3.3 Wind farm successors

In the current section, we briefly present results from our wind-farm large-eddy simulations. Figure 6 shows instantaneous fields of velocity magnitude at the hub height, for both the subcritical and supercritical cases, at the end of each simulation. High wake meandering can be observed in both cases, together with the gradual formation of a wind farm wake when moving downstream of the first row. However, from this figure it is not possible to appreciate the differences existing between the two simulations at the large scale. Therefore, it is instructive to look at the perturbation pressure field, which is evaluated by computing the difference between the pressure experienced in the successor domain at a given height and horizontal location and the horizontally averaged pressure from the concurrent precursor evaluated at the same height. This quantity is shown in Fig. 7, for the physical portion of the successor domain, i.e. excluding the fringe and advection damping regions. Substantial differences can be noticed between the two cases. In particular, the subcritical case (N1) shows a broader region of pressure increase with respect to the supercritical case (N2), with its maximum located around the first-row turbines. While in both cases such a positive pressure anomaly lasts for the first few wind farm rows, it extends for a longer distance upstream of the turbines in case N1. Additionally, the local and individual pressure increase in front of each wind turbine can also be observed. At the end of the wind farm, both cases show a negative pressure anomaly, less substantial than the positive one in front, and again more pronounced in the subcritical case.

Figure 6Instantaneous velocity magnitude at the hub height, sampled at the end of wind farm cases N1 (a) and N2 (b).


Figure 7Perturbation pressure at the hub height for wind farm cases N1 (a) and N2 (b).


Figure 8Perturbation pressure below the inversion layer for wind farm cases N1 (a) and N2 (b). The dashed line indicates the displacement of the streamline located at z=H at the inlet, magnified 5 times.


Inside the wind farm, a generally favorable pressure gradient is observed. The formation of these pressure anomalies can be explained by looking at the wind farm as an obstacle to the approaching flow. First, streamlines are displaced upwards and downwards at the wind farm entrance and exit, respectively. In addition, the boundary layer has to displace upwards after the first turbine rows in order to compensate for the mass flux deficit produced by wind turbine wakes. This results in a displacement of the inversion layer near the wind farm start that is larger than what is observed at the wind farm exit. Then, as the flow is vertically stratified starting from the inversion height, an upward perturbation of the flow particles corresponds to a positive pressure anomaly felt below, as a taller column of dense air is locally overtopping the ground. Furthermore, these vertical perturbations trigger interfacial and internal waves in the inversion layer and atmosphere aloft, respectively. As reported by Lin (2007), the interfacial wave amplitude depends on the Fr number. In subcritical conditions, these waves can propagate upstream, determining an increase of layer depth before the obstacle, thus augmenting the total energy of the flow. This cannot happen in supercritical conditions, as the interfacial wave speed is smaller than the advection speed. In such a condition, upstream propagation is produced by internal waves only, and smaller anomalies of the inversion layer in terms of both the amplitude and spatial scale are observed.

The positive correlation that exists between the boundary layer displacement and pressure perturbations is especially clear in Fig. 8. Here, perturbation pressure is shown in combination with the streamline displacement at the inversion height (magnified 5 times), which is a good estimation of the position of the inversion center locally (Allaerts and Meyers2017).

Moreover, it is clear that the pressure field around the wind farm is more complex than the mere superposition of individual turbine effects but is regulated by the interaction between the wind farm and the stratified boundary layer. It is reasonable to assume that such large-scale pressure disturbances around the farm also perturb the wind speed. Though such an effect is too small to be visualized in Fig. 6, it can affect wind farm power in different ways. First, the adverse pressure gradient upstream of the first row decreases the approaching wind and, depending on thermal stratification, different levels of global blockage will be experienced. Inside the wind farm, despite oscillations in the subcritical case, the pressure gradient is always favorable and counteracts the wind speed reduction. In the wind farm wake, especially for the subcritical case, the gravity waves clearly induce pressure oscillations, which may result in an apparently intermittent or non-monotonic wake recovery.

We emphasize that the two simulations share an almost identical wind profile and only differ in the background thermal stratification. The very different observed behavior, only due to differences in the potential temperature profile, is currently not captured by engineering models. In fact, even if modeling of the wind farm blockage had been attempted, for example, by superimposing individual induction effects via, for example, the vortex cylinder (Branlard and Gaunaa2014) or the Rankine half body (Gribben and Hawkes2019) approach, the models would erroneously produce identical results for wind farm cases N1 and N2. These LES results highlight the importance of accounting for atmospheric stability in engineering models.

4 MSC model results and verification

In this section, results from our wind farm LESs are compared against the original 3LM and the newly developed MSC model. For reference, we also show results obtained using the wake model alone (current industry standard), and the combination of the latter with the vortex cylinder model, to include local blockage effects. Verification of the model is performed by comparing its predictions against LES, focusing on both the flow variables, in Sect. 4.1, and on turbine quantities in Sect. 4.2.

In Table 6, we report the input parameters used to define the two background states for the meso-scale sub-model. The wind angle profile coincides with the one shown in Fig. 5. For the background height-averaged velocities in each layer, or to set the free-stream velocity for the wake model with and without local induction, we use the inflow velocity profile at the end of the advection damping region (the actual start of the meaningful portion of the successor domain). We observed that this velocity profile differs slightly (3.5 % at hub height) from the profile in the concurrent precursor simulation (this was also observed by Luca Lanzilao and Johan Meyers, personal communication, 2023). It is not clear why the fringe region induces such a modification of the input profile, and this is a topic for further investigation.

Table 63LM input parameters for cases N1 and N2.

Download Print Version | Download XLSX

As described in Sect. 2.2.3, wind turbines are mirrored in the MSC model, while in all other models we do not perform such an operation. The effect of including ground effects or not by mirroring wind turbines is shown in Fig. 9 for the wake model without induction and for the full MSC model. Data are normalized by the average power at the first row in order to remove the difference in first row power produced by neglecting blockage effects in the wake model.

Figure 9Effect of turbine mirroring on row-averaged power distribution for subcritical (a) and supercritical (b) cases. We show results for the MSC model and the wake model alone. Row-averaged power is normalized by its value at the first wind farm row to exclude the bias produced by the absence of blockage when using the wake model only.


The rationale for exclusively mirroring wind turbines in the MSC model is that while both this model and LESs incorporate stability effects, the potential temperature profile is not an input for the wake model. Consequently, if used independently, the wake model predicts strictly speaking the flow around a wind farm operating in a fully neutral atmospheric boundary layer, irrespective of whether the model is supplemented with a local induction parameterization or not. In this case, as shown by Lanzilao and Meyers (2024), the pattern of production is similar to the one presented for the wake model with turbine mirroring in Fig. 9. Conversely, when stability is modeled within the LES and MSC framework, it induces a favorable pressure gradient inside the wind farm due to the flow being confined below the inversion layer (Smith2024). This effect, which enhances wake recovery within the wind farm, is implicitly accounted for by the wake model alone if mirroring is not performed since wake effects are reduced. For this reason, the adopted wake model as defined by Niayifar and Porté-Agel (2016) captures the pattern of production with more accuracy when turbine mirroring is not performed. In general, it should be considered that, once a model configuration is selected (i.e. number of rotor quadrature points, superposition method, wake model, mirroring or no mirroring, turbulence intensity model) and model tuning is performed, it is very possible that two different configurations may achieve similar results with distinct sets of tuning coefficients. These considerations uncover a confounding artifact of engineering models: when an already tuned model is supplemented with additional physics or features, it may achieve poorer predictive accuracy and it would need to be re-tuned, since the original coefficients or its configuration choices might be already accounting for some of the missing physics. In fact, as can be noticed from Fig. 9, the trend of row-averaged thrust and power distributions with respect to the first row is better captured in the wake model without mirroring effects. This is not surprising, as this model configuration is identical to Niayifar and Porté-Agel (2016) apart from the turbulence intensity model, whose tuning is not aimed at converging turbine power, but rather at providing the correct TI for the TI–wake expansion relation in the wake model. Conversely, the addition of turbine mirroring produces increased wake effects, leading to thrust and power distributions that are monotonically decreasing with downstream distance from the first few rows. This is a feature of fully neutral conditions, where the absence of stratification does not produce any strong favorable pressure gradient throughout the wind farm. In relation to the LES data of the present study, it can be seen how the wake model overpredicts thrust and power at the first rows, as it neglects blockage effects. On the other hand, the absence of a stratification parameterization leads these quantities to be strongly under predicted towards the wind farm end when using turbine mirroring. Interestingly, we observed that the increased wake effects produced by rotor images tend to balance the overestimation given by the absence of blockage at the first rows, reducing the error in overall wind farm power. This is a lucky situation, where two model issues compensate each other to reduce the global error while completely misrepresenting the actual shape of thrust and power distributions. Regarding the MSC model results of Fig. 9, it can be noticed how the stability-induced favorable pressure gradient between the first and last wind farm rows causes the background velocity to accelerate within the wind farm. Although the magnitude of this effect depends on the specific conditions, it results in a positive shift in power towards the wind farm end if turbine mirroring is not applied. Nevertheless, since the MSC model is now capturing this missing physics, increased wake effects produced by turbine images are not an issue for this case and ground effects can be modeled without affecting the shape of power and thrust distribution.

Regarding the 3LM, we use a numerical grid of 400×203 km in the streamwise and spanwise direction respectively, discretized with a uniform 500×500 m grid. As the NREL 5 MW thrust curve is not available in official literature, we ran several large-eddy simulations of isolated ADMs with uniform, non-turbulent inflow to compute the thrust and power curves used in the reduced models for the present analysis (see Appendix B). In both the MSC model and original 3LM runs, we use five coupling iterations, which are sufficient to bring the pressure residual below 10−4. We noticed that only one iteration is sufficient to capture the gross wind farm–gravity wave interactions (results not shown here), while three iterations provide thrust and power distributions within 0.5 % of the value obtained using five iterations.

4.1 Flow field verification

In this section we compare the flow quantities predicted by the reduced-order models against the LES results.

Figure 10Depth-averaged perturbation velocity in the wind farm layer, obtained after solving Eqs. (1), (2), and (5) for subcritical (a) and supercritical (b) cases. Background perturbation velocity obtained after solving Eqs. (21) and (22) for subcritical (c) and supercritical (d) cases. Velocity profile produced by the average matching procedure at the center of the first (e) and last (f) wind farm rows. Subcritical and supercritical cases are displayed in the same diagram. Dashed lines highlight the hub height hhub, while dotted lines indicate hhub±D/2.


Figure 10 shows, for both subcritical and supercritical conditions, the depth-averaged perturbation velocities in the wind farm layer, calculated after solving the 3LM step – Eqs. (1), (2), and (5) – and the 3LMR step – Eqs. (21) and (22). The background perturbation velocity field u1bk is then compared with vertical profiles of Ub(x), obtained after applying the average matching procedure – Eqs. (25), (26), and (27) – for two locations of interest, i.e. the middle of the first and the last wind farm rows. First, it can be seen from the reconstructed hub-height velocity that the region of maximum flow deceleration is located more upstream in the subcritical than the supercritical case. However, the difference between the two regimes is more evident in the downstream direction, as the supercritical case is characterized by several distinct V-shaped waves generated by the principal sources of boundary layer displacement, i.e. the wind farm entrance and exit. Conversely, in the subcritical case, the wind farm excites resonant lee waves, as the vertical wavenumber m becomes imaginary. These is in agreement with Allaerts and Meyers (2019) and Devesse et al. (2022). Moreover, in both cases the wind farm wake is deflected towards the negative y direction by momentum exchange between the wind farm and upper layers. Another interesting aspect is the complex and non-symmetric wave pattern that arises to either side of the wind farm due to the interaction of two different wave trains departing from the wind farm entrance and exit respectively. This is a consequence of the wind veer produced by the Coriolis force, which causes a misalignment between the geostrophic wind and the hub-height velocity. As a consequence, as interfacial and internal waves are formed around and above H, respectively, there exists a spanwise shift between the two wave trains, approximately given by Lxsin (Φ(H)), where Lx is the wind farm length in the x direction and Φ(H) is the flow angle with respect to the x axis at the inversion center. This phenomenon is very clear in Fig. 10c, where the interaction between resonant lee waves is stronger on the southern rather than the northern side of the wind farm. Moreover, looking at the reconstructed profiles of Ub, it can be noticed how the subcritical case experiences more blockage than the supercritical case at the wind farm entrance (Fig. 10e), while the background velocity for case N1 is higher than case N2 at the exit of the cluster (Fig. 10f).

Figure 11 displays the comparison of the predicted perturbation pressure and velocity in the wind farm layer between the LES and MSC model, spanwise-averaged from y=0.3 to y=2.7 km. Regarding the 3LM step within the MSC (equivalent to the original 3LM formulation), perturbation pressure agrees quite accurately with LES results, while depth-averaged perturbation velocities are overestimated both in the wind farm and upper layer (the latter is not shown).

Figure 11Left: perturbation pressure inside the boundary layer given by the 3LM within the MSC and from LES results. LES data are averaged within the upper layer only, to remove the effect of local pressure increase in front of each turbine, not present in the 3LM. Right: depth-average perturbation velocity in the wind farm layer.


As can be noticed, the model predicts a stronger velocity deficit with increasing downstream distance from wind farm start. The same behavior has been observed also by Allaerts and Meyers (2019), who argued that the 3LM underpredicts shear stress, especially at the interface between the wind farm and upper layer, as it neglects the effect of added turbulence intensity produced by the wind farm. While such limitations in capturing velocity deficit and shear stress are topics for further research, we note that they have almost no influence on the reconstruction step, as only the pressure field is used to compute depth-averaged, large-scale velocity variations. In fact, we observed that pressure-induced wind speed perturbations are small if compared to the ones produced by the wind farm, leading to a negligible magnitude of shear stress terms in the 3LMR step (not shown here). These are mainly retained for numerical reasons, as they improve the conditioning number of the linear system matrix.

Figure 12 compares LES results against the solution obtained at the final step of the MSC model, i.e. the superposition of Ub(x) and the analytical wake and induction models. For reference, the velocity resulting from the sole combination of the wake and vortex cylinder models is also shown. First, it can be noticed how the latter approach produces in practice the same results in both subcritical and supercritical conditions, as they only differ slightly in their inflow profiles. Conversely, in the MSC model, large-scale pressure perturbations greatly influence the micro-scale velocity, making it possible to capture gravity wave flow patterns around and inside the wind farm.

Figure 12Micro-scale flow field obtained from LES (a, d), MSC model (b, e), and wake model with vortex cylinder (c, f). Top and bottom diagrams correspond to the subcritical and supercritical case, respectively.


In Fig. 13, we compare data from Fig. 12 against LES results in more detail, for all reduced models. Specifically, we display the row-averaged (spanwise-averaged between the row start and end locations) hub-height velocity for the subcritical case only (data from the supercritical case are not reported, as they point to the same conclusions). Regarding the original 3LM formulation, such an analysis was not performed by Allaerts and Meyers (2019), but it is interesting as it allows one to visualize what happens inside the wake model and how blockage effects have been introduced in the original 3LM formulation. In the wakes-only approach, the free-stream velocity is constant up to the first wind farm row, as it is not influenced by meso-scale effects. Here it presents a discontinuity due to the presence of the first wake. The same holds for the original 3LM, but the free-stream velocity is now reduced due to the 3LM coupling through the sampling point upstream the wind farm. Both models show velocity jumps at the turbine locations due to the discontinuity at the disk given by the wake model. In general, they both predict a higher velocity than LES at the turbine locations, which results in an overestimation of thrust and power, as it will be shown later. Regarding the wake model coupled with the local induction model, it still underpredicts global wind farm blockage, but the velocity distribution is more accurate inside the wind farm. Nevertheless, it is worth remembering that, when using this approach, turbine self-induction is not accounted for in estimating the free-stream velocity for a given turbine. Furthermore, local induction effects given by the vortex cylinder model at a given turbine from all remaining turbines are very limited (e.g. blockage is ≈1 % at 5 diameters upstream, with a CT of 0.9).

Figure 13Comparison of row-averaged wind between LES and different reduced models. Blue: wake model, orange: wake and local induction model, green: original 3LM, red: MSC model.


For this reason, although including local induction effects improves velocity results, it only produces a negligible improvement in the actual thrust and power estimates if compared with the wakes-only approach. Looking instead at the MSC model, it can be noticed how large- and small-scale blockage effects are both accurately captured, together with velocity oscillations in the wind farm wake, operated by the large-scale gravity-wave-induced pressure gradient. Additionally, the model is also able to satisfactorily capture the velocity distribution inside the wind farm, with the largest error observed just after first row.

Focusing on the region directly upstream the wind farm, Fig. 14 displays the row-averaged percentage velocity reduction for the MSC model and LES results only, from both subcritical and supercritical conditions, highlighting the ability of the proposed model to capture the different amount of blockage produced in each case.

Figure 14Comparison between row-averaged velocity deficit upwind of the wind farm, obtained from LES and using the MSC model, for subcritical and supercritical conditions.


4.2 Turbine data verification

Regarding turbine quantities, Fig. 15 shows comparisons of row-averaged thrust and power distributions between LES and the models. First of all, looking at LES data, it can be noticed that a significant difference exists between subcritical and supercritical conditions. In the first case, pressure disturbances produced by the resonant lee waves induce velocity variations at a comparable scale to the wind farm streamwise length Lx. This results in a fluctuation of thrust and power that is characterized by a wavelength roughly equal to Lx, not observed in the supercritical regime. Interestingly, Fig. 15 highlights how the wake model alone consistently overestimates both thrust and power. Moreover, results are negligibly affected when local induction effects are added. At the same time, the wake model is indeed able to capture the global trend in power and thrust distributions, and it produces quite accurate results if diagrams are normalized by power and thrust experienced at the first row (not shown here). This highlights how absolute values are affected by a global bias that arises from failing to capture the effective velocity at the first row, i.e. the effect of large-scale gravity wave blockage. This concept of “an overprediction bias that pervades the entire wind farm”, previously hypothesized by Bleeg et al. (2018), is further corroborated by our results.

Figure 15Comparison of row-averaged thrust and power distributions for subcritical (a, c) and supercritical (b, d) cases. Blue: wake model, orange: wake and local induction model, green: original 3LM, red: MSC model.


The original 3LM formulation of Allaerts and Meyers (2019) was the very first attempt to solve such structural deficiency of wake models. The approach does capture a certain amount of blockage, as the free-stream velocity used by the wake model is uniformly reduced according to the perturbation velocity obtained from the 3LM at the upstream sampling point. We highlight that, in this approach, the wake model is still run using a horizontally uniform inflow velocity; thus it is not able to capture any trend inside the wind farm. As a result, thrust and power distributions are simply shifted to lower values as compared to the wakes-only approach with or without local induction. Moreover, the amount of blockage that the original 3LM is able to capture is very sensitive to the sampling location, which is clearly too far from the first row for the supercritical case. Finally, looking at results from the MSC model, it is evident how it can accurately capture both large-scale blockage effects and wind farm–gravity wave interactions, showing very good agreement to the LES results. A very interesting aspect to note, which is captured by the MSC model, is the combined effect of the unfavorable pressure gradient upstream the wind farm and the favorable gradient inside. The magnitude of these effects and the prevalence of one over the other depend on many parameters such as the inversion height, strength, and lapse rate. In our simulations, the first wind farm row is characterized by a lower power generation in subcritical than supercritical conditions, as gravity-wave-induced blockage is higher. On the other hand, subcritical conditions lead to stronger favorable pressure gradients than supercritical inside the farm, resulting in a higher power production of waked turbines. As a result, the wind farm produces similar power for the two cases analyzed here, as shown in Table 7. This emphasizes that in addition to modeling the blockage effects, properly capturing the large-scale interactions between wind farms and the atmosphere is essential to accurately estimate wind farm power. In fact, depending on the specific stratification conditions, the combination of these effects may hamper or enhance power production as compared to fully neutral conditions.

Table 7Total wind farm power obtained from LESs and as predicted by reduced-order models.

Download Print Version | Download XLSX

In Fig. 16 we show both the model error with respect to LES (both cases are plotted together) for each individual turbine and the model error on global wind farm power for each case. Looking at the first metric (left panel), it can be noticed how the wake model with and without local induction effects and the original 3LM consistently overestimate wind farm power, both at the first row (isolated points at the top-right of the diagram) and inside the wind farm, with the 3LM best performing among the three. Moreover, all these models are not able to capture the effects of wind farm–gravity wave interaction on the waked turbines, resulting in a circular cloud of points shifted upwards with respect to the zero-error line. Conversely, the MSC model not only substantially improves predictions at the first row, but also captures power and thrust oscillations inside the wind farm. As a consequence, points are more distributed along the zero-error line, proving that the proposed framework can model the actual physical processes produced by thermal stratification. Regarding the second metric, i.e. global wind farm power prediction, the MSC model is substantially more accurate than all other models (2.7 % and 1.6 % error for the subcritical and supercritical case, respectively), additionally underestimating slightly wind farm power (a desirable effect for industrial applications). Conversely, using any of the other models always overestimates wind farm power by at least 10 % for this high thrust scenario, with the wakes-only approach being the least accurate (14 % and 19 % error for the subcritical and supercritical cases, respectively). Figure 16 also confirms that including local induction effects improves wake-model predictions only marginally.

Figure 16(a) Error on power for each turbine considered individually. Data from both subcritical and supercritical cases are shown together. (b) Global wind farm power error of each model in percentage of LES. Blue: wake model, orange: wake and local induction model, green: original 3LM, red: MSC model.


5 Conclusions

The purpose of the present study was to introduce a fast, reduced-order model capable of capturing the interaction between wind farms and the thermally stratified atmosphere, along with local induction and wake effects. We proposed the multi-scale coupled (MSC) model, which combines a heterogeneous background velocity field, computed numerically using a meso-scale model, with analytical results from wake and local induction models.

To verify the model's accuracy, comparisons were made against finite wind farm large-eddy simulations (LESs) characterized by similar velocity and shear stress profiles but different potential temperature structures. Due to the high amount of computational resources required to conduct the above LES studies, model verification has only been performed against two atmospheric states. Specifically, while a neutral stratification within the ABL and a free-atmosphere lapse rate of 1 K km−1 have been prescribed in both cases, the inversion strength has been varied to produce a subcritical regime, where interfacial waves can propagate upstream, and a supercritical regime, where upstream propagation is only performed by internal waves. The background atmospheric states for both the original three-layer and MSC models have been computed from LES results, while turbine thrust and power curves have been evaluated by running several non-turbulent isolated turbine simulations. Furthermore, the turbulence intensity model used in the wake model was adjusted to match the turbulence intensity observed in LES.

The LES results revealed substantial differences in the underlying physics between the two regimes, despite observing a similar total wind farm power. In the subcritical case, characterized by a higher blockage and lower first-row power, gravity waves induced a stronger favorable pressure gradient within the wind farm compared to the supercritical case. Consequently, the waked turbines in the subcritical regime experienced a faster wake recovery. These opposing effects led to comparable overall wind farm power in both cases.

Regarding the MSC model, we found that the 3LM step strongly overpredicts the depth-averaged velocity deficit in the cluster wake. We believe that this is due to a poor modeling of the shear stress terms within the 3LM, suggesting that further research is required in this direction. Nevertheless, this issue has negligible impact on the background velocity field resulting from the reconstruction step. In fact, only the large-scale perturbation pressure, accurately predicted by the 3LM, is used in the 3LMR step to force the three-layer equations. In addition, comparisons of wind turbine thrust, power, and wind velocities demonstrated the MSC model's ability to accurately capture the intricate interaction between wind farms and the atmosphere. In contrast, the wake-only approach, with or without local induction effects, failed to differentiate between the two regimes and produced results solely dependent on the background velocity profile, disregarding global blockage effects. This introduced a positive wind speed bias that propagated downstream, resulting in an overestimation of power throughout the wind farm. Although the original 3LM formulation partially addressed this issue by correcting the free-stream velocity used in the wake model, it failed to capture the beneficial pressure gradient effect within the wind farm due to the uniform and non-local velocity correction. Nevertheless, the model still overestimated power within the wind farm, particularly in supercritical conditions where the wind correction became minimal due to the large upstream sampling distance of 10 rotor diameters, leading to an underestimation of blockage effects.

Overall, if compared against LES, the MSC model underestimates wind farm power by ≈2 %, while all other approaches overestimate power from ≈10 % to ≈19 %, with the original 3LM and the wake-only approach without local induction being the most and least accurate respectively.

We believe that the coupling framework introduced in the MSC model could present a significant contribution towards advancing our understanding of the complex physical phenomena arising in large-scale wind farms. Additionally, it may help refine the power estimation methodologies for existing and future wind farm clusters, enhancing the accuracy and reliability of power predictions.

Appendix A: Near-wake correction for Gaussian wake model

As pointed out in Sect. 2.1.2, the Bastankhah and Porté-Agel (2014) Gaussian wake model is limited in predicting wake deficit in the near wake. In particular, Eq. (9) becomes singular when CT>8σ2/D2. Hence, based on the value of CT, the singularity region may extend more or less behind the rotor, or be completely non-existent. Nevertheless, as values of CT in below-rated turbine operation conditions range from 0.7 to 0.9, a singularity region is to be expected in many cases. Moreover, since wake width σ also depends on CT through Eq. (10) in a non-linear way, the specific distance at which Eq. (9) turns non-singular cannot be evaluated unless a numerical method is used. Additionally, C(x) would be zero at such a location, which is of course non-physical.

To overcome this set of issues, we propose a near-wake correction for the Gaussian model, by smoothly transitioning from the super-Gaussian wake model of Shapiro et al. (2019) to the Bastankhah and Porté-Agel (2014) profile in the far wake through a damping function. The resulting model is fully equivalent to the Gaussian model at downstream distances greater than 4 diameters, while it approaches a top-hat distribution close to the rotor. In addition to the latter property of the super-Gaussian wake model, the resulting deficit function possesses the desirable feature that the maximum deficit is not located directly past the rotor disk, but somewhere in the near wake (Blondel and Cathelain2020).

According to Shapiro et al. (2019), the functions f(r,σ) and C(x) in the super-Gaussian wake model are given by


where Γ is the gamma function. The coefficient n, which is a function of the downstream distance from the rotor, determines the deficit shape. When n=2, the model is fully equivalent to the Gaussian wake model, while C(x) tends to a top-hat distribution for higher values of n. To enhance the performance of the model in the near wake, we evaluate n as follows:


where nn and wn are the exponent and weight associated with the near wake, while nf and wf pertain to the far wake. The parameters xc and δ are set to 2 and 4 D respectively. Figure A1 shows the evolution of nn, nf, and n, as well as the weights associated with the near- and far-wake regions. As can be noticed, n has already fell down to 2 at 4 downstream diameters, while it increases approaching the rotor so that wake deficit tends to a top-hat distribution. On the other hand, the damping functions associated with the weights cross at 2 D and are characterized by a merging region of 4 D, which is controlled by parameters xc and δ.

Figure A1(a) Streamwise evolution of overall, near-wake, and far-wake components (n, nn, and nf respectively) and (b) near- and far-wake weighting function evolution.


Figure A2(a) Flow field at the hub height predicted by the super-Gaussian, Gaussian, and Gaussian with near-wake correction, for a wind speed of 8 m s−1 and a thrust coefficient of 0.9. (b) Streamwise evolution of the wind velocity at the wake centerline for the three models, for three different values of the thrust coefficient. The singularity region of the Gaussian wake model is evident from the diagrams.


Figure A2 compares velocity evolution along the wake centerline, for the super-Gaussian and Gaussian models, the latter with and without the proposed near-wake correction. On the left panel, the whole velocity field is shown for a free-stream wind speed U of 8 m s−1 and a CT coefficient of 0.9, while, in the right panel we show velocity evolution at the wake centerline for three different values of CT, at the same wind speed. As can be noticed, the Gaussian model is singular below roughly 1.5 D, for CT=0.9, whereas a null wind speed is reached slightly below 1 D with CT=0.6. For the super-Gaussian model we used the coefficients proposed by Blondel and Cathelain (2020). These lead to a stronger deficit in the far-wake region and a minimum wind speed located around 2.5 D in all cases.

The Gaussian model with near-wake correction predicts a slightly lower velocity in the near wake than the super-Gaussian model, while the maximum deficit is located around 2 D. After 4 D, the predicted velocity exactly overlaps data obtained with the original Gaussian wake model, highlighting how the two produce equivalent results when turbine spacing is greater than such a distance.

Appendix B: Turbine data for the NREL 5 MW

To model wind turbines in our LESs, we use an actuator disk model that takes into account the detailed rotor information instead of just using prescribed turbine CT from curves, which are used instead to evaluate thrust and power in reduced-order models. This introduces an extra source of uncertainty in comparing low- and high-fidelity data, as a power or thrust difference could be both due to structural model deficiencies or simply due to the fact that different turbines (CT or CP) are used for a given wind speed. Moreover, for the NREL 5 MW turbine, a reliable CT vs. wind speed curve is not available in the literature. For instance, the thrust curve provided in Jonkman et al. (2009) also contains rotor inertial forces, resulting in a CT vs. wind speed curve that is biased towards higher values.

These considerations motivated us to evaluate CT and CP curves for the NREL 5 MW turbine using LES, thus removing any source of bias when comparing reduced-order model results against high-fidelity data from LES. Specifically, we run 16 isolated turbine simulations, spanning wind speed magnitudes from 5 to 20 m s−1, using the same turbine model adopted for the wind farm simulations. The inflow profile is characterized by a uniform velocity, while the outlet boundary condition is set to zero gradient on all fields. Top and bottom and spanwise boundaries are set to periodic. The domain extends 10 D from the rotor center in all directions and has a uniform spacing of 10 m. The simulations have been advanced in time for 600 s, and data have been averaged for the last 300 s. The Gaussian projection width has been set to 18.75 m.

Figure B1 shows turbine quantities regulated by the control system, i.e. rotor omega, blade pitch, and generator torque, alongside the resulting CT and CP curves.

Data from the right panel of Fig. B1 have been used to update CT and CP values at each MSC model coupling iteration, as well as in all other reduced model runs reported in the present paper.

Figure B1(a) Rotor omega, blade pitch, and generator torque as a function of wind speed, resulting from the control system action. (b) Thrust and power coefficient curves as a function of wind velocity.


Appendix C: 3LM/3LMR solution

In order to solve Eqs. (1), (2), and (5), we first transform them in Fourier space. This leads to a 7×7 linear system for each k and l wavenumber couple, with k and l the wavenumbers associated with the x and y directions, respectively. After some manipulation, displacements η1 and η2 can be eliminated from the equations, and pressure can be directly expressed as a function of the perturbation velocity components in the wind farm and upper layers, reducing the number of equations to 5 for each wavenumber couple. To speed up the numerical solution, linear systems associated with each wavenumber couple are combined to form a unique system, characterized by a block-diagonal solution matrix A3LM of size 5NkNl×5NkNl, where Nk and Nl are the wavenumber sizes in the streamwise and spanwise directions respectively. We use a sparse matrix format to express A3LM, while the complex linear system is solved using the generalized minimum residual (GMRES) (Saad and Schultz1986) iterative method. The full linear system, corresponding to the 3LM solution, is given by

(C1) A k 1 , l 1 A k 1 , l 2 A k N k , l N l s k 1 , l 1 s k 1 , l 2 s k N k , l N l = b k 1 , l 1 b k 1 , l 2 b k N k , l N l ,


(C2) A k i , l j = i σ 1 + ν t , 1 | κ | 2 + D 11 + C 11 / H 1 - D 11 / H 1 f c + D 21 + C 21 / H 1 - D 21 / H 1 - D 11 / H 2 i σ 2 + ν t , 2 | κ | 2 + D 11 / H 2 - D 21 / H 2 f c + D 21 / H 2 k i H 1 Φ / σ 1 k i H 2 Φ / σ 2 - f c + D 12 + C 12 / H 1 - D 12 / H 1 i k i / ρ i σ 1 + ν t 1 | κ | 2 + D 22 + C 22 / H 1 - D 22 / H 1 i l j / ρ - D 12 / H 2 - f c + D 12 / H 2 i k i / ρ - D 22 / H 2 i σ 2 + ν t , 2 | κ | 2 + D 22 / H 2 i l j / ρ l j H 1 Φ / σ 1 l j H 2 Φ / σ 2 1 / ρ


(C3) s k i , l j = u ^ 1 k i , l j u ^ 2 k i , l j v ^ 1 k i , l j v ^ 2 k i , l j p ^ k i , l j , b k i , l j = f ^ x k i , l j f ^ y k i , l j 0 0 0 .

In the above equations, f^x and f^y are the Fourier transforms of the Gaussian-filtered wind farm force defined in Eq. (18), Φ is the complex stratification coefficient, defined as the complex quantity between square brackets in Eq. (5), |κ|2 is equal to (ki2+lj2), and σ1 and σ2 can be evaluated as (kiU1+ljV1) and (kiU2+ljV2) respectively and are referred to as the intrinsic frequencies in the wind farm and upper layers (Smith2010).

For the 3LMR step, i.e. Eqs. (21) and (22), the procedure is identical upon elimination of the last row and column from A(ki,lj). By doing so, the size of A3LMR becomes 4NkNl×4NkNl. The unknown and right-hand side of the linear system associated with the 3LMR step are defined as

(C4) s k i , l j = u ^ 1 b k k i , l j u ^ 2 b k k i , l j v ^ 1 b k k i , l j v ^ 2 b k k i , l j , b k i , l j = - i k i p ^ k i , l j / ρ - i l j p ^ k i , l j / ρ - i k i p ^ k i , l j / ρ - i l j p ^ k i , l j / ρ .

After solving the two linear system associated with matrices A3LM and A3LMR in cascade, u^1bk components are transformed back into physical space, and Eqs. (25)–(27) are used to reconstruct the background velocity field Ub(x). The latter is then interpolated at turbine quadrature points from the 3LM grid by means of bilinear interpolation formulas.

Code availability
Data availability

The dataset can be made available from the authors upon request.

Author contributions

Conceptualization: SeS, DA, and JB; methodology: SeS; software: SeS and AA; validation: SeS; formal analysis: SeS; investigation: SeS; computational resources: JB; data curation: SeS; writing – original draft preparation: SeS; writing – review and editing: JB and DA; visualization: SeS; supervision: JB and DA; project administration: JB; funding acquisition: JB. All authors have read and agreed to the published version of the article.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Computational resources provided by the Digital Research Alliance of Canada and Advanced Research Computing at the University of British Columbia are gratefully acknowledged.

Financial support

This research has been supported by the Natural Sciences and Engineering Research Council of Canada (grant no. 556326).

Review statement

This paper was edited by Sukanta Basu and reviewed by Javier Sanz Rodrigo and one anonymous referee.


Ainslie, J.: Calculating the flowfield in the wake of wind turbines, J. o Wind Eng. Indust. Aerodynam., 27, 213–224,, 1988. a

Allaerts, D. and Meyers, J.: Boundary-layer development and gravity waves in conventionally neutral wind farms, J. Fluid Mech., 814, 95–130,, 2017. a, b

Allaerts, D. and Meyers, J.: Gravity Waves and Wind-Farm Efficiency in Neutral and Stable Conditions, Bound.-Lay. Meteorol., 166, 269–299,, 2018.  a

Allaerts, D. and Meyers, J.: Sensitivity and feedback of wind-farm-induced gravity waves, J. Fluid Mech., 862, 990–1028,, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x

Allaerts, D., Broucke, S. V., van Lipzig, N., and Meyers, J.: Annual impact of wind-farm gravity waves on the Belgian-Dutch offshore wind-farm cluster, J. Phys.: Conf. Ser., 1037, 072006,, 2018. a

Bastankhah, M. and Porté-Agel, F.: A new analytical model for wind-turbine wakes, Renew. Energy, 70, 116–123,, 2014. a, b, c, d, e, f, g, h, i

Bleeg, J., Purcell, M., Ruisi, R., and Traiger, E.: Wind Farm Blockage and the Consequences of Neglecting Its Impact on Energy Production, Energies, 11, 1609,, 2018. a, b

Blondel, F. and Cathelain, M.: An alternative form of the super-Gaussian wind turbine wake model, Wind Energy Science, 5, 1225–1236,, 2020. a, b

Branlard, E. and Gaunaa, M.: Cylindrical vortex wake model: right cylinder, Wind Energy, 18, 1973–1987,, 2014. a, b, c, d

Branlard, E. and Meyer Forsting, A. R.: Assessing the blockage effect of wind turbines and wind farms using an analytical vortex model, Wind Energy, 23, 2068–2086,, 2020. a, b, c, d, e

Branlard, E., Quon, E., Forsting, A. R. M., King, J., and Moriarty, P.: Wind farm blockage effects: comparison of different engineering models, J. Phys.: Conf. Ser., 1618, 062036,, 2020. a, b

Brogna, R., Feng, J., Sørensen, J. N., Shen, W. Z., and Porté-Agel, F.: A new wake model and comparison of eight algorithms for layout optimization of wind farms in complex terrain, Appl. Energy, 259, 114189,, 2020. a

Businger, J. A.: Transfer of Momentum and Heat in the Planetary Boundary Layer, in: Proceedings of the Symposium on Arctic Heat Budget and Atmospheric Circulation, Lake Arrowhead, The Rand Corporation, 305–332, (last access: 3 May 2024), 1966. a

Calaf, M., Meneveau, C., and Meyers, J.: Large eddy simulations of fully developed wind-turbine array boundary layers, Phys. Fluids, 22, 015110,, 2010. a, b

Centurelli, G., Vollmer, L., Schmidt, J., Dörenkämper, M., Schröder, M., Lukassen, L. J., and Peinke, J.: Evaluating Global Blockage engineering parametrizations with LES, J. Phys.: Conf. Ser., 1934, 012021,, 2021. a

Crespo, A. and Hernandez, J.: Turbulence characteristics in wind-turbine wakes, J. Wind Eng. Indust. Aerodynam., 61, 71–85,, 1996. a

Devesse, K., Lanzilao, L., Jamaer, S., van Lipzig, N., and Meyers, J.: Including realistic upper atmospheres in a wind-farm gravity-wave model, Wind Energ. Sci., 7, 1367–1382,, 2022. a

Devesse, K., Lanzilao, L., and Meyers, J.: A meso-micro atmospheric perturbation model for wind farm blockage, arXiv [preprint],, 2023. a

Etling, D.: Modelling the vertical ABL structure, World Scientific, 45–86,, 1996. a

Farrell, A., King, J., Draxl, C., Mudafort, R., Hamilton, N., Bay, C. J., Fleming, P., and Simley, E.: Design and analysis of a wake model for spatially heterogeneous flow, Wind Energ. Sci., 6, 737–758,, 2021. a, b

Frandsen, S., Barthelmie, R., Pryor, S., Rathmann, O., Larsen, S., Højstrup, J., and Thøgersen, M.: Analytical modelling of wind speed deficit in large offshore wind farms, Wind Energy, 9, 39–53,, 2006. a

Gill, A. E.: Chapter Four – Equations Satisfied by a Moving Fluid, in: Atmosphere–Ocean Dynamics, vol. 30 of International Geophysics, Academic Press, 63–94,, 1982. a

Gribben, B. J. and Hawkes, G. S.: A potential flow model for wind turbine induction and wind farm blockage, Systems and Engineering Technology, (last access: 3 May 2024), 2019. a, b

Jensen, N.: A note on wind generator interaction, no. 2411 in Risø-M, Risø National Laboratory, ISBN 87-550-0971-9, 1983. a

Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5MW Reference Wind Turbine for Offshore System Development, NREL – National Renewable Energy Laboratory,, 2009. a, b

Katić, I., Højstrup, J., and Jensen, N.: A Simple Model for Cluster Efficiency, in: European Wind Energy Association Conference and Exhibition, EWEC '86, edited by: Palz, W. and Sesto, E., A. Raguzzi, 407–410, (last access: 3 May 2024), 1987. a, b

Lanzilao, L. and Meyers, J.: A new wake-merging method for wind-farm power prediction in the presence of heterogeneous background velocity fields, Wind Energy, 25, 237–259,, 2022. a, b

Lanzilao, L. and Meyers, J.: An Improved Fringe-Region Technique for the Representation of Gravity Waves in Large Eddy Simulation with Application to Wind Farms, Bound.-Lay. Meteorol., 183, 567–593,, 2023. a, b, c

Lanzilao, L. and Meyers, J.: A parametric large-eddy simulation study of wind-farm blockage and gravity waves in conventionally neutral boundary layers, J. Fluid Mech., 979, A54,, 2024. a

Larsen, G.: A simple wake calculation procedure, Tech. Rep. Risø-M-2760 Risø, Tech. rep., Risø DTU National Laboratory for Sustainable Energy, (last access: 3 May 2024), 1988. a

Lin, Y.-L.: Mesoscale Dynamics, Cambridge University Press, ISBN 9780511619649,, 2007. a, b

Lissaman, P. B. S.: Energy Effectiveness of Arbitrary Arrays of Wind Turbines, J. Energy, 3, 323–328,, 1979. a

Meneveau, C.: The top-down model of wind farm boundary layers and its applications, J. Turbul., 13, N7,, 2011. a, b

Meyers, J. and Meneveau, C.: Large Eddy Simulations of Large Wind-Turbine Arrays in the Atmospheric Boundary Layer, ARC, ISBN 978-1-60086-959-4,, 2010. a

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

Nappo, C. J. (Ed.): Copyright, in: vol. 102 ofInternational Geophysics, Academic Press,, 2012. a

Niayifar, A. and Porté-Agel, F.: Analytical Modeling of Wind Farms: A New Approach for Power Prediction, Energies, 9, 741,, 2016. a, b, c, d, e, f, g, h, i, j

Nieuwstadt, F. T. M.: On the solution of the stationary, baroclinic Ekman-layer equations with a finite boundary-layer height, Bound.-Lay. Meteorol., 26, 377–390,, 1983. a, b

Nishino, T. and Dunstan, T. D.: Two-scale momentum theory for time-dependent modelling of large wind farms, J. Fluid Mech., 894, A2,, 2020. a

Nygaard, N. G., Steen, S., Poulsen, L., and Pedersen, J. G.: Modelling cluster wakes and wind farm blockage, J. Phys.: Conf. Ser., 1618, 062072,, 2020. a

Nygaard, N. G., Poulsen, L., Svensson, E., and Grønnegaard Pedersen, J.: Large-scale benchmarking of wake models for offshore wind farms, J. Phys.: Conf. Ser., 2265, 022008,, 2022. a

OpenWind: Theoretical basis and validation, AWS Truepower, Version 1.3., Tech. rep., UL Renewables, (last access: 3 May 2024), 2010. a

Ørsted: Ørsted presents update on its long-term financial targets, (last access: 3 May 2024), 2019. a

Panofsky, H. A.: Determination of stress from wind and temperature measurements, Q. J. Roy. Meteorol. Soc. 89, 85–94,, 1963. a

Rampanelli, G. and Zardi, D.: A Method to Determine the Capping Inversion of the Convective Boundary Layer, J. Appl. Meteorol., 43, 925–933,<0925:AMTDTC>2.0.CO;2, 2004. a, b, c

Saad, Y. and Schultz, M. H.: GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM J. Scient. Stat. Comput., 7, 856–869,, 1986. a

Segalini, A.: An analytical model of wind-farm blockage, J. Renew. Sustain. Energ., 13, 033307,, 2021. a

Shapiro, C., Starke, G., Meneveau, C., and Gayme, D.: A Wake Modeling Paradigm for Wind Farm Design and Control, Energies, 12, 2956,, 2019. a, b, c

Smith, R. B.: Linear theory of stratified hydrostatic flow past an isolated mountain, Tellus, 32, 348–364,, 1980. a, b

Smith, R. B.: Interacting Mountain Waves and Boundary Layers, J. Atmos. Sci., 64, 594–607,, 2007. a

Smith, R. B.: Gravity wave effects on wind farm efficiency, Wind Energy, 13, 449–458,, 2010. a, b, c, d, e

Smith, R. B.: The wind farm pressure field, Wind Energ. Sci., 9, 253–261,, 2024. a

Starke, G. M., Meneveau, C., King, J. R., and Gayme, D. F.: The area localized coupled model for analytical mean flow prediction in arbitrary wind farm geometries, J. Renew. Sustain. Energ., 13, 033305,, 2021. a

Stevens, R., Gayme, D., and Meneveau, C.: Coupled wake boundary layer model of wind-farms, J. Renew. Sustain. Energ., 7, 023115,, 2015. a, b

Stevens, R. J. A. M.: Dependence of optimal wind turbine spacing on wind farm length, Wind Energy, 19, 651–663,, 2016. a

Stipa, S., Ajay, A., and Brinkerhoff, J.: Toolbox fOr Stratified Convective Atmospheres (TOSCA), OSF [code],, 2023.  a

Stipa, S., Ajay, A., Allaerts, D., and Brinkerhoff, J.: TOSCA – an open-source, finite-volume, large-eddy simulation (LES) environment for wind farm flows, Wind Energ. Sci., 9, 297–320,, 2024. a, b, c, d

Troldborg, N. and Meyer Forsting, A.: A simple model of the wind turbine induction zone derived from numerical simulations, Wind Energy, 20, 2011–2020,, 2017. a

Voutsinas, S., Rados, K., and Zervos, A.: On the Analysis of Wake Effects in Wind Parks, Wind Eng., 14, 204–219, 1990. a

Wu, K. L. and Porté-Agel, F.: Flow Adjustment Inside and Around Large Finite-Size Wind Farms, Energies, 10, 2164,, 2017. a

Short summary
This paper introduces the multi-scale coupled (MSC) model, an engineering framework aimed at modeling turbine–wake and wind farm–gravity wave interactions, as well as local and global blockage effects. Comparisons against large eddy simulations show that the MSC model offers a valid contribution towards advancing our understanding of the coupled wind farm–atmosphere interaction, helping refining power estimation methodologies for existing and future wind farm sites.
Final-revised paper