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

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 (Meneveau, 2011). 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 Gaunaa, 2014; Branlard and Meyer Forsting, 2020), 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 Forsting, 2020). 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é-Agel, 2017; Allaerts and Meyers, 2017, 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; Ainslie, 1988; Larsen, 1988; OpenWind, 2010; Bastankhah and Porté-Agel, 2014; Niayifar and Porté-Agel, 2016). 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 Meneveau, 2011, and Stevens, 2016) 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=U_{\mathrm{b}}/\sqrt{g^{\prime }H}$, where U_b is the bulk velocity inside the boundary layer, g^′ is the reduced gravity, and H is the inversion height (Smith, 2010). 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 1 Conceptual 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

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{U}_{\mathrm{1}}\frac{\partial u_{\mathrm{1}}}{\partial x}+V_{\mathrm{1}}\frac{\partial u_{\mathrm{1}}}{\partial y}+\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial x}=-f_{\mathrm{c}}{v}_{\mathrm{1}}+{\mathit{\nu }}_{t,\mathrm{1}}\left(\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{1}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{1}}}{\partial y^{\mathrm{2}}}\right)\\ +\frac{D_{\mathrm{11}}}{H_{\mathrm{1}}}\left(u_{\mathrm{2}}-u_{\mathrm{1}}\right)+\frac{D_{\mathrm{12}}}{H_{\mathrm{1}}}\left(v_{\mathrm{2}}-v_{\mathrm{1}}\right)-\frac{C_{\mathrm{11}}}{H_{\mathrm{1}}}{u}_{\mathrm{1}}-\frac{C_{\mathrm{12}}}{H_{\mathrm{1}}}{v}_{\mathrm{1}}\\ -\frac{f_x}{H_{\mathrm{1}}}\\ U_{\mathrm{1}}\frac{\partial v_{\mathrm{1}}}{\partial x}+V_{\mathrm{1}}\frac{\partial v_{\mathrm{1}}}{\partial y}+\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial y}=f_{\mathrm{c}}{u}_{\mathrm{1}}+{\mathit{\nu }}_{t,\mathrm{1}}\left(\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{1}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{1}}}{\partial y^{\mathrm{2}}}\right)\\ +\frac{D_{\mathrm{21}}}{H_{\mathrm{1}}}\left(u_{\mathrm{2}}-u_{\mathrm{1}}\right)+\frac{D_{\mathrm{22}}}{H_{\mathrm{1}}}\left(v_{\mathrm{2}}-v_{\mathrm{1}}\right)-\frac{C_{\mathrm{21}}}{H_{\mathrm{1}}}{u}_{\mathrm{1}}-\frac{C_{\mathrm{22}}}{H_{\mathrm{1}}}{v}_{\mathrm{1}}\\ -\frac{f_y}{H_{\mathrm{1}}}\\ U_{\mathrm{1}}\frac{\partial {\mathit{\eta }}_{\mathrm{1}}}{\partial x}+V_{\mathrm{1}}\frac{\partial {\mathit{\eta }}_{\mathrm{1}}}{\partial y}+H_{\mathrm{1}}\left(\frac{\partial u_{\mathrm{1}}}{\partial x}+\frac{\partial v_{\mathrm{1}}}{\partial y}\right)=\mathrm{0}\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}{U}_{\mathrm{2}}\frac{\partial u_{\mathrm{2}}}{\partial x}+V_{\mathrm{2}}\frac{\partial u_{\mathrm{2}}}{\partial y}+\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial x}=-f_{\mathrm{c}}{v}_{\mathrm{2}}+{\mathit{\nu }}_{t,\mathrm{2}}\left(\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{2}}}{\partial y^{\mathrm{2}}}\right)\\ -\frac{D_{\mathrm{11}}}{H_{\mathrm{2}}}\left(u_{\mathrm{2}}-u_{\mathrm{1}}\right)-\frac{D_{\mathrm{12}}}{H_{\mathrm{2}}}\left(v_{\mathrm{2}}-v_{\mathrm{1}}\right)\\ U_{\mathrm{2}}\frac{\partial v_{\mathrm{2}}}{\partial x}+V_{\mathrm{2}}\frac{\partial v_{\mathrm{2}}}{\partial y}+\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial y}=f_{\mathrm{c}}{u}_{\mathrm{2}}+{\mathit{\nu }}_{t,\mathrm{2}}\left(\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{2}}}{\partial y^{\mathrm{2}}}\right)\\ -\frac{D_{\mathrm{21}}}{H_{\mathrm{2}}}\left(u_{\mathrm{2}}-u_{\mathrm{1}}\right)-\frac{D_{\mathrm{22}}}{H_{\mathrm{2}}}\left(v_{\mathrm{2}}-v_{\mathrm{1}}\right)\\ U_{\mathrm{2}}\frac{\partial {\mathit{\eta }}_{\mathrm{2}}}{\partial x}+V_{\mathrm{2}}\frac{\partial {\mathit{\eta }}_{\mathrm{2}}}{\partial y}+H_{\mathrm{2}}\left(\frac{\partial u_{\mathrm{2}}}{\partial x}+\frac{\partial v_{\mathrm{2}}}{\partial y}\right)=\mathrm{0},\end{array}\right.\end{array}$$

where subscripts 1 and 2 indicate the wind farm and upper layer respectively, η is the layer displacement, H is the layer height, $\mathbf{U}=(U,V)$ and $\mathit{u}=(u,v)$ are the depth-averaged background and perturbation velocities respectively, ν_t is the eddy viscosity (see Sect. 2.3), and f_c=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 (Gill, 1982). The term $\mathit{f}/H_{\mathrm{1}}$ 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 C_ij and D_ij describe the perturbation of friction at the ground and at the top of the wind farm layer respectively, defined as

$$\begin{array}{}\text{(3)}& C_{ij}={\displaystyle \frac{∥\mathit{\tau }{|}_{z=\mathrm{0}}∥}{{∥{\mathbf{U}}_{\mathrm{1}}∥}^{\mathrm{3}}}}\left[\begin{array}{cc}\mathrm{2}{U}_{\mathrm{1}}^{\mathrm{2}}+V_{\mathrm{1}}^{\mathrm{2}}& U_{\mathrm{1}}{V}_{\mathrm{1}}\\ V_{\mathrm{1}}{U}_{\mathrm{1}}& U_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{V}_{\mathrm{1}}^{\mathrm{2}}\end{array}\right]\end{array}$$

and

$$\begin{array}{}\text{(4)}& D_{ij}={\displaystyle \frac{∥\mathit{\tau }{|}_{z=H_{\mathrm{1}}}∥}{{∥\mathrm{\Delta }\mathbf{U}∥}^{\mathrm{3}}}}\left[\begin{array}{cc}\mathrm{2}\mathrm{\Delta }{U}^{\mathrm{2}}+\mathrm{\Delta }{V}^{\mathrm{2}}& \mathrm{\Delta }U\mathrm{\Delta }V\\ \mathrm{\Delta }V\mathrm{\Delta }U& \mathrm{\Delta }{U}^{\mathrm{2}}+\mathrm{2}\mathrm{\Delta }{V}^{\mathrm{2}}\end{array}\right],\end{array}$$

where $\mathrm{\Delta }\mathbf{U}=(\mathrm{\Delta }U,\mathrm{\Delta }V)=(U_{\mathrm{2}}-U_{\mathrm{1}},V_{\mathrm{2}}-V_{\mathrm{1}})$ 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

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\widehat{p}=\left[{\displaystyle \frac{i\left(N^{\mathrm{2}}-{\mathrm{\Omega }}^{\mathrm{2}}\right)}{m}}+g^{\prime }\right]\widehat{\mathit{\eta }}\end{array}$$

based on the linear, three-dimensional, non-hydrostatic, non-rotating gravity wave theory of Smith (1980) (see also Lin, 2007, and Nappo, 2012). The top-hat symbol is used to indicate that Eq. (5) is expressed in Fourier space, $\mathit{\eta }={\mathit{\eta }}_{\mathrm{1}}+{\mathit{\eta }}_{\mathrm{2}}$ is the total vertical ABL displacement, $N=\sqrt{\mathit{\gamma }g/{\mathit{\theta }}_{\mathrm{0}}}$ is the Brunt–Väisälä frequency, where the free-atmosphere lapse rate $\mathit{\gamma }=\mathrm{d}\mathit{\theta }/\mathrm{d}z$ is assumed constant in this study, $g^{\prime }=g\mathrm{\Delta }\mathit{\theta }/{\mathit{\theta }}_{\mathrm{0}}$ is referred to as the reduced gravity, where Δθ is the potential temperature jump across the inversion layer, and θ₀ is the reference potential temperature. The quantity m identifies the vertical wavenumber, evaluated using the dispersion relation

$$\begin{array}{}\text{(6)}& m^{\mathrm{2}}=\left(k^{\mathrm{2}}+l^{\mathrm{2}}\right)\left({\displaystyle \frac{N^{\mathrm{2}}}{{\mathrm{\Omega }}^{\mathrm{2}}}}-\mathrm{1}\right),\end{array}$$

in which $\mathit{\kappa }=(k,l)$ is the horizontal wavenumber and $\mathrm{\Omega }=-\mathit{\kappa }\cdot {\mathbf{U}}_{\mathrm{3}}$ 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 e_k 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

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{U_{\mathrm{\infty }}-U\left(\mathit{x}\right)}{U_{\mathrm{\infty }}}}=W\left(\mathit{x}\right)=C\left(x\right)f\left(r\right(\mathit{x}),\mathit{\sigma }(x\left)\right)\mathcal{H}\left(x\right),\end{array}$$

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 Meyers, 2022, 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

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }f(r,\mathit{\sigma })=\exp \left(-{\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right)\text{(9)}& {\displaystyle }& {\displaystyle }C\left(x\right)=\mathrm{1}-\sqrt{\mathrm{1}-{\displaystyle \frac{C_{\mathrm{T}}}{\mathrm{8}(\mathit{\sigma }/D{)}^{\mathrm{2}}}}},\end{array}$$

where C_T 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 $C_{\mathrm{T}}>\mathrm{8}{\mathit{\sigma }}^{\mathrm{2}}/D^{\mathrm{2}}$. Specifically, depending on the value of the C_T 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 $C_{\mathrm{T}}<\mathrm{8}{\mathit{\sigma }}^{\mathrm{2}}/D^{\mathrm{2}}$. 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 C_T 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

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\mathit{\sigma }\left(x\right)=k^{*}x+c_{\mathrm{s}}\sqrt{\mathit{\beta }}D\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T}}}}{\sqrt{\mathrm{1}-C_{\mathrm{T}}}}}\text{(12)}& {\displaystyle }& {\displaystyle }{k}^{*}=a_{\mathrm{s}}\mathrm{TI}+b_{\mathrm{s}},\end{array}$$

where, according to Niayifar and Porté-Agel (2016), a_s=0.3837, b_s=0.003678, and c_s=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 C_T, 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

$${\mathrm{TI}}_k=\sqrt{{\mathrm{TI}}_{\mathrm{\infty }}^{\mathrm{2}}+{\mathrm{TI}}_{+,k}^{\mathrm{2}}},$$

where TI_∞ is the background turbulence intensity provided as a user-defined input, while TI${}_{+,k}$ is the added turbulence intensity, evaluated as

$$\begin{array}{}\text{(13)}& {\mathrm{TI}}_{+,k}=\underset{j}{\max }\left(\mathcal{H}\left(x_k{|}_j\right)I_j\left(x_k{|}_j\right){\displaystyle \frac{A_{\mathrm{w}}}{\mathit{\pi }{R}_k^{\mathrm{2}}}}\right),\end{array}$$

where ℋ is again the Heaviside function, x_k|_j is the distance between turbine k and turbine j rotor center, in the wake reference frame of turbine j, R_k is rotor radius of turbine k, and A_w is an area whose value goes from 0 to $\mathit{\pi }{R}_k^{\mathrm{2}}$, depending on whether turbine j is waking k or not, with 0 corresponding to the non-waking situation (see Niayifar and Porté-Agel, 2016). 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 I_j is an empirical function which, according to Crespo and Hernandez (1996), is given by

$$\begin{array}{}\text{(14)}& I_j\left(x\right)=d_{\mathrm{s}}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T},j}}\right)\right)}^{e_{\mathrm{s}}}{\mathrm{TI}}_{\mathrm{\infty }}^{f_{\mathrm{s}}}{\left({\displaystyle \frac{x}{D_j}}\right)}^{g_{\mathrm{s}}}.\end{array}$$

In the present study, we use the values proposed by Niayifar and Porté-Agel (2016) for e_s, f_s, and g_s, equal to 0.8325, 0.0325, and −0.32 respectively, while we vary d_s in order to minimize $∥{\mathrm{TI}}_{k,\mathrm{LES}}-{\mathrm{TI}}_{k,\mathrm{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 d_s=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:

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }U\left(\mathit{x}\right)=U_{\mathrm{\infty }}-\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{u}_k{W}_k\left(\mathit{x}\right)\text{(16)}& {\displaystyle }& {\displaystyle }U\left(\mathit{x}\right)=U_{\mathrm{\infty }}-\sqrt{\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{\left(u_k{W}_k\left(\mathit{x}\right)\right)}^{\mathrm{2}}},\end{array}$$

where k denotes the kth wind turbine, N_t is the number of turbines in the array, u_k 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 e_k 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 u_k=U_∞, leading to a simple summation of W(x) or W(x)². This produces deficits that are too high and has been addressed by Niayifar and Porté-Agel (2016) and Voutsinas et al. (1990), who used u_k=U(x_k) with x_k, 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 N_q=16 quadrature points, evaluating u_k as

$$\begin{array}{}\text{(17)}& u_k={\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}U\left({\mathit{x}}_{k,q}\right),\end{array}$$

where x_k,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 2 Comparison 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. f_k=f_k(U_∞). In the coupling method proposed by Allaerts and Meyers (2019), f_k is corrected for blockage effects based on a perturbation velocity u^up so that f_k=f_k(U), where $U=U_{\mathrm{\infty }}+u^{\mathrm{up}}$. The point-wise perturbation velocity u^up is calculated by evaluating the meso-scale perturbation velocity u₁(x) at x^up, 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 $x^{\mathrm{up}}=\min \left(x_k\right)-\mathrm{10}$ D and $y^{\mathrm{up}}=\mathrm{1}/N_{\mathrm{t}}\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{y}_k$, where x_k and y_k are the rotor center coordinates of turbine k. Finally, f_k 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

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathit{f}(x,y)=& \underset{\mathrm{0}}{\overset{L_x}{\int }}\underset{\mathrm{0}}{\overset{L_y}{\int }}G(x-x^{\prime },y-y^{\prime })\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{\mathit{f}}_k\\ & \mathit{\delta }\left(x^{\prime }-x_k,y^{\prime }-y_k\right)\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime },\end{array}\end{array}$$

with

$$\begin{array}{}\text{(19)}& G(x,y)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }{L}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{x^{\mathrm{2}}+y^{\mathrm{2}}}{L^{\mathrm{2}}}}\right),\end{array}$$

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 f_k's dependency on u^up by expanding f_k(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 u^up 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 u_k, which is calculated at each turbine location using Eqs. (7) and (15). At this point, turbine thrust is evaluated as

$$\begin{array}{}\text{(20)}& {\mathit{f}}_k={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{u}_k^{\mathrm{2}}{\displaystyle \frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}}}{C}_{\mathrm{T},k}{\mathit{e}}_k,\end{array}$$

where e_k is defined in Sect. 2.1.2, and C_T,k is the thrust coefficient of the kth turbine. After filtering f_k 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 $U^{i+\mathrm{1}}=U_{\mathrm{\infty }}+u^{\mathrm{up}}$. The wake model is run again with the updated free-stream wind by setting $U_{\mathrm{\infty }}=U^{i+\mathrm{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={∥p^{i+\mathrm{1}}-p^i∥}_{\mathrm{2}}/{∥p^{i+\mathrm{1}}∥}_{\mathrm{2}}$, 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 3 Conceptual 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/C_T 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

$$\begin{array}{}\text{(21)}& {\displaystyle }\left\{\begin{array}{l}{U}_{\mathrm{1}}\frac{\partial u_{\mathrm{1}}^{\mathrm{b}k}}{\partial x}+V_{\mathrm{1}}\frac{\partial u_{\mathrm{1}}^{\mathrm{b}k}}{\partial y}=-f_{\mathrm{c}}{v}_{\mathrm{1}}^{\mathrm{b}k}+{\mathit{\nu }}_{t,\mathrm{1}}\left(\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{b}k}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{b}k}}{\partial y^{\mathrm{2}}}\right)\\ +\frac{D_{\mathrm{11}}}{H_{\mathrm{1}}}\left(u_{\mathrm{2}}^{\mathrm{b}k}-u_{\mathrm{1}}^{\mathrm{b}k}\right)+\frac{D_{\mathrm{12}}}{H_{\mathrm{1}}}\left(v_{\mathrm{2}}^{\mathrm{b}k}-v_{\mathrm{1}}^{\mathrm{b}k}\right)-\frac{C_{\mathrm{11}}}{H_{\mathrm{1}}}{u}_{\mathrm{1}}^{\mathrm{b}k}\\ -\frac{C_{\mathrm{12}}}{H_{\mathrm{1}}}{v}_{\mathrm{1}}^{\mathrm{b}k}-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial x}\\ U_{\mathrm{1}}\frac{\partial v_{\mathrm{1}}^{\mathrm{b}k}}{\partial x}+V_{\mathrm{1}}\frac{\partial v_{\mathrm{1}}^{\mathrm{b}k}}{\partial y}=f_{\mathrm{c}}{u}_{\mathrm{1}}^{\mathrm{b}k}+{\mathit{\nu }}_{t,\mathrm{1}}\left(\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{1}}^{\mathrm{b}k}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{1}}^{\mathrm{b}k}}{\partial y^{\mathrm{2}}}\right)\\ +\frac{D_{\mathrm{21}}}{H_{\mathrm{1}}}\left(u_{\mathrm{2}}^{\mathrm{b}k}-u_{\mathrm{1}}^{\mathrm{b}k}\right)+\frac{D_{\mathrm{22}}}{H_{\mathrm{1}}}\left(v_{\mathrm{2}}^{\mathrm{b}k}-v_{\mathrm{1}}^{\mathrm{b}k}\right)-\frac{C_{\mathrm{21}}}{H_{\mathrm{1}}}{u}_{\mathrm{1}}^{\mathrm{b}k}\\ -\frac{C_{\mathrm{22}}}{H_{\mathrm{1}}}{v}_{\mathrm{1}}^{\mathrm{b}k}-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial y}\end{array}\right.\text{(22)}& {\displaystyle }\left\{\begin{array}{l}{U}_{\mathrm{2}}\frac{\partial u_{\mathrm{2}}^{\mathrm{b}k}}{\partial x}+V_{\mathrm{2}}\frac{\partial u_{\mathrm{2}}^{\mathrm{b}k}}{\partial y}=-f_{\mathrm{c}}{v}_{\mathrm{2}}^{\mathrm{b}k}+{\mathit{\nu }}_{t,\mathrm{2}}\left(\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{b}k}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{b}k}}{\partial y^{\mathrm{2}}}\right)\\ -\frac{D_{\mathrm{11}}}{H_{\mathrm{2}}}\left(u_{\mathrm{2}}^{\mathrm{b}k}-u_{\mathrm{1}}^{\mathrm{b}k}\right)-\frac{D_{\mathrm{12}}}{H_{\mathrm{2}}}\left(v_{\mathrm{2}}^{\mathrm{b}k}-v_{\mathrm{1}}^{\mathrm{b}k}\right)-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial x}\\ U_{\mathrm{2}}\frac{\partial v_{\mathrm{2}}^{\mathrm{b}k}}{\partial x}+V_{\mathrm{2}}\frac{\partial v_{\mathrm{2}}^{\mathrm{b}k}}{\partial y}=f_{\mathrm{c}}{u}_{\mathrm{2}}^{\mathrm{b}k}+{\mathit{\nu }}_{t,\mathrm{2}}\left(\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{2}}^{\mathrm{b}k}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{v}_{\mathrm{2}}^{\mathrm{b}k}}{\partial y^{\mathrm{2}}}\right)\\ -\frac{D_{\mathrm{21}}}{H_{\mathrm{2}}}\left(u_{\mathrm{2}}^{\mathrm{b}k}-u_{\mathrm{1}}^{\mathrm{b}k}\right)-\frac{D_{\mathrm{22}}}{H_{\mathrm{2}}}\left(v_{\mathrm{2}}^{\mathrm{b}k}-v_{\mathrm{1}}^{\mathrm{b}k}\right)-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial y}\end{array}\right.,\end{array}$$

where subscripts 1 and 2 again refer to the wind farm and upper layers, respectively, and ${\mathit{u}}^{\mathrm{b}k}=(u^{\mathrm{b}k},v^{\mathrm{b}k})$ 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, u^bk 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 U_b(x) with the layer average ${\mathbf{U}}_{\mathrm{1}}+{\mathit{u}}_{\mathrm{1}}^{\mathrm{b}k}$ 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

$$\begin{array}{}\text{(23)}& ∥{\mathbf{U}}_{\mathrm{b}}\left(z\right)∥={\displaystyle \frac{u^{*}}{\mathit{\kappa }}}\left(\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-\mathrm{\Psi }\left({\displaystyle \frac{z}{L}}\right)\right),\end{array}$$

where u^* is the friction velocity, z₀ 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

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\sqrt{{\left(U_{\mathrm{1}}+u_{\mathrm{1}}^{\mathrm{b}k}\right)}^{\mathrm{2}}+{\left(V_{\mathrm{1}}+v_{\mathrm{1}}^{\mathrm{b}k}\right)}^{\mathrm{2}}}=& {\displaystyle \frac{u^{*}(x,y)}{\mathit{\kappa }\left(H_{\mathrm{1}}-z_{\mathrm{0}}\right)}}\\ & \underset{z_{\mathrm{0}}}{\overset{H_{\mathrm{1}}}{\int }}\left(\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)\right)\mathrm{d}z,\end{array}\end{array}$$

which yields

$$\begin{array}{}\text{(25)}& u^{*}(x,y)={\displaystyle \frac{\mathit{\kappa }\left(H_{\mathrm{1}}-z_{\mathrm{0}}\right)\sqrt{{\left(U_{\mathrm{1}}+u_{\mathrm{1}}^{\mathrm{b}k}\right)}^{\mathrm{2}}+{\left(V_{\mathrm{1}}+v_{\mathrm{1}}^{\mathrm{b}k}\right)}^{\mathrm{2}}}}{H_{\mathrm{1}}\left[\ln \left(\frac{H_{\mathrm{1}}}{z_{\mathrm{0}}}\right)-\mathrm{1}\right]+z_{\mathrm{0}}}}.\end{array}$$

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

$$\begin{array}{}\text{(26)}& {\mathit{\varphi }}^{\prime }(x,y)=\arctan \left({\displaystyle \frac{V_{\mathrm{1}}+v_{\mathrm{1}}^{\mathrm{b}k}}{U_{\mathrm{1}}+u_{\mathrm{1}}^{\mathrm{b}k}}}\right)-\arctan \left({\displaystyle \frac{V_{\mathrm{1}}}{U_{\mathrm{1}}}}\right)\end{array}$$

so that the overall flow angle at a given location is $\mathit{\varphi }\left(\mathit{x}\right)=\mathrm{\Phi }\left(z\right)+{\mathit{\varphi }}^{\prime }(x,y)$. The background velocity profile inside the wind farm layer can finally be expressed as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathbf{U}}_{\mathrm{b}}\left(\mathit{x}\right)& ={\displaystyle \frac{u^{*}(x,y)}{\mathit{\kappa }}}\left(\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)\cos \right(\mathit{\varphi }\left(\mathit{x}\right))\mathit{i}\\ & +\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)\sin \left(\mathit{\varphi }\right(\mathit{x}\left)\right)\mathit{j}).\end{array}\end{array}$$

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 U_b(x) from ${\mathit{u}}_{\mathrm{1}}^{\mathrm{b}k}$.

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 U_b(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_{\mathrm{\infty }}=u_k=∥{\mathbf{U}}_{\mathrm{b}}\left(\mathit{x}\right)∥$. 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\left(\mathit{x}{|}_k\right)$, 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_∞=∥U_b(x)∥. With the above variations, Eqs. (15 and (16) become

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }{U}_{\mathrm{w}}\left(\mathit{x}\right)=∥{\mathbf{U}}_{\mathrm{b}}\left(\mathit{x}\right)∥-\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{u}_{k}W\left(\mathit{x}{|}_k\right)\text{(29)}& {\displaystyle }& {\displaystyle }{U}_{\mathrm{w}}\left(\mathit{x}\right)=∥{\mathbf{U}}_{\mathrm{b}}\left(\mathit{x}\right)∥-\sqrt{\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{\left(u_k{W}_k\left(\mathit{x}{|}_k\right)\right)}^{\mathrm{2}}},\end{array}$$

where U_w(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 u_k is equal to the average of U_w(x) among quadrature points of turbine k, rather then the free-stream velocity at location x. Moreover, it is worth noticing that U_w(x) is, in reality, a velocity magnitude, as the directional information contained in U_b(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 U_b(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 Gaunaa, 2014; Branlard and Meyer Forsting, 2020). 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

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{u}_k^{vc}(r,x^{\prime })& ={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{t}}}{\mathrm{2}\mathit{\pi }}}[{\displaystyle \frac{R-r+|R-r|}{\mathrm{2}|R-r|}}+{\displaystyle \frac{x^{\prime }k(r,x^{\prime })}{\mathrm{2}\mathit{\pi }\sqrt{rR}}}(K\left(k^{\mathrm{2}}(r,x^{\prime })\right)\\ & +{\displaystyle \frac{R-r}{R+r}}\mathrm{\Pi }\left(k^{\mathrm{2}}(r,\mathrm{0}),k^{\mathrm{2}}(r,x^{\prime })\right)\left)\right],\end{array}\end{array}$$

where $x^{\prime }=x-x_k$ is the axial distance from the rotor and $r=\sqrt{(y-y_k{)}^{\mathrm{2}}+(z-z_k{)}^{\mathrm{2}}}$ is the radial distance from the wake axis (with x_k, y_k, and z_k 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

$$\begin{array}{}\text{(31)}& k(r,x^{\prime })={\displaystyle \frac{\mathrm{4}rR}{(R+r{)}^{\mathrm{2}}+{x^{\prime }}^{\mathrm{2}}}}.\end{array}$$

The momentum theory calibration relation (Branlard and Meyer Forsting, 2020) is used to express γ_t as

$$\begin{array}{}\text{(32)}& {\mathit{\gamma }}_{\mathrm{t}}=-∥{\mathbf{U}}_{\mathrm{b}}\left({\mathit{x}}_k\right)∥\left(\mathrm{1}-\sqrt{\mathrm{1}-C\mathrm{T}}\right).\end{array}$$

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

$$\begin{array}{}\text{(33)}& U\left(\mathit{x}\right)=U_{\mathrm{w}}\left(\mathit{x}\right)+\sum _{k=\mathrm{1}}^{N_{\mathrm{t}}}{u}_k^{vc}(r,x^{\prime }).\end{array}$$

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 C_T, and power coefficient C_P. 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, C_T, and C_P 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 U_disk for each turbine. Then, the hypothetical free-stream wind speed is calculated as $U_{\mathrm{\infty }}=U_{\mathrm{disk}}/(\mathrm{1}-a)$, where $a=(\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}})/\mathrm{2}$ is the axial induction factor. This free-stream wind speed U_∞ is then used to update C_T and C_P 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 C_T 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 u₁, u₂, η, 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 u^bk, and Eqs. (25), (26), and (27) are used to reconstruct the background velocity field U_b(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 U_w(x) at the quadrature points. In particular, combining Eqs. (17) and (28) yields

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}{U}_{\mathrm{w}}\left({\mathit{x}}_{i,q}\right)& +\sum _{\begin{array}{c}j\\ j\ne i\end{array}}^{N_{\mathrm{t}}}{\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}{U}_{\mathrm{w}}\left({\mathit{x}}_{j,q}\right)W\left({\mathit{x}}_{i,q}{|}_j\right)\\ & ={\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}∥{\mathbf{U}}_{\mathrm{b}}\left({\mathit{x}}_{i,q}{|}_j\right)∥,\end{array}\end{array}$$

which produces the linear system Au^W=b, where

$$\begin{array}{}\text{(35)}& {\displaystyle }& {\displaystyle }{u}_i^W={\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}{U}_{\mathrm{w}}\left({\mathit{x}}_{i,q}\right)\text{(36)}& {\displaystyle }& {\displaystyle }{A}_{ij}=\left\{\begin{array}{ll}\mathrm{1}& i=j\\ \sum _{\begin{array}{c}j\\ j\ne i\end{array}}^{N_{\mathrm{t}}}\frac{\mathrm{1}}{N_q}\sum _{q=\mathrm{1}}^{N_q}W\left({\mathit{x}}_{i,q}{|}_j\right)& i\ne j\end{array}\right.\text{(37)}& {\displaystyle }& {\displaystyle }{b}_i={\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}∥{\mathbf{U}}_{\mathrm{b}}\left({\mathit{x}}_{i,q}{|}_j\right)∥.\end{array}$$

In Eqs. (34), (35), (36), and  (37), the notation ${\mathit{x}}_{i,q}{|}_j$ denotes quadrature point q of turbine i, seen in the local wind coordinate system of turbine  j.

Figure 4 Sketch 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 $u_i^W$ without considering turbine self-induction. Using Eq. (30) and averaging on quadrature points yields

$$\begin{array}{}\text{(38)}& u_i^{W+I}=u_i^W+\sum _{\begin{array}{c}j\\ j\ne i\end{array}}^{N_{\mathrm{t}}}{\displaystyle \frac{\mathrm{1}}{N_q}}\sum _{q=\mathrm{1}}^{N_q}{u}_j^{vc}\left(r_{i,q},x_{i,q}{|}_j\right),\end{array}$$

where r_i,q and $x_{i,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 C_T coefficient of each turbine is updated using $u_i^{W+I}$ and the input C_T curves, while turbine thrust is evaluated with Eq. (20), where we set $u_k=u_k^{W+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⁻⁴.

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 1 Input parameters for the MSC model.

In the context of the present study, we chose H₁ as twice the average turbine hub height, while $H_{\mathrm{2}}=H-H_{\mathrm{1}}$. 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=H₁ 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 Nieuwstadt, 1983), or from numerical simulations. The same holds for background friction velocity and ambient turbulence intensity at the hub height, where the condition $u{*}^{\mathrm{2}}=∥\mathit{\tau }{|}_{z=\mathrm{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

$$\begin{array}{}\text{(39)}& {\mathit{\nu }}_t\left(z\right)=\left\{\begin{array}{ll}\mathit{\kappa }{u}^{*}z{\left(\mathrm{1}-\frac{z}{H}\right)}^{\mathrm{2}}& z<H\\ \mathrm{0}& z\mathit{>=}H\end{array}\right..\end{array}$$

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 Meyers, 2023) 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 ${\mathit{\lambda }}_z=\mathrm{2}\mathit{\pi }{U}_{\mathrm{g}}/N$, where N is the Brunt—Väisälä frequency and U_g 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⁻¹ and a fringe damping coefficient of ν_FR=0.03 s⁻¹ 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 2 Fringe and advection damping region information.

The size of the successor domain is 40 km × 21 km × 28 km in the streamwise, spanwise, and vertical direction respectively, discretized with $\mathrm{1554}\times \mathrm{1194}\times \mathrm{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 3 Mesh information for the wind farm cases.

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×10⁵ 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 4 ABL parameters used for the finite wind farm simulation presented in this section.

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 f_c 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 10⁵ s; then data are averaged for 2×10⁴ 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 h_ref, 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 θ₀, 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 5 ABL 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 q_min within the boundary layer, and geostrophic wind angle.

Figure 5 Velocity 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 6 Instantaneous velocity magnitude at the hub height, sampled at the end of wind farm cases N1 (a) and N2 (b).

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

Figure 8 Perturbation 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 Meyers, 2017).

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 Gaunaa, 2014) or the Rankine half body (Gribben and Hawkes, 2019) 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 6 3LM input parameters for cases N1 and N2.

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 9 Effect 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 (Smith, 2024). 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⁻⁴. 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 10 Depth-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 h_hub, while dotted lines indicate $h_{\mathrm{hub}}\pm D/\mathrm{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 $u_{\mathrm{1}}^{\mathrm{b}k}$ is then compared with vertical profiles of U_b(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 L_xsin (Φ(H)), where L_x 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 U_b, 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 11 Left: 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 U_b(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 12 Micro-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 C_T of 0.9).

Figure 13 Comparison 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 14 Comparison 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 L_x. This results in a fluctuation of thrust and power that is characterized by a wavelength roughly equal to L_x, 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 15 Comparison 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 7 Total wind farm power obtained from LESs and as predicted by reduced-order models.

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⁻¹ 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 $C_{\mathrm{T}}>\mathrm{8}{\mathit{\sigma }}^{\mathrm{2}}/D^{\mathrm{2}}$. Hence, based on the value of C_T, the singularity region may extend more or less behind the rotor, or be completely non-existent. Nevertheless, as values of C_T 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 C_T 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 Cathelain, 2020).

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

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }f(r,\mathit{\sigma })=\exp \left(-{\displaystyle \frac{r^n}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}{D}^{\mathrm{2}-n}\right)\text{(A2)}& {\displaystyle }& {\displaystyle }C\left(x\right)={\mathrm{2}}^{\mathrm{2}/n-\mathrm{1}}-\sqrt{{\mathrm{2}}^{\mathrm{4}/n-\mathrm{2}}-{\displaystyle \frac{n{C}_{\mathrm{T}}}{\mathrm{16}\mathrm{\Gamma }(\mathrm{2}/n)(\mathit{\sigma }/D{)}^{\mathrm{4}/n}}}},\end{array}$$

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:

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }n\left(x\right)=w_{\mathrm{n}}\left(x\right)n_{\mathrm{n}}\left(x\right)+w_{\mathrm{f}}\left(x\right)n_{\mathrm{f}}\left(x\right)\text{(A4)}& {\displaystyle }& {\displaystyle }{n}_{\mathrm{n}}\left(x\right)=\mathrm{2}\exp (-\mathrm{0.68}x/D)+\mathrm{2}\text{(A5)}& {\displaystyle }& {\displaystyle }{n}_{\mathrm{f}}\left(x\right)=\mathrm{2}\text{(A6)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{n}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}-\tanh \left(\mathrm{7}\left(x-x_{\mathrm{c}}\right)/\mathit{\delta }\right)\right]\text{(A7)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{f}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}+\tanh \left(\mathrm{7}\left(x-x_{\mathrm{c}}\right)/\mathit{\delta }\right)\right],\end{array}$$

where n_n and w_n are the exponent and weight associated with the near wake, while n_f and w_f pertain to the far wake. The parameters x_c and δ are set to 2 and 4 D respectively. Figure A1 shows the evolution of n_n, n_f, 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 x_c and δ.

Figure A1 (a) Streamwise evolution of overall, near-wake, and far-wake components (n, n_n, and n_f 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⁻¹ 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⁻¹ and a C_T coefficient of 0.9, while, in the right panel we show velocity evolution at the wake centerline for three different values of C_T, at the same wind speed. As can be noticed, the Gaussian model is singular below roughly 1.5 D, for C_T=0.9, whereas a null wind speed is reached slightly below 1 D with C_T=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 C_T 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 (C_T or C_P) are used for a given wind speed. Moreover, for the NREL 5 MW turbine, a reliable C_T 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 C_T vs. wind speed curve that is biased towards higher values.

These considerations motivated us to evaluate C_T and C_P 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⁻¹, 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 C_T and C_P curves.

Data from the right panel of Fig. B1 have been used to update C_T and C_P 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 η₁ and η₂ 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 A^3LM of size 5N_kN_l×5N_kN_l, where N_k and N_l are the wavenumber sizes in the streamwise and spanwise directions respectively. We use a sparse matrix format to express A^3LM, while the complex linear system is solved using the generalized minimum residual (GMRES) (Saad and Schultz, 1986) iterative method. The full linear system, corresponding to the 3LM solution, is given by

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& \left[\begin{array}{cccc}A\left(k_{\mathrm{1}},l_{\mathrm{1}}\right)& & & \\ & A\left(k_{\mathrm{1}},l_{\mathrm{2}}\right)& & \\ & & \mathrm{\ddots }& \\ & & & A\left(k_{N_k},l_{N_l}\right)\end{array}\right]\\ & \left[\begin{array}{c}s\left(k_{\mathrm{1}},l_{\mathrm{1}}\right)\\ s\left(k_{\mathrm{1}},l_{\mathrm{2}}\right)\\ \mathrm{⋮}\\ s\left(k_{N_k},l_{N_l}\right)\end{array}\right]=\left[\begin{array}{c}b\left(k_{\mathrm{1}},l_{\mathrm{1}}\right)\\ b\left(k_{\mathrm{1}},l_{\mathrm{2}}\right)\\ \mathrm{⋮}\\ b\left(k_{N_k},l_{N_l}\right)\end{array}\right],\end{array}\end{array}$$

where

$$\begin{array}{}\text{(C2)}& \begin{array}{rl}& A\left(k_i,l_j\right)=\left[\begin{array}{cc}i{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\nu }}_{t,\mathrm{1}}|\mathit{\kappa }{|}^{\mathrm{2}}+\left(D_{\mathrm{11}}+C_{\mathrm{11}}\right)/H_{\mathrm{1}}& -D_{\mathrm{11}}/H_{\mathrm{1}}\\ f_{\mathrm{c}}+\left(D_{\mathrm{21}}+C_{\mathrm{21}}\right)/H_{\mathrm{1}}& -D_{\mathrm{21}}/H_{\mathrm{1}}\\ -D_{\mathrm{11}}/H_{\mathrm{2}}& i{\mathit{\sigma }}_{\mathrm{2}}+{\mathit{\nu }}_{t,\mathrm{2}}|\mathit{\kappa }{|}^{\mathrm{2}}+D_{\mathrm{11}}/H_{\mathrm{2}}\\ -D_{\mathrm{21}}/H_{\mathrm{2}}& f_{\mathrm{c}}+D_{\mathrm{21}}/H_{\mathrm{2}}\\ k_i{H}_{\mathrm{1}}\mathrm{\Phi }/{\mathit{\sigma }}_{\mathrm{1}}& k_i{H}_{\mathrm{2}}\mathrm{\Phi }/{\mathit{\sigma }}_{\mathrm{2}}\end{array}\right.\\ & \left.\begin{array}{ccc}-f_{\mathrm{c}}+\left(D_{\mathrm{12}}+C_{\mathrm{12}}\right)/H_{\mathrm{1}}& -D_{\mathrm{12}}/H_{\mathrm{1}}& i{k}_i/\mathit{\rho }\\ i{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\nu }}_{t_{\mathrm{1}}}|\mathit{\kappa }{|}^{\mathrm{2}}+\left(D_{\mathrm{22}}+C_{\mathrm{22}}\right)/H_{\mathrm{1}}& -D_{\mathrm{22}}/H_{\mathrm{1}}& i{l}_j/\mathit{\rho }\\ -D_{\mathrm{12}}/H_{\mathrm{2}}& -f_{\mathrm{c}}+D_{\mathrm{12}}/H_{\mathrm{2}}& i{k}_i/\mathit{\rho }\\ -D_{\mathrm{22}}/H_{\mathrm{2}}& i{\mathit{\sigma }}_{\mathrm{2}}+{\mathit{\nu }}_{t,\mathrm{2}}|\mathit{\kappa }{|}^{\mathrm{2}}+D_{\mathrm{22}}/H_{\mathrm{2}}& i{l}_j/\mathit{\rho }\\ l_j{H}_{\mathrm{1}}\mathrm{\Phi }/{\mathit{\sigma }}_{\mathrm{1}}& l_j{H}_{\mathrm{2}}\mathrm{\Phi }/{\mathit{\sigma }}_{\mathrm{2}}& \mathrm{1}/\mathit{\rho }\end{array}\right]\end{array}\end{array}$$

and

$$\begin{array}{}\text{(C3)}& s\left(k_i,l_j\right)=\left[\begin{array}{c}{\widehat{u}}_{\mathrm{1}}\left(k_i,l_j\right)\\ {\widehat{u}}_{\mathrm{2}}\left(k_i,l_j\right)\\ {\widehat{v}}_{\mathrm{1}}\left(k_i,l_j\right)\\ {\widehat{v}}_{\mathrm{2}}\left(k_i,l_j\right)\\ \widehat{p}\left(k_i,l_j\right)\end{array}\right],b\left(k_i,l_j\right)=\left[\begin{array}{c}{\widehat{f}}_x\left(k_i,l_j\right)\\ {\widehat{f}}_y\left(k_i,l_j\right)\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right].\end{array}$$

In the above equations, ${\widehat{f}}_x$ and ${\widehat{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), $|\mathit{\kappa }{|}^{\mathrm{2}}$ is equal to ($k_i^{\mathrm{2}}+l_j^{\mathrm{2}}$), and σ₁ and σ₂ can be evaluated as (k_iU₁+l_jV₁) and (k_iU₂+l_jV₂) respectively and are referred to as the intrinsic frequencies in the wind farm and upper layers (Smith, 2010).

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

$$\begin{array}{}\text{(C4)}& \begin{array}{rl}& s\left(k_i,l_j\right)=\left[\begin{array}{c}{\widehat{u}}_{\mathrm{1}}^{\mathrm{b}k}\left(k_i,l_j\right)\\ {\widehat{u}}_{\mathrm{2}}^{\mathrm{b}k}\left(k_i,l_j\right)\\ {\widehat{v}}_{\mathrm{1}}^{\mathrm{b}k}\left(k_i,l_j\right)\\ {\widehat{v}}_{\mathrm{2}}^{\mathrm{b}k}\left(k_i,l_j\right)\end{array}\right],\\ & b\left(k_i,l_j\right)=\left[\begin{array}{c}-i{k}_i\widehat{p}\left(k_i,l_j\right)/\mathit{\rho }\\ -i{l}_j\widehat{p}\left(k_i,l_j\right)/\mathit{\rho }\\ -i{k}_i\widehat{p}\left(k_i,l_j\right)/\mathit{\rho }\\ -i{l}_j\widehat{p}\left(k_i,l_j\right)/\mathit{\rho }\end{array}\right].\end{array}\end{array}$$

After solving the two linear system associated with matrices A^3LM and A^3LMR in cascade, ${\widehat{\mathit{u}}}_{\mathrm{1}}^{\mathrm{b}k}$ components are transformed back into physical space, and Eqs. (25)–(27) are used to reconstruct the background velocity field U_b(x). The latter is then interpolated at turbine quadrature points from the 3LM grid by means of bilinear interpolation formulas.