Analytical solution for the cumulative wake of wind turbines in wind farms

This paper solves an approximate form of conservation of mass and momentum for a turbine in a wind farm array. The solution is a fairly simple explicit relationship that predicts the streamwise velocity distribution within a wind farm with an arbitrary layout. As this model is obtained by solving flow governing equations directly for a turbine that is subject to upwind turbine wakes, no ad hoc superposition technique is needed to predict wind farm flows. A suite of large-eddy simulations (LES) of wind farm arrays is used to examine self-similarity as well as validity of the so-called conservation of momentum deficit for turbine wakes in wind farms. The simulations are performed with and without the presence of some specific turbines in the wind farm. This allows us to systematically study some of the assumptions made to develop the analytical model. A modified version of the conservation of momentum deficit is also proposed to provide slightly better results at short downwind distances, as well as in the far wake of turbines deep inside a wind farm. Model predictions are validated against the LES data for turbines in both full-wake and partial-wake conditions. While our results highlight the limitation in capturing the flow speed-up between adjacent turbine columns, the model is overall able to acceptably predict flow distributions for a moderately sized wind farm. Finally, the paper employs the new model to provide insights on the accuracy of common wake superposition methods.


Introduction
Fast running engineering (i.e., control-oriented) wake models are arguably still the most popular tools in the wind energy industry for the design and active control of wind farms due to their simplicity and low computational cost (Stevens & Meneveau 2017;Porté-Agel et al. 2020). In order to improve wind farm flow modelling without increasing computational costs, we need to develop a better theoretical understanding of wind flow physics in wind farms (Meneveau 2019). Modelling of wind farm flows based on engineering wake models often consists of two steps: (i) modelling wakes of each individual wind turbine separately and then (ii) using wake superposition techniques to take into account cumulative wake effects in wind farms. The most common superposition techniques are linear (Lissaman 1979) and root-sum-square (Katić et al. 1987), denoted by A.I and A.II in table 1, respectively. For each method, the table shows how the total velocity deficit Δ at a given position depends on the velocity deficit Δ caused by the i ℎ wind turbine (WT ), where changes from 1 to , and is the number of wind turbines upstream  of that position in a wind farm. In addition, different methods are used in the literature to compute Δ , as shown in table 1. By definition, the velocity deficit for WT is equal to the difference between the incoming velocity and the wake velocity denoted by . The incoming velocity, however, can be defined either based on (i) the velocity at the inlet of the wind farm denoted by These different superposition methods are often claimed to conserve flow properties. For instance, the linear superposition method (A.I) is perceived to conserve momentum deficit (Lissaman 1979), while the root-sum-square method (A.II) is thought to conserve the energy deficit (Katić et al. 1987). However, as described later in §7, the validity of these relationships has not been proven, and so these superposition methods should be viewed as empirical relationships. A result-driven approach was mostly adopted to develop such methods. For instance, a superposition method consisting of methods A.I and B.I used by Lissaman (1979) is known to result in negative velocity values in large wind farms. Katić et al. (1987) overcame this issue by proposing a new method including A.II in conjunction with B.I. More recently, Zong & Porté-Agel (2020) have developed a more physics-based superposition method based on the mean convection velocity for each turbine wake. The convection velocity for the combined wake is obtained through an iterative approach. Other recent studies have examined the accuracy of existing superposition models (e.g., Vogel & Willden 2020;Gunn et al. 2016, among others). Overall, there is not unanimous agreement in the literature on which method is most accurate over the wide variety of operating conditions of wind farms.
In this study, an approach that does not rely on the superposition of single turbine wake models is adapted based on a holistic view of turbines in a wind farm. Unlike prior studies that assume each turbine can be treated as a single turbine and separated from the rest of the wind farm, we directly solve governing equations for wind turbines within a wind farm. This eliminates the need to introduce any superposition method. In the following, §2 develops the integral form of governing equations for turbine wake flows in wind farms. The large eddy simulation (LES) setup is described in §3. §4 derives the analytical solution of conservation of momentum deficit for a turbine within a wind farm. Model predictions are presented in §5. A modified version of the conservation of momentum deficit is proposed and solved in §6. In §7, we examine common wake superposition techniques in the literature. Finally, a summary and a discussion about possible future research directions are provided in §8.

Integral form of governing equations for turbine wakes within a wind farm
Let us assume a wind farm with an arbitrary layout of wind turbines (WT 1 ,WT 2 , ... WT , ... WT ) immersed in a turbulent boundary layer flow with a velocity profile denoted by 0 . The position of WT is denoted by X = ( , , ), where , and are the streamwise, spanwise, and vertical directions in the coordinate system, respectively. Turbines are labelled with respect to their streamwise positions such that −1 , where = {2, 3, ... }. The Reynolds-averaged Navier-Stokes (RANS) equation in the streamwise direction at high Reynolds numbers (neglecting viscosity effects) is given by (Shapiro et al. 2018(Shapiro et al. , 2019b) where , and are the time-averaged streamwise ( ), lateral ( ) and vertical ( ) velocity components, respectively. Turbulent velocity fluctuations are represented by , and and the overbar denotes time averaging. Also, is the time-averaged static pressure and is the air density. The term represents the effect of the thrust force of WT on the horizontal momentum and is given by where is the magnitude of the thrust force of WT in the streamwise direction, is the turbine radius, ( ) is the Dirac delta function, and ( ) is the Heaviside step function. Using the incoming boundary-layer profile 0 ( ), (2.1) can be written as From the continuity equation, we know Multiplying (2.4) by ( 0 − ) and adding the product to the left-hand side of (2.3) yields (2.5) Next, (2.5) is integrated from to with respect to , from to with respect to and from to with respect to , where ≪ < , ≪ ≪ and 0 ≪ < ≪ . Note that ≫ 0 to ensure that the assumption of negligible viscous forces is valid. The value of can be equal to zero only if the Reynolds shear stresses in (2.5) are replaced with the total shear stresses (i.e., sum of turbulent and viscous shear stresses). Without loss of generality, we assume that the integration box includes WT and an arbitrary number of upwind turbines as shown in figure 1. Integrating (2.5) yields where is a member of set if WT is inside the integration box. Also d is d d and d is d d d . In (2.6) and hereafter, any velocity or pressure term with a subscript denotes the value of the given variable in the presence of WT 1 ,WT 2 , ... WT . Now, we perform the same integration once more but this time in the absence of WT . This leads to where set ′ is equal to set excluding (i.e., \ { } = { : ∈ and ∉ { }}). As ≪ , surface integrals at = provide the same results in both (2.6) and (2.7). By subtracting (2.7) from (2.6), we obtain Note that the convective terms, which include crossstream velocity components, as well as the lateral Reynolds shear stress term in (2.5), vanish in equations written after (2.5), according to the fundamental theorem of calculus. For instance, , and therefore the difference in this term with and without WT is expected to be zero for ≪ ≪ (Tennekes & Lumley 1972;Pope 2000). However, the term including is kept due to the presence of the boundary-layer flow. LES data are used in the next section to quantify the value of each term in (2.8) for wind turbines in a wind farm.

Large eddy simulation (LES) setup
We use results from LES to compare against predictions of the analytical model presented in this paper. The LES is performed using the Simulator fOr Applications (SOWFA). SOWFA solves the incompressible filtered Navier-Stokes equations using a finite-volume formulation (Churchfield et al. 2012  to represent the shear stresses at the wall (Schumann 1975;Grötzbach 1987). The simulation is run for 20,000 [s] of simulated time for the turbulence to develop and then data on a boundary is sampled for 5,000 [s] to use in the simulation with the turbines as inflow. Figure 2 shows the spanwise averaged velocity and turbulence profiles as a function of height for the precursor simulation. All simulations presented use the same precursor simulation. A slight logarithmic layer mismatch is observed in the surface layer shown in figure 2b, which is common in LES of atmospheric boundary layers (ABLs) (Brasseur & Wei 2010).
The simulations with the turbines use the precursor as an inflow boundary condition and initial condition. The boundary conditions on the sides are periodic, and the outflow boundary condition adjusts the velocity field to conserve mass. The domain is a sub-sample of the precursor domain with 5 [km] by 1. A conventional variable-speed and variable blade-pitch-to-feather control system is used to control the turbine (Jonkman et al. 2009). The power and thrust coefficients for the turbine were determined by running a set of simulations with uniform inflow at different wind speeds. Figure 3 shows the thrust and power coefficients for the turbine model used in the simulations.
We simulate cases of aligned wind farm array with streamwise inter-turbine spacing of = 5 and spanwise inter-turbine spacing of = 3 , where is the rotor diameter. The inter-turbine spacing is intentionally chosen to be small to create strong turbine interactions within wind farms and thereby make it easier to examine capabilities and limitations of the developed analytical model. All simulated aligned wind farms consist of three turbine columns, but the number of turbine rows in simulations varies from one to five. For each case, the simulations are performed with and without the turbine located in the last row and the middle column, shown by the red colour in figure 4a (e.g., WT 6 for the 3 × 2 array or WT 15 for the 3 × 5 array). Removing turbines from the domain allows us to systematically study the magnitude of different terms in (2.8). We simulate another wind farm array with a slanted layout as shown in figure 4b. This wind farm, hereafter called slanted wind farm, is chosen to study model predictions in partial wake conditions. The slanted wind farm was constructed from the aligned wind farm by laterally shifting each row by 0.75 with respect to the upstream row. Unlike the aligned wind farm, the simulations for the   Figure 4: Schematic top view of the simulated (a) aligned wind farm cases with the number of rows varying from one to five, and (b) the slanted wind farm with five rows. For each case, simulations were performed with and without the last middle turbine shown by the red colour (e.g., WT 6 for the 3 × 2 aligned wind farm or WT 15 for the 3 × 5 slanted wind farm).
slanted wind farm are only performed for the 3 × 5 array with and without WT 15 shown by the red colour in figure 4b.

Derivation of wind farm analytical solution
The results from LES of the aligned wind farms are used to compute the integral quantities of (2.8) in a cubic box surrounding WT . The width of the box ( − ) is 500 [m], which is wide enough to ensure that the wake of WT is included. The vertical extent goes from = 20 [m] to = 300 [m]. The streamwise extent of the integration box is from = − 2 to = > . The first three terms on the right-hand side of (2.8) are surface integrals on a − plane at = . The Reynolds shear stress term is a surface integral over horizontal planes at = and , while the mean flow shear term is the only volume integral in the equation.
We note that the filtered variables from LES are similar to those derived in (2.8), however, there are some differences present, for example: -the LES velocities are spatially filtered quantities, -the unresolved part of the Reynolds stresses is neglected, -the pressure includes the deviatoric part of the Reynolds stress tensor. These differences are expected to cause a small difference in the momentum balance, which is captured in the residual term. See Sullivan et al. (1994);Juneja & Brasseur (1999); Ghate et al.
(2018) among others for a detailed discussion on the differences between pressure and velocity modelled by LES and their true values. Figure 5 shows the integral terms in (2.8) computed from the LES of the aligned wind farm cases with different numbers of rows. For each case, the results are shown for the last turbine in the middle column shown by the red colour in figure 4a (e.g., WT 6 for Row 2 ( = 6)). All terms are normalised by the value of the thrust force term. Figure 5 shows that the pressure term is the dominant term in the near wake, which is expected, as the sudden change in pressure at the rotor generates the turbine thrust force. The value of the momentum deficit term is clearly smaller than the thrust force in the near wake region. However, it increases with an increase of downstream distance and becomes comparable to the thrust force in the far wake region, especially for turbines in the first three rows. This is consistent with recent laboratory results reported by Hulsman et al. (2020). The next important term in the momentum equation is the Reynolds normal stress. This term is expected to be non-negligible for any wake flows, as also shown in prior studies of bluff-body wakes (Terra et al. 2017). Next, we examine the two terms that are introduced by the incoming boundary layer. The value of the Reynolds shear stress term is relatively small for turbines in the first and second rows. It seems that its value slightly increases for the one in the last row, although it is still smaller than other dominant terms. Note that values of the terms, including in (2.6) and −1 in (2.7), are expected to be significant due to the presence of the boundary layer (Cal et al. 2010;Bastankhah & Porté-Agel 2017). However, our LES data suggest that the difference in their values does not seem to be large, at least for a finite-size wind farm with five rows of wind turbines. The mean flow shear term is the last term on the right-hand side of (2.8). This term is the product of the vertical mean incoming flow shear and the vertical velocity component, mainly induced by the wake rotation. Figure 5 shows that this term is small but not negligible and its variation is relatively similar for turbines at different rows. We note that the effect of incoming wind veer is neglected in this analysis, given that its effect on the balance of momentum in the streamwise direction is expected to be insignificant. Indeed, veer becomes more important for the balance of momentum in the spanwise direction, which is not considered here. In the case of strong veer, the shape of the wake cross-section is skewed (Abkar et al. 2018), which is out of the scope of the current study.
It is important to note that for a very large wind farm that asymptotes to a so-called fully developed flow regime, changes in the streamwise direction can be almost neglected. In this case, the energy is mainly transported vertically from the top of the wind farm (Calaf et al. 2010; Abkar & Porté-Agel 2013; Meneveau 2012). However, even for this asymptotic case, the momentum deficit term in (2.8) is far from being negligible as this term is the difference of the momentum deficit flux, with and without the turbine at the same downwind position, not the difference of the one before and after a turbine. Future work is indeed needed to systematically examine the significance of different terms in (2.8) (especially the momentum deficit, pressure and Reynolds shear stress terms) for very large wind farms.
Based on the above discussion on the magnitude of different terms in (2.8), the below equation seems to be a reasonable approximation in the far wake of WT (at least for moderate values of ): (4.1) Figure 5: Variation of different terms in (2.8), normalised by the thrust force term, for the last turbine in the middle column of aligned wind farms with different number of rows. The modified momentum deficit term is discussed later in §6.
Equation (4.1) can be loosely named as conservation of momentum deficit. It states that the presence of a turbine induces a rise in the value of momentum deficit flux, and the magnitude of this rise is equal to the turbine thrust force. It is important to bear in mind that the conservation of momentum deficit is not an intrinsic flow governing equation like those for mass and momentum that should always hold true. As shown earlier, (4.1) is just an approximate relation deduced from conservation of mass and momentum (2.8). Also note that (4.1) is not the same as the one used for single turbine wakes (e.g., Bastankhah & Porté-Agel 2014). See §7 for further discussion.
In the following, we aim to determine by solving (4.1). To achieve this goal, we need to use the assumption of self-similarity for velocity-deficit profiles in turbine wakes. Unlike single isolated turbines, the definition of velocity deficit with respect to the incoming flow (i.e., , ) is not suitable for turbines within wind farms. This definition can even lead to negative values of velocity deficit at the centre of the wake for very large downwind distances, where the wake velocity becomes greater than the incoming flow (i.e., as ( − ) → ∞, → 0 , ). Instead, we define the velocity deficit at a given position X = ( , , ) downwind of WT as the difference of the streamwise velocity in the absence and presence of WT at X; i.e., (4.2) Figure 6 shows that with this definition of velocity deficit, the wake of a turbine in a wind farm exhibits a good degree of self-similarity, akin to a stand-alone turbine. As seen in the figure, velocity-deficit profiles collapse to a single curve for different downwind positions if the velocity deficit is normalised by the maximum velocity deficit and the distance from the wake centre is normalised by the characteristic wake half width . The results are shown in figure 6 for wakes of three turbines: (a) a stand-alone wind turbine, (b) the middle turbine in the last row of the aligned wind farm (i.e., 15 ), and (c) the middle turbine in the last row of the slanted wind farm (i.e., 15 ). It is also worth mentioning that a slight lateral wake deflection observed in the figure for the LES data is likely due to the interaction of rotating wake with the vertical inflow shear discussed by prior studies (e.g., Fleming et al. 2014;Gebraad et al. 2014). The developed analytical model in this paper does not take into account this slight wake deflection for zero-yawed

turbines. As
where is the self-similar function. Shifting the index in (4.3) to − 1, − 2, ... 1 leads to a set of equations as follows: (4.5) Using (4.3) and (4.5) to rearrange (4.1), we obtain , which leads to a slight flow speed-up, especially at short downwind distances. This flow speed-up can cause lower than expected or even negative values of velocity deficit at the wake edges for a stand-alone turbine as shown in figure 6a. In a wind farm, a more pronounced flow speed-up may occur between adjacent turbine columns due to flow blockage effects. For instance, this is the case for the simulated aligned wind farm as shown in figure 6b.
On the other hand, our results suggest that the flow speed-up for the last turbine in the slanted wind farm (figure 6c) seems to be less significant. The magnitude of flow speed-up around wind turbines may depend on several factors such as turbine and inflow properties as well as wind farm layout geometries (Garrett & Cummins 2007;Nishino & Draper 2015). The accurate estimation of flow speed-up is out of the scope of this study and so we use a Gaussian distribution for its simplicity. However, we acknowledge that this assumption can introduce errors in flow prediction at wake edges. See §5 for more discussion. With an assumption of the Gaussian distribution for the wake velocity deficit, Note that although the focus of this paper is turbines with no yaw angles, the developed model can be extended to cases with yawed turbines by replacing and in (4.7) with lateral and vertical positions of the wake centre at each streamwise position, respectively.
Inserting (4.7) into (4.6), computing the surface integral with respect to and (from −∞ to +∞, neglecting ground effects) and approximating 0 with (4.9) Note that for simplicity we assume in this work that the wake width is the same in both lateral and vertical directions. However, different values of lateral ( ) and vertical ( ) wake half widths can be used in (4.7). In this case, any 2 in (4.8) and hereafter should be replaced with the product of and . The quadratic (4.8) has two solutions for . The physically acceptable solution that decays with an increase of is where <> ( , ) denotes spatial averaging over the frontal projected area of WT at = , and the thrust coefficient of WT is given by , ) . (4.11) Values of 1 , 2 , ... determined from (4.10) are inserted into (4.5) to evaluate . While we use < −1 > , to relate , to in (4.11), one may approximate < −1 > , with −1 ( , , ) in (4.10) and (4.11) for simplicity. Note that (4.10) is a recursive sequence, so can be explicitly computed as a function of (for = 1, ... − 1). To compute 1 (i.e., = 1), by setting = 0, (4.10) is reduced to the one developed by Bastankhah & Porté-Agel (2014) for a stand-alone turbine.
The dimensionless coefficient of in (4.10) has an interesting physical interpretation. It quantifies the contribution of WT on the value of (i.e., the velocity deficit associated with WT ). From (4.9), depends on the wind farm layout and inflow conditions, and its value can vary from 0 to 2. As | − | tends to infinity, tends to zero. This means that if the turbines are laterally distant from one another, they do not have any interaction, so (4.10) is reduced to the solution derived for a single turbine. As approaches (i.e, partial-wake conditions), the value of increases. Ultimately, at = (i.e., full-wake conditions) and assuming = , becomes equal to 2 2 /( 2 + 2 ), which can take any value between one and two. In these conditions, tends to 2 if goes to infinity, which occurs when the turbine is immersed in a reasonably large upwind wake. On the other hand, the value of tends to unity if the upwind turbine wake has a size comparable to that of the turbine wake (i.e., ≈ ). The inflow properties also affect the value of . An increase in the level of atmospheric turbulence enhances the wake recovery rate, which in turn leads to an increase in the value of . This occurs particularly for turbines in front rows. Therefore, the coefficient obtained purely based on an analytical approach is analogous to empirical methods used in the literature to quantify the effects of upwind turbine wakes, such as finding the areas of wake overlap with each turbine. For instance, see the so-called 'mosaic-tile' approach used by Rathmann et al. (2006), among others.
The second exponential term in the right-hand side of (4.9) is equal to unity for = . This is the common case in wind farms given that the turbines usually have the same hub height. We, however, leave this term in its original form because of (i) its potential use in imaging techniques to simulate ground effects (Crespo et al. 1999) as well as (ii) its potential application in studying wind farms with variable hub heights (i.e., vertical staggering) (Zhang et al. 2019).
It is also worth mentioning that in the case of 0 = 0 ( , , ), ℎ in (4.10) is substituted with 0 ( , , ). However, note that the developed solution is expected to provide acceptable estimations as long as values of d 0 /d and d 0 /d are not significant. Otherwise, the incoming flow heterogeneity induces non-negligible additional terms in (2.8), akin to the vertical mean flow shear term.

Model predictions
In this section, the developed model is used to predict the flow distribution in each of the two wind farms shown in figure 4, and the results are compared with the LES data. In order to compute predictions using the analytical model, an estimation of the wake recovery rate is required as the only empirical input of the model. Prior studies have suggested that the wake recovery rate for a single turbine is directly proportional to the turbulence intensity of the incoming flow (Niayifar & Porté-Agel 2016;Carbajo Fuertes et al. 2018;Shapiro et al. 2019b;Zhan et al. 2020). In the following, we examine the validity of this assumption for turbines of the aligned wind farm, for which LES were performed in the presence and absence of turbines at different rows.
The lateral velocity deficit profiles at different downstream distances of the middle turbine in each row are analysed using the definition of velocity deficit stated earlier as Δ (X) = −1 (X) − (X) for WT (i.e., the difference between the flow with and without the presence of WT ). A Gaussian distribution (4.7) is fitted to each velocity deficit profile in the horizontal plane at hub height in order to estimate the corresponding wake half width for each downstream position. Figure 7a shows the variation of wake half width with downwind distance for the aligned wind farms with different number of rows. Results show that the wake width is clearly smaller for the turbine in the first row, which is expected because of the lower level of turbulence intensity in the incoming flow. Most of the turbine wakes appear to have a linear expansion. The results for the third and fifth rows appear slightly less linear, which may be due to some uncertainty in the estimation of the wake width using a Gaussian curve fitting as discussed in prior studies (e.g., Quon et al. 2020). However, a linear curve still provides satisfactory agreement. Subsequently, the value of for WT is determined by fitting a linear curve with where is the normalised initial wake half width given by 0.2 √ and = 1 + 1 − , ) / 2 1 − , (Bastankhah & Porté-Agel 2014). Fitting the curve in this way ensures that the initial wake width is consistent throughout the model to give values of that are a fair representation of the wake recovery rate.
Following this, the proportionality of in relation to the incoming turbulence intensity is investigated. For , the value of is defined as the value of turbulence intensity at the rotor position in the absence of the rotor, i.e., < > ( −1, ) /< > (0, ) . Figure 7b shows that / is not exactly identical for all turbines. However, the deviation of this ratio between different rows is not significant, and so the mean value of / = 0.31 can be taken to approximate the relationship between the wake recovery rate and the incoming turbulence intensity . This value is lower than those suggested in previous works, such as 0.  (2020). The discrepancy in the coefficient that relates to may be due to the wake recovery rate having some weak dependencies on incoming flow characteristics other than the incoming turbulence intensity (e.g., the integral length scale). It is not in the scope of this work to establish a complete relationship between the wake recovery rate and incoming turbulence intensity, as well as potentially other important factors. Nevertheless, this is an important topic and needs to be studied more rigorously in future research. This is further discussed in §8.
The final step in estimating the wake recovery rate for the analytical model requires an estimation of the incoming turbulence intensity for each turbine in a wind farm. The empirical relationship suggested by Crespo & Hernández (1996) is used here for simplicity, although other, more detailed, empirical relations could be used instead, such as the one proposed by Ishihara & Qian (2018). The magnitude of incoming turbulence intensity for each turbine is taken as 0 2 + Δ 2 , where 0 is the ambient turbulence intensity, and Δ is the turbulence intensity added by upwind turbines . The value of Δ for WT due to the upwind turbine WT , where < , is estimated with the relationship proposed by Crespo & Hernández (1996) (2016). This involves finding the fraction of overlap area between the turbine rotor and upstream wakes to determine the added turbulence intensity due to each upstream turbine, followed by taking the maximum value of Δ (i.e., the upstream turbine with the largest impact on the incoming flow). By using this method, predictions show that turbines in the second row and subsequent rows experience the same incoming turbulence intensity. However, this is not completely valid for this wind farm case as the LES data show that the value of incoming turbulence intensity for the second row is less than that of subsequent rows. Additionally, this method slightly overestimates for turbines deep inside a wind farm. In order to account for this offset, the Crespo model is employed for the analytical model with a slightly modified constant term to obtain the relationship Δ = 0.66 0.83 0.03 0 [( − )/ ] −0.32 . As demonstrated in Figure 7c, this modified version gives closer predictions of relative to the LES data. Next, the wind farm flow distribution is estimated using (4.5), (4.9) and (4.10). Note that the model can only provide reliable predictions in the far-wake region of wind turbines, where velocity-deficit profiles are self-similar, and also the assumptions made to develop the approximate relation of (4.1) are acceptable. At very short downwind distances, the term under the square root in (4.10) becomes negative and, as a result of this, (4.10) provides complex values. Prediction of the flow in the near wake region is not the objective of this work because it is unlikely for a turbine to be in the near wake of another turbine in realistic situations. However, the complex output of (4.10) might pose a practical issue in implementing the model to determine the incoming velocity for downwind turbines in some cases, where two turbines are laterally far from each other but have similar streamwise locations. In order to address this issue, any value of predicted by (4.10), which is either complex or bigger than 0 is replaced by 0 , where 0 is the maximum theoretical velocity deficit based on the Betz theory and is given by et al. 2010). We adopt this approach for its simplicity, but if the goal is to realistically capture the near-wake region, recent models in the literature (e.g., Shapiro et al. 2019b;Schreiber et al. 2020; Blondel & Cathelain 2020, among others) can be consulted.
Finally, we discuss model predictions against the LES data. Figures 8a and b show contours of ( ℎ − ), normalised by the inflow velocity at hub height ℎ for the 3 × 5 aligned wind farm (i.e., the full-wake case). A different colourscale is used for negative values of ( ℎ − ) to indicate the speed-up region between adjacent turbine columns. This speed-up region appears to be evident between the primary rows of the wind farm, diminishing after roughly three rows of wind turbines. As mentioned earlier in §4, this speed-up region is not accounted for in the developed model due to the assumption that the wake velocity deficit profile is Gaussian. For the aligned wind farm, the speed-up region does not largely interact with downwind turbines, so its impact is expected to be insignificant. However, this is not always the case, as discussed later for the case of the slanted wind farm. Overall, figure 8b shows that the model predictions of velocity in the far-wake region are in good agreement with the LES data. This is also confirmed in figure 8c showing vertical profiles of normalised streamwise velocity at different positions downwind of the middle turbine in the last row (i.e., WT 15 ). Although the results show good agreement for the far-wake region, the figure shows that the model slightly underpredicts the velocity in the lower half of the wake, which becomes more apparent at short downwind distances (e.g., at − 15 = 4 ). This discrepancy may be due to the uncertainty in the estimation of the wake recovery rate, as discussed earlier, or due to other terms being neglected in the right-hand side of the momentum equation (2.8). For cases where the momentum deficit term is less than the thrust force, assuming equality between these two terms will result in an overestimate of velocity deficit, which occurs particularly at short downwind distances. This is further discussed in §6. Additionally, it was assumed for simplicity that the wake recovery rate is the same in both the lateral and vertical directions. However, in the vertical direction, the wake flow may be affected by the mean shear of the incoming flow and presence of the ground resulting in a different value of , as suggested by prior studies (Abkar & Porté-Agel 2015;Xie & Archer 2015).
Finally, figure 8d shows the efficiency of each middle turbine at different rows for both the LES data and predictions from the analytical solution. The efficiency of WT is defined as /(0.5 3 ℎ ), where is the power extracted by WT and is the area swept by the turbine blades. The turbine efficiency can be rewritten as < > 3 ( −1, ) / 3 ℎ , where is the power coefficient of the turbine. For analytical model predictions, the value of < > ( −1, ) is estimated by the model, whereas for each turbine is estimated from figure 3. Overall, there is a good agreement between the efficiency predicted by the analytical model and the LES data, despite a slight underprediction of the efficiency for turbines in the last two rows. One possible explanation for this underprediction is the assumption of equality in (4.1), as discussed previously. As shown in figure 5, this assumption is less accurate for turbines at the end of the wind farm, compared to those in primary rows.
For the case of the slanted wind farm, figures 9a and b show contours of ( ℎ − ), normalised by ℎ , in the horizontal plane at hub height for LES data and the analytical model predictions, respectively. From figure 9c, model predictions show fairly good agreement with the LES data in the far wake of the turbine in the last row (i.e., WT 15 ), but it underpredicts the velocity at short downwind distances. Figure 9d shows the comparison between predictions of turbine efficiency from the analytical model and the LES. Model predictions for turbines in rows 4 and 5 show more satisfactory agreement than those for rows 2 and 3. The observed discrepancy in efficiency for rows 2 and 3 are most likely due to the flow speed-up demonstrated in figure 9a. As discussed in §4, the analytical model does not capture flow speed-up between adjacent turbine columns due to the assumption of Gaussian distribution for velocity-deficit profiles. This leads to underprediction of efficiency for turbines in the second and third rows, as shown in figure 9d. The turbines in the last two rows do not seem to considerably benefit from the flow speed-up because they are in wakes of adjacent column. A similar observation is made in figure 6c showing self-similar profiles for the turbine in the last row.
These results highlight the limitation of this model (and essentially most other existing engineering wake models in the literature) for estimating the incoming flow of turbines that are subject to the speed-up between neighbouring columns. However, our LES data suggest that this speed-up effect appears to be important mostly at primary rows of wind farms and becomes less significant deep inside a wind farm. Moreover, the speed-up effect is expected to have less of an impact for wind farms with larger inter-turbine spacing. An important area of future research would be to evaluate the impact of flow speed-up on the power production for wind farms with various inflow conditions and layout configurations.

A modified version of the conservation of momentum deficit
As shown in figure 5, the magnitude of the momentum deficit term in (2.8) is clearly less than the thrust force in the near wake for all turbines. In the far wake, its value is closer to the turbine's thrust force. However, for turbines deep inside the wind farm (i.e., rows 4 and 5), the momentum deficit term is still slightly smaller than the thrust force even in the far wake region. This section attemps to modify the momentum deficit term such that this term's value would be more comparable to the thrust force magnitude over a broader range of streamwise distances and turbine rows. We start with rewriting the left-hand side of (4.1), so we obtain ∫ The second term on the right-hand side of (6.1) is negative as < −1 and −1 0 . Therefore, one can write ∫ The equality in (6.2) holds everywhere for = 1. For > 1, the equality holds only at large downwind distances where → −1 → 0 , while the difference between the two sides of inequality is larger at shorter downwind distances. Figure 5 shows the variation of the modified momentum deficit term with the streamwise distance for turbines at different rows. While the value of this modified term is bigger than the thrust force in the far wake of turbines in primary rows (i.e., rows 2 and 3), its value in the far wake of those in rows 4 and 5 is closer to the thrust force than the original relation. Moreover, the modified term is closer to the thrust force at short downwind distances for all turbines. Therefore, we introduce a modified version of the conservation of momentum deficit as follows: ∫ Based on the LES data presented in figure 5, this equation is expected to better hold the equality between the thrust and momentum deficit term at short downwind distances for all turbines and at long downwind distances for those deep inside a wind farm. Equation (6.3) can be solved in the same way that we earlier solved (4.1) in §4 to find . By doing so, we find that the solution of (6.3) is the same as (4.10). The only difference with the original solution is that , based on the modified version, is half of the one obtained for the original equation. Therefore, using (4.10) in conjunction with (4.5) and the modified definition of (6.4) forms the solution of the modified version of conservation of momentum deficit (6.3). Figure 10a shows the variation of the centreline ( ℎ − )/ ℎ for turbines in the middle column of the aligned wind farm, while figure 10b shows the efficiency of these turbines. It is worth remembering that model predictions in the upwind induction region and immediately behind the turbines are not valid. The figure shows that far-wake predictions, based on the original and modified versions, do not largely deviate from each other (e.g., compare flow predictions downwind of the last turbine at ( − 15 > 5 ). This is a good sign showing that model predictions are not highly sensitive to inaccuracies arising from using the approximate form of (2.8). However, a closer inspection of figure 10 confirms what we inferred earlier based on the results of figure 5. The original momentum deficit works better in the far wake of turbines at primary rows (e.g., see for row 3 in figure 10b). On the other hand, the modified version provides better predictions at short downwind distances for all turbines. It can also better predict the efficiency of turbines deep inside a wind farm (e.g., see for row 5 in figure 10b). This comparison suggests that, overall, the modified version of the conservation of momentum deficit might be a slightly better choice, at least for turbines deep inside a wind farm. Future studies can perform similar analyses using more numerical and experimental data for larger wind farms to examine the accuracy of this conclusion.

Superposition of single turbine wake (STW) models
In this section, we examine the validity of wake superposition techniques commonly used in the literature. To achieve this goal, we discuss approximations and assumptions that need to be made in order to derive superposition of STW models from the solution of conservation of momentum deficit derived in §4 for a wind turbine in a wind farm.
As discussed in §1, the common approach in existing superposition techniques is to treat each turbine separately and find the value of its wake velocity deficit at a given location. The values of velocity deficit caused by all turbines at that given location are then combined (either linearly or non-linearly) to find the total velocity deficit. By doing this, it is inadvertently assumed that the wind velocity downwind of a turbine subtracted from the one at the same location in the absence of the turbine can be expressed in a form similar to (4.3). As shown in figure 6, this is a valid assumption for a turbine within a wind farm. Next, we try to develop the STW model from the modified version of the conservation of momentum deficit (6.3). If we approximate −1 in this equation with < −1 > ( , ) at each , (4.3), (4.7) and (4.11) can be used to obtain ) . (7.1) Next, we neglect the velocity ratio under the square root on the right-hand side of (7.1). Moreover, < −1 > ( , ) (i.e., the first term on the right-hand side of (7.1)) is substituted with ℎ . This leads to which is the original STW model proposed by Bastankhah & Porté-Agel (2014). The velocity ratio of < −1 > ( , ) /< −1 > ( , ) tends to unity only as ( − ) → 0 and is less than one at positive values of ( − ) in most cases. On the other hand, the velocity ratio of < −1 > ( , ) / ℎ tends to unity only as ( − ) → ∞ and is less than one at definite values of ( − ). Therefore, the two approximations made to develop (7.2) from (7.1) lead to an overprediction of velocity deficit. Flow predictions based on (7.2) and the linear superposition method (i.e., A.I in table 1) are shown in figure 10. One can observe that using this STW superposition model clearly leads to overpredictions of velocity deficit (see figure 10a) and, consequently, underprediction of turbine efficiency (see figure 10b). Similar observations were made by prior studies (e.g., Niayifar & Porté-Agel 2016; Zong & Porté-Agel 2020). In (7.2), is proportional to 0 , which is equivalent to the superposition method B.I shown in table 1. To develop the method B.II in table 1, we follow the same approach to derive (7.2) from (7.1), except < −1 > ( , ) in (7.1) is now replaced with < −1 > ( , ) (i.e., local incoming velocity), so we obtain As often < −1 > ( , ) < −1 > ( , ) ℎ for , the magnitude of velocity-deficit overprediction from (7.3) is significantly less than that from (7.2), which is also seen in figure 10.
In comparison with the solution of the original conservation of momentum deficit (4.10), (7.3) underpredicts the velocity deficit at short downwind distances and overpredicts it at large ones, but overall its predictions do not largely deviate from those of (4.10). However, it is important to note that this is not because (7.3) is based on flow physics for the wake of a turbine within a wind farm. Several approximations made to develop (7.3) from the original solution have opposing effects on the prediction of wake velocity deficit, so their effects are cancelled out by each other to some extent. Figure 10b shows that using (7.3) in conjunction with the linear superposition method induces an error of about 3% − 4% in predictions of turbine efficiency for the studied wind farm. Zong & Porté-Agel (2020) stated that the error of this superposition model might be more significant for larger wind farms. Therefore, the use of the analytical wind farm solution (4.10) is recommended because it is not computationally more expensive than empirical STW superposition models. Nevertheless, if one tends to use superposition of STW models, (7.3) certainly provides more acceptable estimations than (7.2). Another important aspect of superposition techniques is the method used to superpose the velocity deficit caused by each turbine. Bearing in mind (4.5), it is evident that the linear superposition method (i.e., A.I in table 1) is in better agreement with the analytical solution. However, the root-sum-square method (A.II in table 1) is more likely to provide better predictions if the velocity deficit caused by each single turbine is overestimated. Overestimation of velocity deficit for turbine wakes in wind farms mainly occurs for two reasons: (i) using the global incoming velocity ℎ as the reference velocity to compute the velocity deficit for each turbine (7.2), and (ii) estimating the wake recovery rate based only on inflow atmospheric conditions. The latter is expected to underestimate the wake recovery rate and thereby overpredict the velocity deficit, because it does not take into account a faster wake recovery due to the turbulence added by upwind turbines. This may explain why the root-sum-square method has been customarily popular in the literature as either one or both of the above-mentioned assumptions can be commonly found in prior wake modelling studies (e.g, Katić et al. 1987, among others).
Finally, we discuss the recent model developed by Zong & Porté-Agel (2020). This model is based more on wake flow physics than other common superposition methods discussed earlier.
They developed an iterative method to satisfy the following two equations: (7.5) While (7.5) satisfies the conservation of momentum deficit for the whole wind farm, strictly speaking, (7.4) is not equal to the conservation of momentum deficit (4.1) for a turbine within a wind farm. In fact, one can show that The second term on the right-hand side of (7.6) is not negligible for a turbine that experiences wakes of upwind turbines. However, this discrepancy does not seem to cause a noteworthy error in model predictions in the far wake. Figure 11 shows variation of the centreline (1 − 2 / ℎ ) with downwind distance for a turbine that is subject to the full wake of another turbine. For a fair comparison, the same value of obtained from figure 7 is used for all results shown in figure 11. The figure shows that Zong & Porté-Agel (2020)'s model provides good far-wake predictions, which lie between results of the original and modified versions of conservation of momentum deficit. However, a drawback of this model could be the fact that it is not represented in an explicit form. To use this model, one needs to compute surface integrals of wake velocity deficit multiple times at each downstream location in order to obtain converged results through an iterative process. This may increase the computational cost of this superposition model compared with simpler superposition models as well as the explicit wind farm solution developed in the current study.

Summary and future research
The main aim of this paper is to address the puzzling question of 'how should one estimate cumulative wake effects in wind farms based on engineering wake models?', given the abundance of wake superposition methods in the literature. This aim is achieved by directly solving flow governing equations for a turbine that experiences upwind turbine wakes. The developed model can therefore predict the flow distribution in a wind farm without the need for any specific superposition method. The LES data for wind farms with different number of rows are used to perform a budget analysis of the integral form of the mass and momentum equations for turbine wakes in wind farms. Results show that there are important non-negligible terms (e.g., pressure and Reynolds normal and shear stress terms) other than the momentum deficit and thrust force in this equation. However, the so-called conservation of momentum deficit (i.e., the equality of the momentum deficit term and thrust force) seems to be fairly valid in the far wake, at least for a moderately sized wind farm. A modified version of the conservation of momentum deficit is also introduced in an attempt to provide slightly better results at short downwind distances as well as in the far wake of turbines deep inside a wind farm. Performing LES with and without the presence of turbines allows us to properly examine some model assumptions. The data for both full-wake and partial-wake conditions show that wake velocity-deficit profiles are self-similar if they are defined with respect to the wind flow at the same position but in the absence of the turbine. While a Gaussian profile can successfully represent the self-similar profile of the wake at the centre, it falls short in capturing the flow speed-up at wake edges. This results in an underprediction of efficiency for certain turbines in primary rows that experience flow speed-up between adjacent columns.
To quantify the wind farm flow distribution using this analytical model, the only necessary empirical input is the wake recovery rate. Provided that this is properly estimated for turbines within wind farms, our results show, overall, that the proposed model is able to provide acceptable predictions. Based on this model, the velocity distribution downwind of a wind farm consisting of wind turbines (WT 1 ,WT 2 , ... WT ) is given by where the wake half width for WT is estimated from (5.1). The value of (i.e., wake-centre velocity deficit associated with WT ) is determined by where , is the thrust coefficient of WT . For = 1, the equation is reduced to the one for the single turbine wake by setting / ℎ equal to zero. For > 1, the value of in (8.2) is given by where is equal to 2 for the original equation of conservation of momentum deficit (4.1) and is equal to 1 for its modified version (6.4). The value of depends on wind farm turbine layout and inflow conditions. This coefficient, obtained analytically, functions similarly to empirical models (e.g., the 'mosaic-tile' approach) that aim to quantify the net contribution of overlapping wakes. The validity of the common single turbine wake (STW) superposition models is also analysed in this study. It is shown that common superposition models can be derived by making approximations to the analytical wind farm flow solution presented in this work. There are still some important limitations that need to be better understood and addressed in future research if we want to successfully implement engineering wake models for a variety of inflow conditions and wind farm layout configurations: -Neglected terms of the governing equation (2.8): The applicability of the so-called conservation of momentum deficit (4.1) (and its modified version (6.3)) should be examined for wind farm arrays with various inflow conditions and layout geometries. Of special interest is to investigate the significance of different terms in (2.8) for a very large wind farm that asymptotes to a so-called fully developed case.
-Turbine blockage effects: One of the limitations of the proposed model is its inability to capture speed-up effects between adjacent turbine columns caused by turbine flow blockage. Future research can aim at using a more realistic relation (instead of Gaussian) to represent the self-similar velocity-deficit profiles. Furthermore, several recent studies (Bleeg et al. 2018;Segalini & Dahlberg 2020) have demonstrated another important aspect of flow blockage in wind farms. These studies, among others, have shown that flow blockage caused by wind turbines within wind farms can decrease the efficiency of upstream turbines by reducing their effective incoming velocity. For more accurate prediction of power production in large wind farms, the effect of downwind turbines on their upwind counterparts needs to be included in wind farm flow engineering models.
-Wake recovery rate : Model predictions are quite sensitive to the value of . Therefore, accurate estimation and modelling of is of great importance. To achieve this goal, a better understanding of the turbulence distribution in wind farms and how this affects the turbine wake recovery rate is essential. Moreover, additional research is needed on the possible effects of any inflow and turbine operating conditions on wake recovery, other than the incoming turbulence intensity.
Declaration of Interests. The authors report no conflict of interest.