The Jensen wind farm parameterization

Wind farm power production is known to be significantly affected by turbine wakes. When mesoscale numerical models are used to predict power production, the turbine wakes cannot be resolved directly because they are sub-grid features, and therefore their effects need to be parameterized. Here we propose a new wind farm parameterization that is based on the Jensen model, a well-known analytical wake model that predicts the expansion and wind speed of an ideal wake. The Jensen parameterization is implemented and inserted into two commonly used atmospheric numerical models: the Weather Research and Forecasting (WRF) model (herein referred to as just “WRF”) and the Model for Prediction Across Scales (MPAS). In addition, the internal variability in wind speed and direction within a wind farm, the wind direction uncertainty, and the superposition of multiple wakes are taken into account with an innovative approach. The proposed approach and parameterization are tested against observational data at two offshore wind farms: Lillgrund (small in size and tightly spaced) and Anholt (large and widely spaced). Results indicate that power production is predicted more accurately with the Jensen wind farm parameterization than with the Fitch wind farm parameterization, which is the only one available in WRF. Power predictions with the Jensen parameterization are similar in WRF and MPAS. The sensitivity to grid resolution is small, and the bias is generally low and negative. In conclusion, we recommend that the Jensen wind farm parameterization be used in WRF and MPAS, especially for coarse resolution, high turbine density, and wind directions aligned with the turbine columns.

1 Introduction

As modern wind farms increase in quantity and size and as wind turbine rotors expand in diameter, understanding their aerodynamic wakes becomes more critical. When a wind turbine wake hits a downstream turbine, it can cause a significant reduction in its power production; these wake losses negatively impact wind farm power production (Archer et al., 2018). Barthelmie et al. (2009) report that, in large offshore wind farms, wake losses are 10 % to 20 % of the total power output. These power losses are even more significant for wind farms with tightly spaced turbines (Dahlberg, 2009). Therefore, despite the progress made in understanding wakes and wake losses (Fleming et al., 2019; Archer and Vasel-Be-Hagh, 2019; Simley et al., 2020; Stevens and Meneveau, 2020; Johlas et al., 2020; Zong and Porté-Agel, 2020; Nouri et al., 2020; Wu and Archer, 2021), an improved understanding of wakes and a more accurate modeling of their impacts are still of significant practical interest for predicting wind farm power production, developing optimal layouts and control strategies, and quantifying the potential unintended impacts of wind farms on the surrounding environment.

Due to the increase in computational power in recent years, it has been possible to study turbine wakes with numerical simulations. Mainly two numerical approaches are employed: computational fluid dynamics (CFD) and mesoscale modeling. Examples of CFD are Reynolds-averaged Navier–Stokes (RANS) models with various levels of sophistication, from 3D with actuator disks (van der Laan et al., 2015) to parabolic (Iungo et al., 2018), linearized, unsteady, or 2D (see Göçmen et al., 2016, for a review), and large-eddy simulation (LES) with actuator disks or lines. High-resolution LES is the most accurate because it solves the fine-scale details of the wakes around the turbines at a grid resolution of 10 m or less (Lu and Porté-Agel, 2011; Archer et al., 2013; Xie and Archer, 2015; Wu and Porté-Agel, 2015; Xie and Archer, 2017). However, LES is highly computationally expensive and is conducted in an idealized environment by prescribing the inflow characteristics, which prohibit its applications to non-idealized, real-world simulations that span over multiple days or over large wind farms.

Mesoscale modeling, on the other hand, is much less computationally demanding than CFD and is often applied to real-world cases where the two-way interaction between the atmospheric boundary layer (ABL) and the wind farms is taken into account. It therefore has more practical applications, such as estimation of the annual energy production of a wind farm or prediction of temperature changes caused by the wakes. In mesoscale numerical simulations, usually with horizontal resolution of the order of kilometers and vertical resolution of the order of tens of meters within the ABL, a wind turbine is often parameterized as an elevated momentum sink (Volker et al., 2015) or as an elevated momentum sink and a source of turbulence within the vertical levels of the turbine rotor disk (Fitch et al., 2012; Abkar and Porté-Agel, 2015; Pan and Archer, 2018).

One of the most widely used wind farm parameterizations is the Fitch scheme (Fitch et al., 2012) in the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008), hereafter referred to as “Fitch-WRF”, which has been used widely to investigate the wakes of wind farms and their impacts (Fitch et al., 2013a, b; Byrkjedal et al., 2014; Fitch, 2015; Jiménez et al., 2015; Eriksson et al., 2015; Volker et al., 2017; Eriksson et al., 2017; Pryor et al., 2018; Pan et al., 2018; Lundquist et al., 2019). In two recent studies, Pryor et al. (2019) and Shepherd et al. (2020) conducted high-resolution mesoscale simulations to analyze the performance of Fitch-WRF in modeling the downstream wake effects and impact of wind turbine arrays on near-surface conditions. Xia et al. (2019) used Fitch-WRF to understand the underlying mechanisms of wind-farm-induced changes in near-surface temperature over west central Texas. Lee and Lundquist (2017) evaluated the performance of the Fitch-WRF in various atmospheric conditions for a wind farm of 200 wind turbines of 1.5 MW in central Iowa. They reported that meteorological conditions and vertical grid resolution significantly affect the performance of the model. We note that a recent paper by Archer et al. (2020) identified two issues with Fitch-WRF (a code bug in the way the Fitch parameterization is inserted in the WRF model and the excessive value of a coefficient used to calculate the turbine-added turbulence) that likely affected past studies that used WRF v3.6 to v4.2.

Despite its wide adoption and use, Fitch-WRF has been found to generally underestimate wake losses and overpredict the power output of wind farms, especially when the wind is aligned with the turbine columns (Jiménez et al., 2012; Pan and Archer, 2018). As discussed extensively in Pan and Archer (2018), the underlying issue with Fitch-WRF is that, just like most other wind farm parameterizations that have been proposed in the literature, it treats all the turbines that are positioned in the same grid cell in the same way. The individual turbine coordinates are just used to assign each turbine to the center of a grid cell, effectively neglecting the turbine layout within each grid cell. Hence, all turbines in a grid cell are subject to the same inflow wind speed (i.e., spatially averaged and ensemble-averaged velocity in the grid cell), regardless of their position, and their wakes, which are sub-grid features, are neglected. Eriksson et al. (2015) suggested increasing the horizontal grid resolution to 333 m or finer to improve the wake results with Fitch-WRF, which is unfeasible for large-scale or long-term simulations (Rai et al., 2019). Abkar and Porté-Agel (2015) attempted to account for the wind farm layout by introducing a correction parameter ξ. This parameter ξ, however, is sensitive to the wind farm density and wind farm layout (e.g., aligned vs. staggered) and has to be obtained from ad hoc LES results.

To explicitly take into account the layout of the wind farm, here we propose to incorporate an analytical wake loss model into the wind farm parameterization. Analytical wake loss models are simplified representations of the wakes that are based on analytical equations for the wind speed deficit (Archer et al., 2018). The only prior study that incorporated an explicit expression for the wake deficit of each wind turbine was Pan and Archer's 2018 hybrid parameterization, which was based on the geometric model by Ghaisas and Archer (2016). Although the hybrid parameterization performed very well when coupled with the WRF model (herein referred to as just “WRF”), the applicability of the geometric model to any wind farm and any wind turbine is questionable, as it was calibrated based on one specific wind farm (Lillgrund in Sweden) and one wind turbine (Siemens 2.3 MW).

In this study we develop and apply a new wind farm parameterization, based on the Jensen model (Jensen, 1983; Katic et al., 1986), also known as Park model (Peña et al., 2014), to two mesoscale models: WRF and the Model for Prediction Across Scales (MPAS; Skamarock et al., 2012). The Jensen wake model was selected for this parameterization because it is possibly the most widely used analytical wake loss model (Katic et al., 1986; Kirchner-Bossi and Porté-Agel, 2018; Staid et al., 2018; Rivas et al., 2009; Vasel-Be-Hagh and Archer, 2017) and because it performs reasonably well regardless of the wind turbine layout or wind direction (Gaumond et al., 2014; Keane et al., 2016; Tian et al., 2017; Ritter et al., 2017; Archer et al., 2018; Ge et al., 2019). In particular, Archer et al. (2018) evaluated the performance of six popular analytical wake loss models – namely Jensen, Larsen (Larsen, 1988), Frandsen (Frandsen et al., 2006), two Gaussian models (Xie and Archer, 2015; Bastankhah and Porté-Agel, 2014), and the geometric model (Ghaisas and Archer, 2016) – using field data collected at three utility-scale wind farms: Lillgrund in Sweden, which is a mid-sized, closely spaced, offshore wind farm with a regular layout; Nørrekær in Denmark, which is a small, moderately spaced, inland wind farm with a regular layout; and Anholt in Denmark, which is a large, widely spaced, offshore wind farm with an irregular layout. Every analytical wake loss model's performance varied from farm to farm and even wind direction to wind direction. While the Jensen model was not the best model all the time, it stood out for its consistently strong performance and for rarely ranking last for all directions and all farms. The Jensen model assumes a top-hat distribution of the velocity deficit in every turbine wake and then applies superposition methods to account for the interaction among multiple wakes.

Like all other analytical wake loss models, the Jensen model was developed based on one implicit common assumption: that the upstream undisturbed wind speed and direction are the same for all turbines within the wind farm. With an increase in the size of modern wind farms, however, significant variability in the distribution of wind speed and direction within a wind farm is expected, due to surface heterogeneity and mesoscale weather systems (van der Laan et al., 2017; Peña et al., 2018). Neglecting horizontal variability within large wind farms could introduce inaccuracies into the annual energy production and power density predictions. This study addresses this issue and accounts for the variability of wind speed and direction using an innovative approach. Additionally, four wake superposition methods, including one proposed here for the first time, are examined in combination with the new Jensen wind farm parameterization.

The impacts of turbines on the flow are still parameterized as an elevated momentum sink and a turbulence source, like those in Fitch-WRF. This study validates the Jensen parameterization's performance against observational data collected at two offshore commercial wind farms: Lillgrund (small in size and tightly spaced) and Anholt (large and widely spaced).

2 Framework of the wind farm parameterization in WRF and MPAS

2.1 The Fitch wind farm parameterization

In the Fitch parameterization, wind turbines are represented as elevated drag elements that reduce the wind speed and produce turbulent kinetic energy (TKE) at each vertical level k that intersects the rotor. The momentum sink and TKE source terms induced by the turbines in a grid cell are proportional to the fractional rotor area contained in that level (A_k) and to the grid-cell horizontal wind speed at that level (U_k):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u_k}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{N_{\mathrm{t}}}{A_{\mathrm{cell}}}}{\displaystyle \frac{A_k{C}_{\mathrm{T}}{U}_k{u}_k}{\left(z_{k+\mathrm{1}}-z_k\right)}},\end{array}$$

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial v_k}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{N_{\mathrm{t}}}{A_{\mathrm{cell}}}}{\displaystyle \frac{A_k{C}_{\mathrm{T}}{U}_k{v}_k}{\left(z_{k+\mathrm{1}}-z_k\right)}},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{TKE}}_k}{\partial t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{N_{\mathrm{t}}}{A_{\mathrm{cell}}}}{\displaystyle \frac{A_k{C}_{\mathrm{TKE}}{U}_k^{\mathrm{3}}}{\left(z_{k+\mathrm{1}}-z_k\right)}}.\end{array}$$

The power P generated by the turbines in a grid cell is estimated as

$$\begin{array}{}\text{(4)}& P={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{N}_{\mathrm{t}}\mathit{\rho }A{C}_{\mathrm{P}}{U}_{\mathrm{h}}^{\mathrm{3}}.\end{array}$$

In these equations, u_k and v_k are the horizontal wind components, z_k is the height of vertical level k, N_t is the number of turbines in the grid cell, A_cell is the horizontal cross-sectional area of the grid cell (note that A_cell is calculated differently in WRF and MPAS with different grid shapes), ρ is the air density (set to a constant, 1.23 kg m⁻³), A is the turbine rotor area (equal to $\frac{\mathit{\pi }}{\mathrm{4}}{D}^{\mathrm{2}}$, where D is the rotor diameter), U_h is the hub-height wind speed that is interpolated from the horizontal wind speed at the vertical levels surrounding the hub height, C_P and C_T are the power and thrust coefficients (prescribed functions of U_h and dependent on the turbine model), and $C_{\mathrm{TKE}}=C_{\mathrm{T}}-C_{\mathrm{P}}$ is the so-called TKE coefficient. We note that C_TKE as defined in the Fitch parameterization results in an overestimate of TKE (Abkar and Porté-Agel, 2015; Pan and Archer, 2018); Archer et al. (2020) proposed that C_TKE should be revised to one-quarter of the original value. This revised C_TKE has been added in WRF v4.2.1 and is used in this study.

The turbine-induced TKE term from Eq. (3) is directly added to the Mellor–Yamada–Nakanishi–Niino level 2.5 (MYNN) planetary boundary layer (PBL) scheme (Nakanishi and Niino, 2009) in WRF, as Fitch-WRF only works in combination with this particular PBL scheme. The Fitch wind farm parameterization has not been implemented in MPAS v7.0 (the latest version). We inserted it in MPAS in this study. It should be noted that the added TKE is improperly advected in WRF with the original Fitch parameterization, due to a code bug in versions v3.6 to v4.2, but it is fixed in v4.2.1 (Archer et al., 2020). We use the latest version of the models, WRF v4.2.1 and MPAS v7.0 (with a modification in the added TKE advection), in this study.

2.2 The Jensen wind farm parameterization

The Jensen wind farm parameterization consists of two steps: the estimation of the wind speed deficit of a single turbine wake (via the Jensen wake model) and the wake superposition method to account for the interaction and overlapping of multiple wakes.

2.2.1 The Jensen wake model

The normalized wind speed deficit δ is defined as

$$\begin{array}{}\text{(5)}& \mathit{\delta }\left(x\right)={\displaystyle \frac{\mathrm{\Delta }U\left(x\right)}{U_{\mathrm{\infty }}}}={\displaystyle \frac{U_{\mathrm{\infty }}-U\left(x\right)}{U_{\mathrm{\infty }}}},\end{array}$$

where ΔU is the wind speed deficit, U_∞ is the undisturbed hub-height wind speed, and U is the (reduced) wind speed in the wake that is a function of the downwind distance x from the turbine. The Jensen wake model assumes a top-hat distribution of the wind speed deficit in the wake, meaning that the (reduced) wind speed in the wake is assumed to be uniform along y and along z within the edges of the wake itself for each downstream distance x (i.e., the wake is three-dimensional but axisymmetric). The hub-height wind speed at turbine i caused by the wake of turbine j is expressed as

$$\begin{array}{}\text{(6)}& U_{ij}=U_{\mathrm{\infty }}\left(\mathrm{1}-{\mathit{\delta }}_{ij}\right)=U_{\mathrm{\infty }}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{2}a}{{\left(\mathrm{1}+\mathrm{2}{k}_{\mathrm{w}}\frac{x_{ij}}{D}\right)}^{\mathrm{2}}}}\right],\end{array}$$

where x_ij is the along-wind distance between turbines i and j; k_w is the rate of wake expansion (equal to 0.075 and 0.04 for onshore and offshore wind farms, respectively, Archer et al., 2018); and, after applying momentum theory, the induction factor a can be related to the thrust coefficient C_T by

$$\begin{array}{}\text{(7)}& a={\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-C_t}}{\mathrm{2}}}.\end{array}$$

If turbine i with rotor area A_i is not perfectly aligned with the upstream turbine j along the wind direction, the Jensen model needs a modification to account for the fact that only a portion of its rotor (A_0,ij) is affected by the wake of turbine j, while the rest of the rotor ($A_i-A_{\mathrm{0},ij}$) experiences the undisturbed wind speed U_∞ (Archer et al., 2018):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{U}_{ij}& =U_{\mathrm{\infty }}{\displaystyle \frac{A_i-A_{\mathrm{0},ij}}{A_i}}+U_{\mathrm{\infty }}\left(\mathrm{1}-{\mathit{\delta }}_{ij}\right){\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\\ & =U_{\mathrm{\infty }}-{\mathit{\delta }}_{ij}{U}_{\mathrm{\infty }}{\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}.\end{array}\end{array}$$

The wind speed provided by Eq. (8) is effectively a rotor-average wind speed (also known as rotor-equivalent or rotor-layer wind speed), which is a better representation of the wind speed experienced by the entire rotor than the wind speed at the exact location of the hub, and ultimately provides a better estimate of the power production of the turbine (Choukulkar et al., 2016; St. Pé et al., 2018).

Lastly, the wake predicted by the Jensen model is conical and three-dimensional, but, in its original formulation, the only directions that matter are x, along which the wind speed decreases linearly in the wake, and y, along which the wind speed deficit is constant within the lateral wake edges and zero outside of them. This two-dimensional horizontal slice of the wake is then repeated at all vertical levels within the wake cone. In addition to being axisymmetric, the Jensen model was originally formulated for a shearless wind flow, in which the undisturbed wind speed did not vary with height but was equal to U_∞ at all levels. However, a shearless flow hardly ever occurs in reality. To properly account for vertical wind shear, the wake can no longer be treated as axisymmetric. Here the wind speed deficit as a function of z is obtained as follows. The undisturbed hub-height wind speed U_∞ is used to calculate the hub-height (reduced) wind speed in the wake U_i, and then the ratio $U_i/U_{\mathrm{\infty }}$ is used to multiply the wind speed at each vertical level k within the cone cross section. By doing so, the proposed Jensen model is effectively three-dimensional and accounts for vertical wind shear; thus, it can be perfectly integrated in WRF.

2.2.2 The wake superposition methods

When multiple wakes from multiple turbines j (j=1 … N) overlap at turbine i, the incoming flow speed for turbine i is calculated by a wake superposition method. A review of different methods is given in Porté-Agel et al. (2020). The two most common superposition methods of wind speed deficits are a linear superposition (Lissaman, 1979), hereafter referred to as “M1”,

$$\begin{array}{}\text{(9)}& \mathrm{M}\mathrm{1}:U_i=U_{\mathrm{\infty }}-\sum _{j=\mathrm{1}}^N\left[{\mathit{\delta }}_{ij}{U}_{\mathrm{\infty }}{\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\right],\end{array}$$

and a squared superposition (Katic et al., 1986), hereafter referred to as “M2”,

$$\begin{array}{}\text{(10)}& \mathrm{M}\mathrm{2}:U_i=U_{\mathrm{\infty }}-\sqrt{\sum _{j=\mathrm{1}}^N\left[{\mathit{\delta }}_{ij}^{\mathrm{2}}{U}_{\mathrm{\infty }}^{\mathrm{2}}{\left({\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\right)}^{\mathrm{2}}\right]}.\end{array}$$

At the upwind turbine j, the incoming flow speed U_j is not necessarily equal to the undisturbed wind speed U_∞, since turbine j itself may also be affected by wakes. As such, an alternative method (Voutsinas et al., 1990), hereafter referred to as “M3”, can be used to estimate U_i based on a slight modification of Eq. (10):

$$\begin{array}{}\text{(11)}& \mathrm{M}\mathrm{3}:U_i=U_{\mathrm{\infty }}-\sqrt{\sum _{j=\mathrm{1}}^N\left[{\mathit{\delta }}_{ij}^{\mathrm{2}}{U}_j^{\mathrm{2}}{\left({\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\right)}^{\mathrm{2}}\right]}.\end{array}$$

We propose a new wake overlap method, hereafter referred to as “M4”, based not on the superposition of the wind speed deficits but the resulting wind speeds as follows:

$$\begin{array}{}\text{(12)}& \mathrm{M}\mathrm{4}:U_i=\sqrt{{\displaystyle \frac{\sum _{j=\mathrm{1}}^N{U}_{ij}^{\mathrm{2}}}{N}}},\end{array}$$

where U_ij is given by Eq. (8).

Figure 1 Example of how wake overlapping is accounted for with method M4. The four turbines T₁–T₄ are in different grid cells; T₄ is affected by three wakes, each affecting a different portion of its rotor (the thick black line is the affected portion of the rotor, and the thick green line is the undisturbed portion) and causing a different wind speed; the resulting wind speed is proportional to the summation of the wind speeds caused by each upstream turbine i, i.e., V_4,i, using Eq. (12).

Note that these four methods are not derived from fundamental principles, and therefore none of them conserves mass, energy, or momentum. Methods that are based on the linear or squared superposition of the wind speed deficits, like M1–M3, may, at least in principle, cause negative wind speeds because each additional wake further reduces the resulting wind speed. The advantages of method M4 are that it will never cause negative wind speeds and that it is well suited for overlapping wakes that come from different grid cells with different upstream wind speeds and wind directions. Looking at Fig. 1 as an example, turbine 4 is affected by three wakes, each of which causes its own wind speed (e.g., V_4,1 is caused by turbine 1 and is a partial wake case, with some of the rotor affected by the wake and the rest unaffected, thus experiencing V_4,0). The resulting wind speed at turbine 4 is the result of the mixing of the three individual wind speeds. With M4 (Eq. 12), the resulting wind speed is basically the average of the (partially recovered and therefore relatively high) wind speed from the farthest turbine and that (fresh and relatively low) of the nearest turbine, plus all the others in between. In a sense, M4 can be thought of as a way to indirectly account for partial wake recovery due to the added turbulence from multiple wakes and entrainment.

On the other hand, method M4 may underestimate the deficit in the case of perfectly aligned turbines because it tends to dilute the wakes of the nearest turbines with the partially recovered wakes of those further upstream. Even in the aligned case, however, M4 is effective because it avoids the unrealistic continuous drop in wind speed as more and more turbines are aligned, which other wake superposition methods suffer from. While every wake superposition strategy introduces some level of underestimation or overestimation, comparison with observational data, discussed in the next sections, indicates that M4 is generally the most accurate.

In the Fitch parameterization, there is no wake overlapping because each turbine in a grid cell is treated the same, i.e., as a front-row turbine with the same incoming flow speed U_∞, and because there are no wakes at all inside the grid cell. Thus, in Eqs. (1) through (4), the tendency and power at a wind turbine are simply multiplied by the number of turbines in the grid cell N_t to give the overall contributions of the turbines. But in the Jensen parameterization, a slight modification is needed because each wind turbine in a grid cell is affected by the upstream wakes differently. Specifically, the wind speed at each vertical level that intersects the rotor (U_k) is not the same at all turbines; thus, U_k in Eqs. (1) to (3) is replaced by $U_i/U_{\mathrm{h}}\times U_k$,

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle \frac{\partial u_k}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{t}}}{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{cell}}}}{\displaystyle \frac{A_k{C}_{\mathrm{T}}\left(\frac{U_i}{U_{\mathrm{h}}}{U}_k\right)u_k}{\left(z_{k+\mathrm{1}}-z_k\right)}},\text{(14)}& {\displaystyle }& {\displaystyle \frac{\partial v_k}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{t}}}{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{cell}}}}{\displaystyle \frac{A_k{C}_{\mathrm{T}}\left(\frac{U_i}{U_{\mathrm{h}}}{U}_k\right)v_k}{\left(z_{k+\mathrm{1}}-z_k\right)}},\text{(15)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{TKE}}_k}{\partial t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{t}}}{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{cell}}}}{\displaystyle \frac{A_k{C}_{\mathrm{TKE}}{\left(\frac{U_i}{U_{\mathrm{h}}}{U}_k\right)}^{\mathrm{3}}}{\left(z_{k+\mathrm{1}}-z_k\right)}},\end{array}$$

and U_h in Eq. (4) is replaced by the rotor-equivalent wind speed U_i,

$$\begin{array}{}\text{(16)}& P={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{t}}}\mathit{\rho }A{C}_{\mathrm{P}}{U}_i^{\mathrm{3}}.\end{array}$$

Note that the thrust and power coefficients are evaluated using the local U_i and that the grid-cell hub-height wind speed U_h coincides with U_∞,i. Alternative choices could be made for U_∞,i, such as the wind speed at a grid cell located at some distance upstream of the wind farm along one wind direction (among the varying wind directions simulated inside the farm) or the wind speed from offline simulations without the wind farm. However, any alternative would be even more arbitrary and would potentially introduce more errors than using simply the local grid cell U_h. Another advantage is that U_h is straightforward to implement in the codes and fast to execute.

Note that, if the wind turbines in a grid cell are not affected by any upstream turbine wakes, the Jensen wind farm parameterization predicts exactly the same result as the Fitch parameterization.

2.2.3 Treatment of wind speed variability

In their original formulations, methods M1 to M4 are all based on the assumption that the undisturbed wind speed is uniform for all turbines in the entire wind farm, and thus the equations for M1–M4 (Eqs. 9–12) contain only one value of U_∞. When they are implemented in a mesoscale numerical model, such as WRF and MPAS, the grid-cell horizontal wind speed is used as the undisturbed speed U_∞. If the wind farm of interest is so large that it cannot be entirely contained in a single grid cell of the numerical domain or if the grid resolution is fine, the undisturbed speed U_∞ is replaced with U_∞,i in the first term and with U_∞,j in the second term of Eqs. (9)–(12), which are the grid-cell wind speed at the grid cell of turbine i and j at the beginning of the time step. Thus, the multi-cell versions of M1 to M4 are expressed as

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }\mathrm{M}\mathrm{1}:U_i=U_{\mathrm{\infty },i}-\sum _{j=\mathrm{1}}^N\left[{\mathit{\delta }}_{ij}{U}_{\mathrm{\infty },j}{\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\right],\text{(18)}& {\displaystyle }& {\displaystyle }\mathrm{M}\mathrm{2}:U_i=U_{\mathrm{\infty },i}-\sqrt{\sum _{j=\mathrm{1}}^N\left[{\mathit{\delta }}_{ij}^{\mathrm{2}}{U}_{\mathrm{\infty },j}^{\mathrm{2}}{\left({\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\right)}^{\mathrm{2}}\right]},\text{(19)}& {\displaystyle }& {\displaystyle }\mathrm{M}\mathrm{3}:U_i=U_{\mathrm{\infty },i}-\sqrt{\sum _{j=\mathrm{1}}^N\left[{\mathit{\delta }}_{ij}^{\mathrm{2}}{U}_j^{\mathrm{2}}{\left({\displaystyle \frac{A_{\mathrm{0},ij}}{A_i}}\right)}^{\mathrm{2}}\right]},\text{(20)}& {\displaystyle }& {\displaystyle }\mathrm{M}\mathrm{4}:U_i=\sqrt{{\displaystyle \frac{\sum _{j=\mathrm{1}}^N{\left[U_{\mathrm{\infty },i}\frac{A_i-A_{\mathrm{0},ij}}{A_i}+U_{\mathrm{\infty },j}\left(\mathrm{1}-{\mathit{\delta }}_{ij}\right)\frac{A_{\mathrm{0},ij}}{A_i}\right]}^{\mathrm{2}}}{N}}}.\end{array}$$

As discussed later, all wake superposition methods, M1–M4, are inserted in WRF and MPAS, and their performance is compared. Note that the multi-cell wake superposition methods also address the issue of wind speed variability within the wind farm because they directly use the wind speed at each grid cell. Also, Eqs. (13)–(16) do not need any modifications in single cell cases.

2.2.4 Treatment of wind direction variability and uncertainty

When turbines in a wind farm are located in multiple grid cells, the variability in wind direction within the wind farm needs to be taken into account (note that the variability in wind speed was addressed in the previous section with the multi-cell wake superposition methods). Here we propose an innovative approach that considers each individual turbine and all the relevant upstream turbines that might affect it. As displayed in Fig. 2, each grid cell within the wind farm contains one or more wind turbines and is affected by an upstream wind that may come from a different direction than that in the grid cell itself. To enhance computational speed, only the subset of turbines contained within an angle of ±30^∘ around the wind direction at the turbine of interest (red circle) and within a distance < 20 D are considered (the turbines within the dashed red sector). The wind speed deficit caused by these upstream turbines is calculated along the direction that is the average (green arrow) between the wind vector at the upstream cell (black arrow) and that at the cell of the turbine of interest (red arrow). The calculations are conducted starting from the most upstream turbine and going downstream using the grid-cell hub-height wind components.

Figure 2 Example of how the wind direction variability within the farm is treated in the Jensen wind farm parameterization. The wind turbines are black circles, and the turbine of interest is red. The black arrows are the wind vectors at the grid cells; the red arrow is the wind vector at the grid cell of interest, replicated at the relevant upstream grid cells; and the green arrows are the average of the wind direction at the grid cells and that at the grid cell of interest.

The wind direction defines the path of a turbine wake, which determines the wake conditions (e.g., full vs. partial wake) at the downwind turbines. However, wind direction uncertainty exists in both measurements and numerical models (Gaumond et al., 2014). Wind vanes and sonic anemometers have a wind direction accuracy that is typically ±3 and ±2^∘, respectively (Archer et al., 2016). This means that an observation of wind direction α is as likely to be from any direction between α−3^∘ and α+3^∘ if a vane was used. For a mesoscale simulation, the time step is roughly 10 s, and the wind speed and direction are fixed during this time step. The wind at a time step is essentially regarded as spatially averaged and ensemble averaged within a grid cell. In reality, many factors affect the flow, causing uncertainty in the wind direction in this time step (e.g., turbulent fluctuations, turbine wake meandering, sub-grid phenomena, and uncertainty in the numerical model).

To address the wind direction uncertainty, following Gaumond et al. (2014) and Göçmen et al. (2016), a Gaussian wind direction averaging method is used with a standard deviation of 2^∘. We choose seven angles around the flow direction, ±2.5, ±1.5, ±0.5, and 0^∘, to apply the Gaussian weighted averaging, implying that the Jensen wind farm parameterization is called seven times in a time step – once for each of these seven wind directions. As such, the simulation results are effectively averaged over a window of wind directions that is ±2.5^∘ around the intended wind direction. The output power, momentum sink, and added TKE are all averaged using Gaussian weighted averages. This Gaussian wind direction averaging is not added in the Fitch parameterization because it is insensitive to wind direction. We note that Volker et al. (2015) used nine wind directions in a much wider range (±11.25^∘) and ran nine cases with these wind directions separately. The observations in this study are averaged over a ±2.5^∘ window too but with no Gaussian averaging.

3 Methods and data

3.1 Wind power observations at two operational wind farms

The Jensen and Fitch parameterizations are evaluated against data collected at two offshore wind farms: Lillgrund, small in size and tightly spaced with a regular layout (Fig. 3a), and Anholt, large and widely spaced with an irregular layout (Fig. 3b). The Lillgrund wind farm is located in a narrow straight between Denmark and Sweden and consists of 48 Siemens 2.3 MW wind turbines with rotor diameter D=93 m and hub height H=65 m, for a total installed capacity of approximately 110 MW. The Anholt wind farm is located between Djursland and the island of Anholt in Denmark and consists of 111 Siemens 3.6 MW wind turbines with D=120 m and H=82 m, for a total installed capacity of approximately 400 MW. The measurements at Lillgrund were from approximately 16 months at a frequency of 1 min or less, and those at Anholt were from January 2013 to June 2015 at a frequency of 10 min. No atmospheric stability information was available from the measurements at either wind farm. The observational datasets used here are the same as those in Archer et al. (2018); thus, the readers are referred to that article for a detailed description of the data-cleaning procedure, including the yaw bias corrections. As a result of the data-cleaning process, observations are only available for selected columns along and around several directions of alignment at both farms, displayed in Fig. 3. We also analyzed data for a few directions of non-alignment for each farm, shown in red in Fig. 3, selected because they were relatively close to directions of alignment (i.e., within ±20^∘).

Figure 3 Wind turbine layout at Lillgrund (a) and Anholt (b). The Lillgrund wind farm is also plotted in panel (b) with red dots using the same scale as that for the Anholt wind farm. The dashed black lines indicate the wind directions (black arrows) aligned with wind turbine columns. The red arrows indicate the directions of non-alignment. The dashed red lines are the columns analyzed in this study for the non-alignment directions.

3.2 Simulation setup in WRF and MPAS

In the simulations, we use the latest version of the models: WRF v4.2.1 and MPAS v7.0. To focus on the wind farm parameterization, we perform a series of idealized simulations under neutral stability conditions. Although we recognize that atmospheric stability may impact the evolution, shape, length, and duration of the wakes (Ghaisas et al., 2017; Xie and Archer, 2017), neither the Lillgrund nor the Anholt dataset included any measurements of atmospheric stability, and therefore in this study we performed all the simulations under neutral conditions. For the physical options (note that MPAS uses the same physics packages from WRF), only the surface layer scheme and the PBL scheme are turned on. The surface layer scheme sf_sfclay_physics is set to 5, which provides the momentum drag over a sea surface, with the roughness length given by the Charnock relation (Charnock, 1955). The PBL scheme is the MYNN level 2.5 model (bl_pbl_physics = 5; Nakanishi and Niino, 2009), which is the only available option with the Fitch wind farm parameterization.

The simulations are carried out in a rectangular domain of 144 km × 144 km × 10 km in the x, y, and z directions, respectively. The vertical resolution is 4.4 m near the surface and stretched to a 211 m grid spacing at the domain top, with a total of 62 vertical levels (of which four levels are below the rotor and nine levels intersect the rotor). To investigate the sensitivity to the horizontal resolution, we set up three horizontal resolutions of Δx=4 km (for Lillgrund; 24 km for Anholt, as explained later), Δx=2 km, and Δx=1 km but keep the domain size (i.e., the length and width of the domain) the same. The bottom surface is set as water, with surface roughness calculated in the surface layer scheme. A rigid boundary is used at the top of the domain by setting the vertical velocity to zero. A Rayleigh damping layer is applied on the vertical velocity within the top 2000 m of the domain to prevent reflection of gravity waves. Open radiative lateral boundary conditions are used in WRF, and periodic lateral boundary conditions are used in MPAS since the open boundary conditions are not available. Based on our tests, as the wind farm wake extends less than 50 km downwind, the horizontal domain of 144 km × 144 km in MPAS is large enough to minimize the effect of the reentering of the wind farm wake.

The flow in the domain is driven by a pressure gradient that is balanced by a constant and uniform geostrophic wind, with the Coriolis parameter set to $\mathrm{1.11}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ at a latitude of 50^∘ N. With the imposed geostrophic wind (i.e., u_g=7.6 m s⁻¹, v_g=5.88 m s⁻¹), the steady-state flow has a wind speed of about 8.5 m s⁻¹ at 70 m (i.e., about the hub height) from the wind direction of 225^∘. The flow in the whole domain is initialized with a constant potential temperature profile of 300 K from the surface to a height of 900 m, and then it is linearly increased with height at a constant rate of 3 K km⁻¹. The horizontal wind components are set equal to the geostrophic wind at all levels. The humidity is set to zero. All runs are integrated first for 96 h without the turbines to reach a steady state. Then, the simulation is continued for another 4 h with the wind farm placed at the southwest corner of the domain (i.e., wind farm center at one-third of the domain from the south and the west boundaries).

To ensure a proper comparison of the Jensen vs. Fitch parameterizations in WRF and MPAS, the physical schemes and the domain configurations are set equal in the two models to the extent possible. Because MPAS uses hexagonal grids, here we define the horizontal grid resolution as the distance between the center of two adjacent grids. Due to the different grid shapes in WRF and MPAS (rectangle vs. hexagon), turbines may be partitioned into different grid cells even with the same grid resolution. We also note that the turbine density ($N_{\mathrm{t}}/A_{\mathrm{cell}}$) in WRF and MPAS is different because of the difference in the grid shape (the area of a hexagonal grid is about 15 % smaller than that of a rectangular grid with the same grid resolution). Thus, even with the same configurations and the same wind farm parameterization, the predicted power generation at the same wind farm is expected to be slightly different between the two models.

When the flow in the domain reaches a steady state after 96 h, the wind farm is activated. Based on the wind farm layouts, we choose six wind directions for Lillgrund (i.e., 180, 222, 255, 270, 300, and 315^∘) and eight wind directions for Anholt (i.e., 168, 179, 183, 210, 228, 240, 339, and 341^∘) to run the simulations. Here we rotate the wind farm layout and use the same steady flow with a direction of 225^∘ for all runs rather than performing the simulations initialized with different geostrophic winds to obtain the steady flow with the desired wind directions. The simulations are categorized into two groups based on the grid resolution: single-cell and multi-cell cases. For the single-cell cases, the entire wind farm is contained in a single grid cell, while for the multi-cell cases, the turbines of the wind farm are located in multiple grid cells. For the single-cell cases, we use different grid resolutions at Lillgrund and Anholt, 4 and 24 km, respectively, because of the different wind farm sizes. For the multi-cell cases, we use two grid resolutions, 2 and 1 km, at both farms. As mentioned in Sect. 2.2.4, to account for the wind direction uncertainty, the simulations with the Jensen parameterization and the observations are averaged over a wind direction sector of ±2.5, or 5^∘ wide.

4 Results and discussion

4.1 Wake superposition methods

Since one key factor in the Jensen parameterization is the superposition of multiple wakes, we first examine the four wake superposition methods, M1–M4, for the single-cell cases at Lillgrund and Anholt, and then we choose the best two methods to apply in WRF and MPAS for multi-cell cases. For the single-cell cases, only the results from WRF are displayed because the Jensen parameterization in WRF and MPAS produces identical results due to the same undisturbed wind speed and direction used for all the turbines. For both the single-cell and multi-cell cases, power output from the Fitch parameterization is also shown for comparison.

Relative power is used for the comparative analysis, which is defined as the ratio of the power of each turbine in a column over that of the front-row turbine in the same column. The performance statistics used are bias error and root-mean-square error (RMSE), defined as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(X_i-X_{i,\mathrm{obs}}\right)}{N}},\text{(22)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left(X_i-X_{i,\mathrm{obs}}\right)}^{\mathrm{2}}}{N}}},\end{array}$$

where “obs” represents an observation, X_i is the variable (i.e., relative power) of interest for turbine i, and N is the sample size. These statistics are expressed as a percent in this study.

At Lillgrund, we simulate the power output along eight columns for five directions of alignment (i.e., wind directions under full-wake conditions), 180, 222, 225, 270, and 300^∘, and along one column for one direction of non-alignment (partial wake condition), 315^∘ (Fig. 3a). Since Lillgrund is tightly spaced and regular, there is no direction that is truly of non-alignment; even 315^∘ for the column led by turbine T46 is partially affected by wakes. For the wind directions of 180, 222, and 225^∘, two columns are selected because they represent two distinct situations in which one column includes a gap and the other does not.

At Anholt, we simulate the power output for six directions of alignment, 168, 179, 183, 228, 339, and 341^∘, and two directions of non-alignment (no-wake condition), 210 and 240^∘ (Fig. 3b). Because of the irregular wind farm layout, we choose the wind directions of alignment with as many turbines as possible aligned along them to run the simulations, following Archer et al. (2018).

We note that the power measurements at any wind farm are the result of local wind speeds at the turbines, often in strongly non-homogeneous conditions. By contrast, the model results by design are very diffusive in the horizontal direction, with a true resolution that is coarser than the grid size. As such, a direct comparison of the WRF and MPAS results against the observations reported in the next section should be interpreted with this limitation in mind.

4.1.1 Single-cell cases

Figure 4 shows the relative power from superposition methods M1–M4 in the Jensen parameterization at Lillgrund. The relative power from the Fitch parameterization is also shown for comparison, although it is always equal to 1 since the sub-grid turbine wakes are not considered.

Figure 4 Lillgrund (single cell): relative power from observations and from the Jensen wind farm parameterization with the four wake superposition methods, M1–M4, along the directions of alignment (a–h) and the direction of non-alignment (i). The black circles represent the mean observed relative power, and the bars are 1 standard deviation. The simulation results from the Jensen parameterization, and the observations are averaged over ±2.5^∘; the bin size for the observations is 0.5^∘. The wind direction averaging was not applied to the results from the Fitch parameterization because it is insensitive to wind direction.

For alignment and non-alignment directions, methods M1–M4 predict an identical and sudden drop in relative power at the second turbine. The identical relative power can be derived from Eqs. (9) to (12), which give the same expression when N is equal to 1 for the single-cell version of the equations. However, this power drop at the second turbine is overestimated in all the columns compared with the observations, regardless of the distance between the first and the second turbine. A similar trend was also found in previous studies using the Jensen analytical wake loss model with a narrow wind averaging sector of ±2.5^∘ (Simisiroglou et al., 2019). We note that, for the non-alignment wind direction of 315^∘ (Fig. 4i), methods M1–M4 predict an identical relative power at all turbines (except for the last one), because all but the last one experience only one upstream partial wake.

After the second turbine, the four methods differ significantly. For most columns in Fig. 4, an increase in relative power at the third turbine, a phenomenon caused by the enhanced turbulence mixing after the first turbine that makes the second turbine wake recover faster, is reproduced well by M3 and M4, while M1 and M2 do not capture this feature. Not only at the third turbine but generally, M1 and M2 (M1 in particular) substantially underpredict the power output at all the downwind turbines compared with the observations.

As shown in Table 1, M1 produces the largest overall bias (−18.5 %) and the largest overall RMSE (22.5 %), followed by M2 with a bias of −13.2 % and RMSE of 16.1 %. The reason why M1 and M2 perform worse is that they use the undisturbed wind speed in calculating the velocity deficit (the second term in Eqs. 9 and 10).

Table 1 Bias and root-mean-square error (RMSE) in power prediction by the Fitch parameterization and by the Jensen parameterization with the wake overlap methods, M1, M2, M3, and M4, at Lillgrund for single-cell cases (in percent). The wind turbine columns are shown in Fig. 4a–h, and the leading turbine of each column is indicated in parentheses where there is ambiguity. “N.A.” in parentheses represents the non-alignment wind direction. The results are from WRF.

For the two better-performing methods, M3 and M4, the performance of M3 is consistent at all the downwind turbines at Lillgrund, with most of its predictions lower than the observations, in agreement with Archer et al. (2018), but still within the error bars, whereas the performance of M4 is less consistent. As shown in Fig. 4, M4 tends to underpredict the power output at the near-front turbines (i.e., the second and third turbines) and overpredict the power output at the inner turbines (i.e., fourth, fifth, and following turbines). Additionally, M4 predicts the power increase caused by the gap well, while M3 still underpredicts the power there (Fig. 4b and d). The resulting overall bias for M4 is 2.5 %, while that for M3 is −5.7 % (Table 1). In terms of RMSE, the two superposition methods have close values, with a slightly larger overall RMSE for M4 (9.0 % vs. 10.1 % for M3 vs. M4). In general, M3 and M4 perform similarly at Lillgrund for the single-cell cases, with M4 giving a slightly higher power prediction (slightly positive bias).

Figure 5 presents comparisons of the relative power output from M1–M4 at Anholt. Although the turbine spacing and the turbine properties are different between the two wind farms, M1–M4 perform similarly as they did at Lillgrund. Again, M1 and M2 substantially underpredict the power at all the downwind turbines and perform worse than M3 and M4, with large overall RMSEs (23.0 % and 15.5 % for M1 and M2, respectively, Table 2). M3 tends to underpredict the relative power, while M4 slightly overpredicts it. Both M3 and M4 reproduce well the feature of power output remaining unchanged after the fourth turbine in an alignment column (Fig. 5), due to the balance between the momentum extracted by the turbines and that replenished from the boundary layer. The performance of M4, in particular the overprediction at inner turbines, seems to be improved with larger turbine spacing at Anholt.

Figure 5 Same as Fig. 4 but for Anholt (single cell) with the directions of alignment (a–f) and the directions of non-alignment (g, h).

Table 2 As in Table 1 but at Anholt for single-cell cases (in percent). The results are from WRF.

For the non-alignment directions 210 and 240^∘, the Jensen and the Fitch parameterizations predict the same relative power (all equal to 1 except for the third turbine in the column for 240^∘). It is expected since the two parameterizations are essentially the same for the turbines that are not blocked by any upstream turbine wakes, as we pointed out in Sect. 2.2.

In general, M4 slightly outperforms M3 at Anholt, with a smaller overall RMSE and bias. As at Lillgrund, M4 tends to predict higher relative power than M3 also at Anholt but still lower than the observed.

4.1.2 Multi-cell cases

The two better-performing methods, M3 and M4 (multi-cell versions, Eqs. 17 to 20), in the Jensen parameterization are further examined for multi-cell cases at Lillgrund and Anholt. For the multi-cell cases, the turbines in a wind farm are partitioned into different grid cells, implying that the undisturbed wind speed might be different for each turbine (remember that we use the grid-cell hub-height wind speed as undisturbed wind speed for any turbine). Thus, differences in power output are expected due to the difference in grid resolutions and grid shapes (i.e., rectangle in WRF vs. hexagon in MPAS).

Figures 6 and 7 show relative power from M3 and M4 at Lillgrund with grid resolutions of 2 and 1 km, respectively. The results from the Fitch parameterization are also shown to provide a comparison. Comparing the power output by the Jensen parameterization in WRF and MPAS, it appears that the predicted power output from the two numerical models is fairly close, showing similar overall bias and RMSE at Lillgrund (Table 3). No general trend of underprediction or overprediction by one model over the other was found. Small differences can be found between the two models because the wind turbines may be partitioned into different grid cells, depending on the grid resolution, grid shape, and wind direction (remember that we rotate the wind farm layout to represent different wind directions).

Figure 6 Same as Fig. 4 but for Lillgrund multi-cell cases with Δx=2 km.

Figure 7 Same as Fig. 4 but for Lillgrund multi-cell cases with Δx=1 km.

Table 3 As in Table 1 but at Lillgrund for multi-cell cases (in percent). The results are from WRF and MPAS.

The turbine partitioning affects the undisturbed wind speed in the grid cells, regardless of which wind farm parameterization is used, in two ways. First, the more the turbines partitioned into a specific grid cell, the larger the drag force exerted on the flow in that grid cell, resulting in a lower grid-cell wind speed at the end of each time step (which is the undisturbed wind speed for the turbines in that grid cell at the beginning of the next time step), and vice versa. Second, turbine partitioning affects the undisturbed wind speed via the numerical advection of lower wind speeds into downstream grid cells. As discussed later in Sect. 4.2, the combination of the decrease in wind speed due to the turbine drag and the numerical advection of this reduced wind speed is also called the “resolved wake” effect.

With respect to the performance of M3 and M4 in the Jensen parameterization at Lillgrund, we observe that M4 performs more consistently than M3 in the two numerical models (Figs. 6 and 7), suggesting that M4 is less sensitive to the turbine partitioning compared to M3. For example, at the second turbine in each column, the relative power predicted by Jensen with M4 from WRF and from MPAS is nearly the same (although not identical), but the predictions by M3 differ markedly. When compared to the observations, in line with the single-cell cases, M4 always outperforms M3, particularly at finer grid resolutions. For the directions of non-alignment (Figs. 6i and 7i), more differences between the two numerical models and between the simulated and the observed values are found, especially at the last few turbines, due to the partitioning of the turbines into different grid cells. The reduction in relative power with grid resolution is attributable to the resolved wake effect, which will be discussed in the next section.

Figures 8 and 9 show relative power at Anholt with grid resolutions of 2 and 1 km, respectively. Besides the general trend discussed above (i.e., WRF and MPAS give very similar predictions, and M4 predicts slightly higher relative power, which is closer to the observations, than M3), the figures indicate little sensitivity of the relative power to grid resolution along directions of alignment at Anholt.

Figure 8 Same as Fig. 4 but for Anholt multi-cell cases with Δx=2 km.

Figure 9 Same as Fig. 4 but for Anholt multi-cell cases with Δx=1 km.

In summary, the Jensen wind farm parameterization with M4 still outperforms the others at Anholt, as indicated by the smallest overall bias and RMSE (Table 4). Also, increasing the resolution from 2 to 1 km does not improve the power predictions at Anholt, as the bias and RMSE are about the same or slightly worse.

Table 4 As in Table 1 but at Anholt for multi-cell cases (in percent). The results are from WRF and MPAS.

4.2 Sensitivity to grid resolution

The combination of the wind speed reduction due to the extraction of momentum from the flow caused by each additional turbine in a grid cell (i.e., wind turbine drag) and the numerical advection of this reduced wind speed from upstream into downstream grid cells can be considered as the resolved wind farm wake (Jiménez et al., 2015; Eriksson et al., 2015). When the Fitch parameterization is used, this is the only wake effect that is accounted for; when the Jensen parameterization is used, this effect is still present, but it is added on to that of the sub-grid wakes. This causes a potential double-counting issue of the wake effects, which will be explored next by comparing results at various grid resolutions.

To examine the sensitivity of the resolved wakes to grid resolution, we first analyze the results from the Fitch parameterization, where the situation is less complicated without considering sub-grid wakes. In general, the power predictions with the Fitch parameterization decrease as the grid resolution increases. This sensitivity can be observed in both relative power and total power, each of which will be discussed next.

The progressive decrease in relative power from the Fitch parameterization along long columns of turbines can be appreciated, for example, for the 222^∘ case at Lillgrund (Fig. 7c) or the 228^∘ at Anholt (Fig. 9d), none of which is manifest in the observations. Long columns of turbines are more likely to be partitioned over multiple grid cells; thus, the resolved wake effect is further amplified. Comparing relative power from the Fitch parameterization at 2 km vs. 1 km, either at Lillgrund (Figs. 6 and 7) or Anholt (Figs. 8 and 9), the values are generally lower at the finer resolution, for both alignment and non-alignment directions. For example, at 300^∘ at Lillgrund, the relative power at the last turbine is ∼0.6 at 2 km but ∼0.4 at 1 km for WRF and MPAS; at Anholt, at 228^∘ the relative power at the last turbine is >0.7 at 2 km but ∼0.6 at 1 km for WRF and MPAS. Again, this effect is more pronounced at long columns because of the repeated resolved wake effect (i.e., numerical advection and momentum extraction).

The power production of the entire wind farm (Fig. 10) also decreases as the grid resolution gets finer when the Fitch parameterization is used. Note that in Fig. 10 the undisturbed hub-height wind speed is the same for all wind directions and for all methods. Since the results from MPAS are similar to those from WRF, only the WRF results are shown, but the conclusions drawn are still valid for MPAS. The observational data are not displayed in this figure since they are not available for the entire wind farm, as mentioned in Sect. 3.1. However, the power predictions for the selected columns for which data were available, shown in Fig. 11, will also be discussed later.

Figure 10 Total power (MW) at Lillgrund (a) and Anholt (b) from the Jensen wind farm parameterization in WRF with two wake superposition methods, M3 (red) and M4 (blue), and from the Fitch parameterization (grey) from all the turbines in the farm for the selected wind directions. For each method, single-cell cases (SC) and multi-cell cases with grid resolutions of 2 and 1 km are shown. The figures with “N.A.” in the title are the results for the non-alignment wind directions (315^∘ at Lillgrund; 210 and 240^∘ at Anholt).

For all the wind directions at Lillgrund and Anholt, the power output in the single-cell cases is the largest with Fitch. This is consistent with the single-cell results from the previous sections, where the relative power was (incorrectly) always 100 %. When further increasing the grid resolution to 2 or 1 km, notable differences are observed between the single-cell and multi-cell cases with Fitch, with power production being more than 15 % lower at the finest resolution (multi-cell 1 km). This is because the resolved wake effect causes an incremental reduction of wind speed deep within the wind farm as the grid resolution increases. As shown in Fig. 10b, this trend (i.e., power output decreases with increasing grid resolution) is also evident at Anholt, although the magnitude of the decrease going from 2 to 1 km is smaller than at Lillgrund for two reasons. First, the turbine density ($N_{\mathrm{t}}/A_{\mathrm{cell}}$) is much lower at Anholt than at Lillgrund, and, second, the turbine density at Anholt is about the same at 2 and 1 km due to the wide spacing.

Despite the resolved wakes, the Fitch parameterization substantially overpredicts power along directions of alignment, particularly for the coarse-resolution cases (Fig. 11), suggesting that considering only the resolved wakes might not be adequate. Increasing the grid resolution is definitely beneficial with the Fitch parameterization, especially for directions of alignment, because the reduction in wind speed and power caused by the resolved wakes is amplified with finer resolution. By contrast, for directions of non-alignment (e.g., 315^∘ at Lillgrund and 210 and 240^∘ at Anholt in Fig. 11), the power predictions with the Fitch parameterization are slightly underestimated, consistent with the literature (Jiménez et al., 2015; Pan and Archer, 2018). The underestimate for directions of non-alignment occurs because the resolved wakes incorrectly reduce the value of the hub-height wind speed used to calculate power downstream (“incorrectly” because there are no wake losses for non-alignment directions).

Figure 11 Same as Fig. 10 but for selected directions at Lillgrund and Anholt for which observations were available.

The Jensen parameterization includes the sub-grid wake effects, which are less sensitive to grid resolution and tend to decrease the total power output when compared to Fitch's predictions. In particular at Lillgrund, where the sub-grid wake effects are strong, the wind farm power predicted by the Jensen parameterization with M3 and M4 is much smaller than that by Fitch (Fig. 10). Both M3 and M4 predict power outputs that are much less sensitive to grid resolution than Fitch, especially when the turbine density is high, like at Lillgrund. M4 always predicts a higher power output than M3, which is closer to the observations (Fig. 11); this trend is more pronounced at Lillgrund with strong sub-grid wakes.

Note that, for directions of non-alignment, both Fitch and Jensen perform well at both farms (Fig. 11), and the predictions do not improve with finer grid resolution – not even for Fitch.

Unlike in the Fitch parameterization, the resolved wake effect is small in the Jensen parameterization, due to the much larger effect of the sub-grid wakes. Relative power decreases slightly as the grid resolution increases at Lillgrund (e.g., compare Figs. 6 and 7), due to the resolved wake effect. At Anholt, a small decrease in power output can be detected for all directions when the resolution changes from 2 to 1 km (Figs. 10b and 11b). This is the small signature of the resolved wake effect in Jensen's results.

At Lillgrund, because of the regular layout, the trends in total power as a function of grid resolution (Fig. 10a) are very similar to those from the selected columns for which observational data are available (Fig. 11a). At Anholt, however, since the wind farm layout is irregular, only a small fraction of the turbines are truly aligned with any given wind direction, and therefore the displayed wind farm power output in Fig. 10b includes multiple wake conditions (i.e., full, partial, and no wakes; actually, most turbines are under partial or no wake conditions), resulting in a much higher wind farm efficiency. As such, some of the trends of wind farm output differ slightly from those along the columns with data. For example, for the non-alignment direction of 210^∘, the power in Fig. 11b is basically equal at 2 and 1 km and between M3, M4, and Fitch, whereas it exhibits minor differences in Fig. 10b. Despite these minor inconsistencies, the grid sensitivity of the entire wind farm power is generally the same as that at the selected columns with observations at Anholt too.

We conclude that the Jensen parameterization, especially with M4, performs best at both farms. The Jensen parameterization tends to underestimate the power, regardless of the alignment or non-alignment conditions and regardless of the grid resolution. The consistent sign of the bias (negative) in column power output and the minimal sensitivity to grid resolution or wind direction are all desirable properties. By contrast, the Fitch parameterization is sensitive to the grid resolution and to the alignment or non-alignment of the turbine columns with the wind directions. For directions of alignment, the Fitch parameterization largely overestimates the power for single-cell and multi-cell cases at both farms. At Lillgrund, however, at 1 km, the power output predicted by the Fitch parameterization is close to the observed for wind directions of 255 and 270^∘ (those with the largest spacing, 6 and 8.5 D, respectively), but it is still generally overestimated for the other wind directions. For directions of non-alignment, the Fitch parameterization tends to slightly underestimate the power – more so at fine resolution and at Anholt (i.e., widely spaced farm).

With respect to the wake overlapping methods with the Jensen parameterization, we recommend method M4, although M3 also performs satisfactorily for single-cell cases. We note that M4 consistently predicts higher power output than M3.

4.3 Wind speed and TKE distributions

To obtain a better understanding of how the two wind farm parameterizations affect the power production of a wind farm, we compare the simulated wind speed and TKE at Lillgrund and Anholt. Since the two parameterizations perform similarly in WRF and MPAS, we only show the results from WRF here. The wake superposition method in the Jensen parameterization is M4.

Figure 12 shows the vertical profiles of wind speed and TKE from the Fitch and the Jensen parameterizations from a single-cell case at the grid cell containing the wind farm (Lillgrund or Anholt). The wind direction for the selected case at Lillgrund is 222^∘, and that at Anholt is 168^∘. For these wind directions, the wake effects are neither strongest nor weakest, which allows for a representative comparison. Note that the grid resolution for the single-cell case at Lillgrund is 4 km, while that at Anholt is 24 km, resulting in a significant difference in the turbine density between two wind farms (i.e., 3.0 km⁻² vs. 0.2 km⁻²).

Figure 12 Vertical profiles of wind speed (a, c) and TKE (b, d) by the Fitch and Jensen parameterizations from WRF single-cell cases, 222^∘ at Lillgrund (a, b) and 168^∘ at Anholt (c, d).

Comparing the wind speed profiles from the two wind farm parameterizations, we find a much larger wind speed reduction from the Fitch parameterization at Lillgrund, but, surprisingly, the difference between the two parameterizations is almost indistinguishable at Anholt. This different behavior at Lillgrund and Anholt is attributed to the turbine density and the wind farm layout. At Lillgrund, the turbine spacing is small and the wind farm layout is regular, leading to a large grid-cell turbine drag force on the flow. At Anholt, the turbine spacing is large and the wind farm layout is irregular, such that most of the turbines are not (or partially) influenced by sub-grid turbine wakes; thus, the Jensen parameterization behaves similarly to the Fitch parameterization. As a result, when the turbine-induced drag forces are applied over a large grid cell of 24 km × 24 km at Anholt, the difference in the speed reduction between the two parameterizations is almost indistinguishable with respect to the background wind speed.

For the predicted TKE from the two parameterizations, however, differences are observed at both wind farms, with a larger TKE value from the Fitch parameterization. At Anholt, although the wind speeds from the two wind farm parameterizations are nearly the same (Fig. 12c), the TKE profiles are more notably different (Fig. 12d). For example, the turbine-added TKE at hub height at Anholt is about 0.087 m² s⁻² for the Jensen parameterization and 0.098 m² s⁻² for the Fitch parameterization. Compared to the TKE background value of 0.48–0.51 m² s⁻², the value of the turbine-induced TKE source is much larger than the value of the turbine-induced momentum sink, about 20 % vs. 1 %, for both parameterizations.

We note that the TKE profile at Lillgrund is different from that in Fig. 8a of Pan and Archer (2018) with the same grid resolution and hub-height wind speed. Our result from the Fitch parameterization is much lower (e.g., 0.9 m² s⁻² vs. 1.55 m² s⁻² for the maximum TKE from the two), and the maximum TKE occurs near hub height, while their maximum value is located slightly above hub height. These differences are attributable to the fixed bug in the TKE advection in the WRF code (Archer et al., 2020). Along with the bug fix, we also followed Archer et al.'s 2020 suggestion to reduce the TKE coefficient C_TKE to one-quarter of the original value.

At both Lillgrund and Anholt (Figs. 13 and 14), the wind speed at the grid cell of the wind farm is slightly lower and the TKE is larger in the Fitch than in the Jensen parameterization, in line with the vertical profiles from the single-cell case. The differences between the two parameterizations are more pronounced at Lillgrund, where the wake losses are more significant due to the tighter spacing than at Anholt.

Figure 13 Horizontal cross sections of (a, b) hub-height wind speed (m s⁻¹), (c, d) hub-height TKE (m² s⁻²), and (e, f) power production (MW) at Lillgrund from the Jensen (a, c, e) and Fitch (b, d, f) parameterizations. The results are from WRF with a grid resolution of 1 km and wind direction of 222^∘.

Figure 14 Same as Fig. 13 but at Anholt for wind direction of 168^∘. Note that the flow appears to be from the southeast, but in reality it is from 168^∘, and the wind farm is rotated accordingly.

In general, even though a difference in the wind speed between the two wind farm parameterizations is observed, this difference is small (e.g., ∼0.1 m s⁻¹ at hub height for the single-cell case at Lillgrund). Thus, the large difference in the predicted power output between the two parameterizations is mainly determined by the consideration of sub-grid wakes. Our results indicate that the Jensen parameterization, which accounts for the sub-grid wakes, may tend to underpredict the power output in multi-cell cases.

5 Conclusions, limitations, and future work

We describe and implement a new wind farm parameterization in the WRF and MPAS mesoscale models, based on the classical Jensen model, which takes into account the sub-grid wind turbine wakes, the wind speed and direction variability within a wind farm, and the wind direction uncertainty. Four turbine wake superposition methods are examined in the Jensen parameterization, including an innovative method (i.e., M4) that is based on the superposition of the resulting wind speeds in the wake, as opposed to the resulting deficits. The observational data at two offshore wind farms, Lillgrund and Anholt, are used to evaluate the wind farm parameterization. Due to the lack of atmospheric stability information from the measurements, all the simulations were conducted under idealized neutral conditions.

We conclude that the Jensen wind farm parameterization performs well in WRF and MPAS, as the results of the simulations match the observations, particularly with the new proposed wake superposition method, M4. By contrast, the Fitch wind farm parameterization tends to overpredict the power, especially at coarse grid resolution and for directions of alignment. The simulation results at the two wind farms with totally different grid spacing and layout show similar overall bias and RMSE, suggesting the robustness of the Jensen parameterization. Grid sensitivity tests show that, while the Fitch parameterization tends to predict decreasing relative power at downwind turbines and decreasing total power output with increased grid resolution, due to the resolved wakes, this effect is minimal in the Jensen results. The proposed Jensen wind farm parameterization predicts the power generation at each wind turbine, without increasing too much the computational cost (about 2 %), which makes it a valuable tool for practical applications, such as estimation of the annual energy production and optimization of wind farm layout.

There are several limitations to our study. First, we only considered wind speeds that are below the rated wind speed, thus with a high thrust coefficient; when the wind speed is above the rated wind speed, the thrust coefficient decreases. Since the performance of a superposition method depends on the thrust coefficient (Machefaux et al., 2015), our results might be different in such cases. In addition, we did not test the sensitivity of our results to the value of the standard deviation used for Gaussian averaging.

Another limitation is that we focused only on wind power production for validation. Wind farm parameterizations have been extensively used to study in-farm and farm-to-farm wakes. These wake effects on variables like temperature, humidity, turbulence, or heat fluxes are much harder to evaluate in a parameterization because multiple data sources are needed at and downstream of multiple wind farms, including strategically located tall masts and lidar measurements. Improving the power production alone, as we did in this study, does not automatically improve the characterization of the wind farm wakes.

Lastly, when the wind farm is partitioned over multiple grid cells (i.e., in the multi-cell cases) and the Jensen parametrization is used, there is the possibility of both resolved and sub-grid wakes being present simultaneously in the same grid cell, thus potentially double counting some of the wake effects. By contrast, when Fitch-WRF is used, the resolved wake is the only wake effect that is accounted for, but it is generally too weak. When the Jensen parameterization is used, however, the resolved wake is still present, but it is in addition to the sub-grid wakes, which are generally stronger. Overall, we find that accounting for the (strong) sub-grid wakes with the Jensen parameterization, even in the presence of the inevitable (but small) resolved wake, gives more accurate results than relying on just the resolved wake. However, this issue needs to be investigated further, as discussed by Ma et al. (2022).

Future implementations may consider a more sophisticated formulation of the wake expansion coefficient k_w, using for example turbulence intensity and/or atmospheric stability (Stevens et al., 2016); a more advanced wake loss model, such as one based on the Gaussian model (Bastankhah and Porté-Agel, 2014; Xie and Archer, 2015); and a separate model for the TKE in the wake (one that would not rely on advection alone to inject the correct amount of TKE in the wake of the wind farm), to improve the wind farm power and flow prediction in mesoscale numerical simulations.