Design and analysis of a wake model for spatially heterogeneous flow

Methods of turbine wake modeling are being developed to more accurately account for spatially variant atmospheric conditions within wind farms. Most current wake modeling utilities are designed to apply a uniform flow field to the entire domain of a wind farm. When this method is used, the accuracy of power prediction and wind farm controls can be compromised depending on the flow-field characteristics of a particular area. In an effort to improve strategies of wind farm wake modeling and power prediction, FLOw Redirection and Induction in Steady State (FLORIS) was developed to implement sophisticated methods of atmospheric characterization and power output calculation. In this paper, we describe an adapted FLORIS model that features spatial heterogeneity in flow-field characterization. This model approximates an observed flow field by interpolating from a set of atmospheric measurements that represent local weather conditions. The objective of this method is to capture heterogeneous atmospheric effects caused by site-specific terrain features, without explicitly modeling the geometry of the wind farm terrain. The implemented adaptations were validated by comparing the simulated power predictions generated from FLORIS to the actual recorded wind farm output from the supervisory control and data acquisition (SCADA) recordings and large eddy simulations (LESs). When comparing the performance of the proposed heterogeneous model to homogeneous FLORIS simulations, the results show a 14.6 % decrease for mean absolute error (MAE) in wind farm power output predictions for cases using wind farm SCADA data and a 18.9 % decrease in LES case studies. The results of these studies also indicate that the efficacy of the proposed modeling techniques may vary with differing site-specific operational conditions. This work quantifies the accuracy of wind plant power predictions under heterogeneous flow conditions and establishes best practices for atmospheric surveying for wake modeling.


Introduction
Low-fidelity wake modeling utilities such as FLOw Redirection and Induction in Steady State (FLORIS) are typically used for the estimation of wind farm power output or the implementation of wind farm controls that help improve the overall performance of a wind farm. This includes implementing real-time corrective strategies that aid in reducing stress-inducing loads on turbines (Boersma et al., 2017), avoiding operational side effects like noise pollution (Leloudas et al., 2007) or shadow flicker (Clarke, 1991), and maximizing power output through methods of wake steering and power grid optimization (Fleming et al., 2017b). FLORIS, and most other control-oriented wake modeling utilities, implement advanced wake modeling algorithms that are capable of producing accurate results in a uniform set of atmospheric conditions . However, the accuracy of any wake model is highly dependent on its ability to recreate the characteristics present. It is important for these models to be able to emulate the naturally occurring state of the wind farm as closely as possible for the controls processes and power prediction functionalities to operate with reliable accuracy. Most current control-oriented wake modeling utilities use a homogeneous approximation to characterize the initial state of the atmosphere, which can introduce major inaccuracies in the simulation of wind-farm-flow interactions.
The consequences are particularly evident when observing the accuracy of power predictions for wind farms located within complex terrain or wind farms that are otherwise subject to spatially variant conditions in the atmosphere. Because these atmospheres are subject to dramatic changes in the velocity and direction of wind, it is difficult to anticipate how the resulting wakes will form and what kind of power output should be expected. In Yang et al. (2019), an analysis of the impact of spatial heterogeneity in wind farm flow is presented for a site within complex terrain. This study showed that using averaged values of wind conditions caused short-term wind power forecasting to be less accurate, due to spatial heterogeneity within the wind field and the variability of wind turbine power curves. With these effects considered, the current version of FLORIS and many other wake model utilities are not constructed to accurately model fluid flow under these conditions.
It should be noted that there are existing wake models that incorporate elements of heterogeneous wake effects caused by varying atmospheric conditions. For example, one model presented in You et al. (2016) takes a statistical approach in representing heterogeneous power deficit caused by windfarm-flow interactions in spatially variant weather conditions. Another method discussed in Shao et al. (2019) proposes an interaction model used for calculating the turbulence intensity of overlapping wakes and represents the relative positions of wind turbines under arbitrary and varying wind direction conditions. Brogna et al. (2020) presents a technique that superimposes the centerlines of wind turbine wakes in complex terrain by following the streamlines of the background flow field. Clustering methods have also been implemented, such as Katic et al. (1986) and Clifton and Lundquist (2012), where the turbines of a wind farm are sectioned into groups, assigning differing atmospheric characteristics to each cluster of turbines to mimic the heterogeneous conditions observed in natural atmospheres. Additionally, many approaches implement data-driven wake model correction parameters to achieve more accurate solutions, such as those proposed by Schreiber et al. (2019), Shapiro et al. (2019), Howland et al. (2020), and Teng and Markfort (2020).
The aforementioned models present many methods for approximating farm-flow interaction in heterogeneous conditions. As a contribution to this area of research, this article will present a modified version of FLORIS that features an advantageous capability in modeling wind farms with spatially variant weather conditions and complex terrain. This adapted version of FLORIS presents several novel developments within the scope of control-oriented wake modeling research: an interpolation algorithm is implemented, which allows the user to define a gradient of atmospheric characteristics across the flow field, based on several measurements within or adjacent to the wind farm; elements of spatially variant wind direction, wind speed, and turbulence intensity are integrated into wake calculations of the preexisting FLORIS model; and an additional method is introduced to minimize error in power prediction accuracy caused by highturbulence intensity and wind speed variance.
The objective in developing this proposed model is to capture a more accurate representation of the effects of wind farm wake interactions within complex terrain without actually resolving any terrain geometry during simulation. This study aims to analyze the accuracy of power output predictions and wake modeling performance for the proposed wake model, through comparisons to large eddy simulation (LES) wind farm supervisory control and data acquisition (SCADA) records.

Existing FLORIS model
FLORIS (NREL, 2020) is a wake modeling utility that is equipped with tools designed for the control and optimization of wind farms and is being developed at the National Renewable Energy Laboratory (NREL) in collaboration with Delft University of Technology. This tool uses several computational modeling techniques paired with controls algorithms to approximate and optimize wind turbine wake interactions through integration of real-time supervisory control and data acquisition (SCADA) data recorded from wind farms. FLORIS implements the concept of steady-state averaging to simulate the observed dynamic behavior within a wind farm for each iteration in time and can also be used as a simulation tool to compute farm-flow interactions in wind farms under user-defined atmospheric conditions. This section will give an overview of the mathematical theory in which the formulations of the wake models of FLORIS were based. These concepts are also explained in greater detail in Annoni et al. (2018) and Hamilton et al. (2020).

Turbine power output model
The operation and performance of a turbine is modeled with respect to the relationship between the thrust coefficient, C T , and power coefficient, C P . The dependence between these two terms characterizes a turbine's power output and wake propagation, therefore making the understanding of this relationship fundamental to the design and operation of wind farm controls. To model the performance behaviors of a given turbine, a table is constructed inside of FLORIS that tabulates C T and C P with respect to wind speed. This table can be set to a user's self-obtained data, generated independently by NREL's FAST (Jonkman and Buhl, 2005) or by integrating CCBlade (Ning, 2013) with FLORIS. The relationship between C T and C P can also be defined through the concept of actuator disk theory. This theory relates the turbine power output and thrust through the axial induction factor, a, which can be calculated using the definitions from Burton et al. (2002) and Bastankhah and Porté-Agel (2016): (1) (2) From these values, the power can then be calculated for turbines under steady-state and yaw-misalignment conditions, using the following relationship provided by Burton et al. (2002): where ρ is the air density, A is the rotor-swept area, u is the rotor-averaged wind speed, and p is a tuneable parameter that accounts for the power losses due to yaw misalignment seen in simulations (Burton et al., 2002;Fleming et al., 2017a). Thus far, the turbine model discussed in this section does not consider the effects that turbulence may have on the relationship between power output and wind speed. However, Sheinman and Rosen (1992) analyze the effects of turbulence intensity on wind farm power output. In this study, it is shown that turbine power output can be overestimated by more than 10 % if turbulence intensity is not considered. Many empirical and machine-learning methods have been proposed to solve this issue. However, a nonparametric statistical averaging model may be preferred, such as the model developed in Hedevang (2014). In Sect. 3.5, a new method of implementing a turbulence-dependent correction to power will be discussed for FLORIS applications.

Velocity deficit
FLORIS provides an option to select particular models for wake velocity deficit and wake deflection separately to suit the user's performance needs. The variety in modeling capabilities reflects a range of trade-offs between computational efficiency and the number of detailed physics applications applied to calculations. If a model is more computationally expensive, it is likely to implement more sophisticated algorithms as well, in hopes of achieving a more accurate result. These models all have a different approach to modeling turbine wake interactions and offer different strengths and weaknesses in functionality. Most models can be classified as either a velocity deficit or a wake deflection calculation, but there are also the Gaussian and Curl models that incorporate both calculations and extend further into the overall FLORIS wake modeling structure and control tools. For the purposes of this article, only the Gaussian wake model will be explained in depth. See Annoni et al. (2018), Martínez-Tossas et al. (2019), and Bay et al. (2019) for details on additional models in FLORIS.

Gaussian wake
The Gaussian wake model is comprised from a series of papers, including Bastankhah and Porté-Agel (2014), Bastankhah and Porté-Agel (2016), Abkar and Porté-Agel (2015), Niayifar and Porté-Agel (2015), and Dilip and Porté-Agel (2017). This model is a method of calculation that is integrated into the structure of all FLORIS wake modeling and control tools. It integrates the concepts of the Bastankhah and Porté-Agel wake deflection model, the self-similar velocity deficit model, and elements of atmospheric stability into one comprehensive method based off of the concept of a Gaussian wake (Pope, 2000). This section will describe the different concepts that are implemented in this model.

Self-similar velocity deficit
The Gaussian model computes the streamwise velocity deficit at any point in a turbine's wake by using analytical formulations of Reynolds-averaged Navier-Stokes (RANS) equations to an assumed Gaussian wake profile. The Gaussian wake is based on the self-similarity theory used for free shear flows (Pope, 2000) and is developed under the assumption of no pressure gradients within the initial undisturbed free-stream flow and uniform flat terrain (Bastankhah and Porté-Agel, 2014). To calculate the velocity deficit, u(x, y, z), behind the rotor of a turbine, where U ∞ is the freestream velocity; x, y, and z represent the spatial coordinates in the streamwise, spanwise, and vertical directions, respectively; and z h is the turbine hub height. C is the velocity deficit at the wake center; δ represents the wake deflection computed with equations from Bastankhah and Porté-Agel (2016); and σ denotes the wake width in the lateral (y) and vertical (z) directions. The subscript "0" references a term's initial value at the start of the far wake. The wake width in the y and z directions, σ y and σ z , is determined by the thrust coefficient, C T , and the wake expansion rate, which is parameterized by k y and k z : where D is the rotor diameter, u R is the velocity at the rotor, γ denotes the turbine's yaw offset, and u 0 represents the maximum velocity deficit in the wake. Parameters k y and k z are dependent on the value the ambient turbulence intensity, I 0 , as noted in Eq. (8).
The findings of Abkar and Porté-Agel (2015) demonstrate that k y and k z grow at different rates, but in order to simplify the model, k y and k z are usually set as equal. The total velocity deficit at any point in the domain of fluid flow can then be calculated by combining the wakes using the sum-of-squares method described in Katic et al. (1986).
In the scope of this study, it is important to note that the introduction of spatial heterogeneity in initial wind conditions (which is a key principle in the proposed model) violates the original assumption of no pressure gradient for the derivation of the Gaussian wake model. Although this limits the model's ability to conserve key principles that govern the physical dynamics of fluid flow, the results of this study show that the measured improvements in model accuracy outweigh the consequences of incomplete conservation. In Brogna et al. (2020), a modified Gaussian wake model is implemented to simulate wind farms in complex terrain, but the spatial U ∞ evolution is considered only in the superposition of wakes and is omitted for the calculation of the velocity itself. The benefits of an approach similar to this could be investigated in future FLORIS developments to improve overall momentum conservation for the heterogeneous model.

Atmospheric stability
The Gaussian model also implements methods proposed by Abkar and Porté-Agel (2015) and Niayifar and Porté-Agel (2015), which characterize the effects of atmospheric stability by analyzing the levels of veer, shear, and changes to turbulence intensity in the fluid flow. Stull (2012) discusses that an accurate representation of atmospheric stability requires the measurement of many other variables in the atmosphere; but without detailed recordings of elements such as temperature profiles and vertical flux, the three chosen parameters are able to give a rough idea of the state of the atmosphere in the FLORIS model.
To implement the effects of shear, α s , the power-log law of wind is used to define the initial wind speed in the flow field, U init : where a high shear coefficient (α s > 0.2) is indicative of stable atmospheric conditions, and a low shear coefficient (α s < 0.2) characterizes an unstable atmosphere (Stull, 2012). The Gaussian model was designed to avoid the inaccuracies caused by neglecting the effects of turbulence intensity by implementing methods introduced by Niayifar and Porté-Agel (2015). This also includes added turbulence caused by nearby turbine operation to more accurately calculate the rate of wake expansion. Many other linear-flow models use a constant parameter that defines the rate of wake expansion and has no dependency on the operating conditions of the turbine (Jensen, 1983). From the concepts of Niayifar and Porté-Agel (2015), the Gaussian model relates the rate of wake expansion in the lateral and vertical directions directly to the ambient turbulence intensity present at a turbine and two tuned parameters, k a = 0.38371 and k b = 0.003678: The turbulence intensity, I , is calculated by superimposing the initial ambient turbulence intensity (I 0 ) with the sum of the added turbulence caused by the operation of each influencing upstream turbine, j and I + j . The following relationship is used in FLORIS to calculate the ambient turbulence intensity at a given turbine with respect to neighboring turbine wakes: N refers to the number of upstream turbines that create a wake that adds to the ambient turbulence intensity at a downstream turbine's location. In Niayifar and Porté-Agel (2015), this number was assumed to be one, and the closest turbine was only taken into account because it would theoretically give the maximum amount of added turbulence. In the Gaussian model used in FLORIS, all turbines within a distance of 15D upstream and 2D in the spanwise (y) direction are included. Although the saturation effects of turbulence are not yet fully understood in this context, this formulation was shown to be a more accurate method of calculating added turbulence intensity in the findings Chamorro and Porté-Agel (2011), which found that turbulence intensity typically accumulates over two to three turbine rows but then levels off to an equilibrium at this point. Based on the original definition proposed in Crespo and Hernández (1996), the following expression in Eq. (10) has been tuned through comparisons to high-fidelity computational fluid dynamics (CFD) simulations (King et al., 2020b) and several field studies (Fleming et al., , 2020b to accurately calculate the added turbulence due to upstream turbine j : where D j denotes the diameter of turbine j , and A overlap refers to the fraction of the rotor-swept area of the downstream turbine that intersects with the cross-sectional area of the wake from the upstream turbine. The axial induction factor, a j , is evaluated based on the value of C T , as defined in Burton et al. (2002) and Bastankhah and Porté-Agel (2016).
As noted earlier, the Gaussian wake model was developed under the assumption of flat terrain. Since the heterogeneous model was specifically designed to best benefit wind farms located in complex terrain, it is important to know the consequences of violating this assumption. In Fleming et al. (2020b), a field study is presented that focuses on analyzing the performance of the tuned parameters in Eq. 10, by comparing two campaigns located in comparatively simple and complex terrains. The findings of this study indicate that inaccurate tuning of the tuned variables may worsen FLORIS's typical tendency to underpredict wake losses in areas with complex terrain.

Changes to the FLORIS model
Previously, FLORIS derived the initial wind speed, wind direction, and turbulence intensity by using one value to represent the entire flow-field domain. In this article, we describe the modifications to FLORIS to accommodate heterogeneous flows. This section will explain the methods used to calculate wakes based on the gradient of values observed in the undisturbed flow field without wake effects. The motivation behind this development was to create a more detailed characterization of the initial state of the atmosphere, which leads to improvements in the power predictions of a wind farm.

Initializing the heterogeneous flow field
To implement heterogeneity in FLORIS, an interpolation is performed based on several input values assigned to spatially varying coordinates inside or adjacent to the wind farm (see Fig. 2a). These initial inputs are used to approximate the value of atmospheric characteristics at the location of every turbine within the wind farm and at each individual grid point of the FLORIS flow field. FLORIS performs methods of interpolation and extrapolation using software packages provided by SciPy: an open-source scientific computing library for the Python programming language (Virtanen et al., 2020). The packages used in this method include a piecewise linear interpolant and a nearest-neighbor interpolant, which are combined to create an algorithm that calculates a unique value for each x and y coordinate within the flow field. Figure 1 shows a pseudo-code diagram of this process for reference.
The process begins with implementing a piecewise linear interpolation method for all points within the region defined by the input coordinates. First, Delaunay triangulation is performed using the Quickhull algorithm discussed in Barber et al. (1996). This method forms triangular connections between input points, based on their relative coordinates, and defines each triangle by ensuring its circumcircle remains empty. The result of this triangulation generates a mesh of triangular elements called a simplicial complex. Further details on the concept of Delaunay triangulation are explained in depth in Shewchuk (1999) and Barber et al. (1996).
The next step in determining the interpolated values is to use the established triangular elements to perform barycentric interpolation. During this step, the barycentric coordinates of each point of interest are determined relative to the triangular element in which it resides. Based on each set of barycentric coordinates, the interpolated result is calculated using a weighted average of the values defined at the triangle's vertices (Floater, 2015). A visual depiction of the methods utilized in this piecewise interpolation method (Delaunay triangulation and barycentric interpolation) is shown in Fig. 2b. After these processes are complete, FLORIS assigns the interpolated values to each flow-field grid point and turbine location inside the triangulated region bounded by the input coordinates. Any points that fall outside of this region must be determined through additional extrapolation processes.
Linear barycentric interpolation was chosen be implemented for this step because it is relatively efficient in computation and can be easily implemented without requiring any input parameters other than the locations and values of wind measurements. Although, it must be noted that the accuracy of the interpolated values is dependent on the quality of input measurements provided, the complexity of the terrain geometry, and the weather patterns observed in the physical wind farm.
The extrapolation process implements a nearest-neighbor interpolant to calculate all remaining unknown values. Using the recently interpolated point values in addition to the original input values, this method operates by selecting a single value at the nearest location to the point being extrapolated and assigning this nearest value to the extrapolated point. A visualization of this calculation is depicted in Fig. 3.
The nearest-neighbor extrapolation method was chosen because it defines a feasible relationship between input measurements and does not attempt to extrapolate using a formula derived from a curve-fitting or trend-predictive algorithm. Many other extrapolation methods attempt to predict  a rate of change outward of the interpolation domain by implementing a function that approximates a predicted progression of extrapolated values. For example, it was found that the analytic continuation of radial basis functions (RBFs) and fitted polynomial splines outside of the initial domain often produced a non-feasible output that did not respect the physical limitations of the atmospheric characteristic being extrapolated. Although it was speculated that these methods could likely be adjusted with tuning factors to fit extrapolated data within feasible bounds, efforts to do this were not explored in this study. Instead, the nearest-neighbor algorithm was chosen to simplify implementation of realistic extrapolation within the model.
When solving for the interpolated and extrapolated values for turbulence intensity and wind speed, values are easily computed because they are defined by values on a noncyclical scale. Because wind direction is represented using angles in degrees, the interpolation and extrapolation methods must be circular. The issue of interpolating circular data was addressed by simply computing the interpolation twice for each angle of wind direction, : once for the cosine component, α, and again for the sine component, β. The wind direction in a wind farm, , can be defined as where α = cos and β = sin . After is computed, the wind direction interpolation can then be defined for the entire wind farm.
It should be noted that the vertical (z) dimension is not considered when interpolating and extrapolating from the atmospheric inputs. Instead, all input values are assumed to be at the same z location, and the interpolation is performed on a two-dimensional plane at this height. Although this approximation may result in a less accurate result, this approach allows the interpolation and extrapolation algorithm to operate with less computational cost. Differences in wind speed due to variations in turbine hub height are calculated using the power law in Eq. (7), as described in Sect. 3.2.

Heterogeneous wind speed
Before FLORIS performs any calculations for velocity deficit in wakes, it first assigns an initial value of wind speed (U ) to each grid point in the flow-field grid. In a homogeneous case, these grid points would all have the same value across an x-y plane, but in a heterogeneous case, these grid points all have different values, dependent on the initial values that have already been established through interpolation. After U is defined at each grid point, the wind speed values at each x, y, and z coordinate in the flow-field domain are defined as U init , calculated using the power law in Eq. (7). From this point, the calculation of wakes proceeds in the same way as the homogeneous cases, with the exception of a more complex algorithm for accounting for changes in wind direction, as explained in Sect. 3.3. The velocity deficit behind each turbine is calculated by applying Eq. (4) from Sect. 2.3.1, where the free-stream velocity (U ∞ ) in Eq. (4) is defined as the local U init values at each flow-field grid point. Figure 4 shows visualizations of the resulting wakes after subtracting the calculated velocity deficit from the initial free-stream velocity at each flow-field grid point.

Heterogeneous wind direction
Similar to wind speed, an interpolation of wind direction is initially established across the flow-field grid through the methods of interpolation discussed in Sect. 3.1. The input values of wind direction are defined so that 270 • represents wind movement from west to east (see Fig. 5a), and then once FLORIS begins computations with these wind directions the values are converted so that 0 • represents the wind traveling from west to east (see Fig. 5b). Using these wind direction values, the turbine coordinates are rotated about the center of the simulation domain at these angles, as exemplified in Fig. 5b.
Using the rotated turbine map shown in Fig. 5b for reference, the flow field is adjusted to calculate each turbine wake independently, starting with the turbine that is the furthest upstream. To initiate the rotation of the flow-field grid, the grid points are rotated to the angle that is defined at the given turbine. This initial step is exemplified in Fig. 6 for the calculation of the velocity deficit behind turbine T6 only, but this will also be repeated for each turbine in the entire wind farm. This step is necessary to put the original non-rotated grid points in a frame of reference relative to the each specific turbine as their particular wake is calculated.
Next, to calculate the velocity deficit caused by each turbine's wake, all of the grid points in the flow field are rotated to replicate the effects of changing wind direction. These rotated grid points represent the redirection of the flow in response to changing wind direction within the flow field (see Fig. 7a). Once the velocity deficit has been calculated using the rotated grid points, the points are rotated back to their original positions in the flow field. Figure 7b shows the product of the final step, where the calculated velocity deficit is subtracted from the initial free-stream velocity at each flowfield grid point to reveal the resulting shape of the wake.
As discussed in Sect. 2.2, there is a minor computational expense in simulating the flow field independently for each turbine in the wind farm. This is because FLORIS determines a unique set of rotated grid points relative to the wind direction and coordinates of each turbine separately. The grid spacing in the streamwise (x) direction relative to the direction of flow is kept uniform throughout each iteration of the rotated grid, but the spanwise (y) spacing is adjusted with respect to the local wind direction inside the flow field. This allows the model to replicate a gradual change in wind direction throughout the flow field. The resulting flow-field wake calculation is shown in Fig. 8.
The grid point spacing in the x direction must be kept constant to avoid elongation or distortion of wake propagation and placement. Because the grid spacing in the y direction is not kept uniform, it must be noted that this capability of emulating a gradual change in wind direction may prevent the model from conserving momentum in some situations. Methods of enforcing uniform spacing in the y direction for each individual turbine wake have been developed but are not currently implemented because doing so limits the model's ability to create a gradient of wind directions within the flow field. In future work, methods of enforcing momentum conservation in this algorithm will be further investigated.
To further exemplify the applications of this functionality, Fig. 9 shows a more complex simulation of non-constant heterogeneous wind direction simulation in an irregularly spaced wind farm. The steps that FLORIS performs to evaluate this flow condition are identical to the ones displayed in Figs. 5-7, except it is personalized to the more complex variations of the depicted state of flow.
It is important to consider that this model was not designed to calculate the effects of changes in wind direction that are extremely dynamic. A change in wind direction that is too drastic will cause grid points in the rotated flow-field grid (the red points shown in Fig. 7a) to overlap each other within the same coordinate system, which may result in erroneous assignment of velocity deficit to these overlapped points. The limiting amount of wind direction change for the heterogeneous model is therefore the amount that causes the flowfield grid points to overlap in this manner. This limit must be determined for each wind farm independently, since it is dependent on the site-specific layout geometry and wind conditions of each case.
Although it may be possible for the wind direction within a wind farm to change this drastically, these conditions often involve multiple adjacent domains of flow that are separated by a boundary, which are difficult to represent in this model. These weather conditions are also most often observed in in-   stances of lower wind speeds and therefore can be considered not as lucrative in regards to power production. Plans for future developments to FLORIS involve designing a more inclusive model that is capable of mitigating issues concerning rapid changes in wind direction.

Heterogeneous turbulence intensity
The geographic distribution of turbulence intensity is established for the initial state of the flow field through the interpolation methods discussed in Sect. 3.1. This strategy of defining a more detailed variation of turbulence intensity in the flow field makes the approximation of wake dissipation and deflection more accurate, therefore improving the estimation of the effect of nearby turbine operation within a wind farm. The implementation of heterogeneous turbulence intensity and heterogeneous wind speed are similar, in that the initial heterogeneous conditions are established throughout the flow field by interpolating from the input values, and then waked conditions are updated throughout FLORIS computations of flow-field interactions. During the calculation of wakes, the ambient turbulence intensity that is initially defined at each turbine location is continuously recalculated to account for     added turbulence intensity resulting from turbine wakes up to 15D upstream, as previously discussed in Sect. 2.3.2 and in Niayifar and Porté-Agel (2015). A horizontal plane of a FLORIS simulation featuring heterogeneous turbulence intensity can be observed in Fig. 10.
It is important to note that in the interest of conserving computational efficiency, calculations for evaluating the rate of wake expansion and recovery are only dependent on the updated turbulence intensity at the location of the turbine creating the wake.

Turbulence correction
In addition to the heterogeneous features, developments were also made to reduce inaccuracies in power output predictions caused by turbulent operating conditions. As mentioned in Sect. 2.1, the accuracy of the zero-turbulence power curve is compromised in conditions of varying turbulence intensity. The revised power calculation, presented in this section, includes a parameter that approximates the effect of turbulence intensity on the power output of a turbine in a wind farm.
Specifically, this approach adjusts the power output with respect to the level of turbulence intensity at a turbine. The adjusted power is calculated by using distribution of the wind speed fluctuations at the turbine, based on calculations that consider the original wind speed and the standard deviation in wind speed. The first step in this algorithm is to create a normalized probability density function, f (x), of wind speeds, x, evenly distributed within the domain of 1 standard deviation from the mean wind speed, µ. The standard deviation, σ , is determined by multiplying the turbulence intensity at the turbine by the mean wind speed, µ. Wind speeds that are greater than the cutout wind speed are omitted.
The value of the power coefficient, C P , in the power table is also determined at each wind speed, x i , and at the original wind speed (µ). The ratio of the adjusted power (P adj ) to the original value of power (P 0 ) is referred to as the turbulence parameter, . The turbulence parameter can be calculated by summing the weighted adjusted values of power in the following expression, for each wind speed, x i , in the domain of the probability density function, f (x i ): where the integral of f (x i ) is approximated by taking 100 samples of the f (x i ). The resulting power curves depending on turbulence intensity are shown in Fig. 11. As the turbulence intensity increases, the power output increases in region 2 and decreases across region 3.
The following expression may be used to calculate the final value of adjusted power output, P adj , with respect to the current turbulence intensity at a turbine: where γ is the yaw angle of the turbine, and represents the turbulence parameter. The value of must always be greater than zero. In future work, this turbulence-correction model could be improved by implementing a similar consideration of the thrust coefficient, C T . Because the velocity deficit computations in this model rely on the value of C T , it may be advantageous to expand this method to calculate an adjustment parameter for the effects of turbulence on rotor thrust.
It is important to note that similar models have been developed that incorporate methods of turbulence renormalization based on machine-learned or empirically derived data (Clifton and Wagner, 2014). The proposed method discussed in this section was developed to attempt to represent the variation of power output due to turbulence effects, while using a simple strategy that is not dependent on the availability of data other than the current wind farm atmospheric measurements and the power curve provided by the turbine manufacturer. In future work, it may be advantageous to incorporate more complex techniques that are able to capture the effects of turbulence intensity with greater detail and accuracy. Figure 11. Adjusted power curve for the NREL 5 MW reference turbine for different turbulence intensities. The dashed lines denote the cut-in, rated, and cutout wind speeds and also represent the boundaries of the first, second, and third regions, respectively.

Model validation and analysis
Two validation studies are presented to analyze the effectiveness of the adapted FLORIS wake model. In Sect. 4.1, a 38-turbine wind farm is simulated using the heterogeneous FLORIS model and compared to results from large eddy simulations to evaluate the accuracy of wind farm power predictions. Additionally, Sect. 4.2 presents a study where a large wind farm is simulated using the heterogeneous FLORIS model and turbulence power calculations. The results of these FLORIS simulations are compared to the wind farm's SCADA data records to further evaluate the model's power prediction accuracy and wake model performance.

Comparisons to LES
This section presents a validation study that evaluates the power prediction accuracy of the proposed heterogeneous FLORIS model in comparison to large eddy simulations of a 38-turbine array, calculated using NREL's tool Simulator fOr Wind Farm Applications (SOWFA) (Fleming et al., 2013). The simulated wind farm contains 38 turbines modeled with NREL's 5 MW reference turbine design criteria (Jonkman et al., 2009) and arranged in a concentric circular layout similar to Thomas et al. (2019), Fleming et al. (2020a), and King et al. (2020a).
Twelve test cases were evaluated for this study; each was simulated using different wind directions, varying from 10 to 340 • in 30 • intervals. Spatially homogeneous inputs were used to simulate wind direction and turbulence intensity, where turbulence intensity was at 9 % for all cases. The freestream wind speed remained close to 8 m s −1 for all cases, with minor spatial variations. FLORIS wind speed inputs were obtained by extracting the free-stream velocity from LES results at locations upwind of turbines in undisturbed flow. These extracted input values create a velocity gradient in the direction normal to the wind direction in the heterogeneous FLORIS model, as seen on the right in Fig. 12. The wind speed input for the homogeneous FLORIS model was obtained by taking an average of the heterogeneous input values for each case.
To analyze the effects of wake losses for the simulated wind farm, additional heterogeneous and homogeneous FLORIS simulations were conducted excluding all FLORIS wake calculations. Figure 13 shows the total wind farm output predictions from all four FLORIS simulation models and compares them to the LES case result (shown in black). The trends observed Fig. 13 indicate that the models which neglected wake loss calculations drastically overestimated total wind farm power output, with the heterogeneous model reporting more accurate power predictions overall. Alternatively, the FLORIS models that did perform wake loss calculations produced a much more accurate estimation of wind farm power output.
The absolute error of the total wind farm power output was calculated for each model to analyze power prediction accuracy at every wind direction evaluated in the case study. The results of these calculations are shown in Fig. 14, and each model's average absolute error over all wind directions is recorded in Table B1 for reference. In Table B1, an 18.9 % decrease in average absolute error is reported when using the heterogeneous wake model compared to the homogeneous wake model.
The mean absolute error (MAE) of power predictions at individual turbines in the wind farm was also calculated using Eq. (14).
where n is the number of turbines in the wind farm, P actual,i is the LES-generated wind farm power output for turbine i, and P model,i represents the predicted power output from the FLORIS model for turbine i. MAE was calculated for each FLORIS model at every wind direction in the case study. The resulting MAE values are shown in Fig. 15, and the average   MAE for all wind directions is reported in Table B2 for each FLORIS model. In comparison to Fig. 14, Fig. 15 indicates a more observable and consistent disparity in power output accuracy between the results of heterogeneous and homogeneous FLORIS wake simulations. Table B2 validates this observation by reporting a 19.5 % decrease in MAE at each turbine when using the heterogeneous model. The data provided in this comparison confirm that the proposed heterogeneous model offers substantial advancements in the generation of accurate power predictions at individual turbines within a wind farm.

Comparisons to wind farm SCADA data
This section summarizes a validation study presenting comparisons of FLORIS power predictions to SCADA-recorded power outputs from an observed operational wind farm. A large, utility-scale wind farm located within mountainous terrain was chosen for this study because it is often subject to unpredictable and dramatic shifts in weather conditions. More information regarding the physical layout and characteristics of this wind farm can be found in Appendix A. The motivation behind performing these simulations was to quantify the effect of the recent developments to FLORIS in reducing the error in power output predictions for wind farms in complex terrain. FLORIS simulations were performed using heterogeneous inputs of wind direction, turbulence intensity, and wind speed, which were taken from the wind farm's SCADA records. These inputs include four wind measurement values for each atmospheric characteristic, derived from meteorological (MET) tower measurements placed in various locations throughout the wind farm. Similar simulations were performed using an identical FLORIS model but with a singular homogeneous input for wind speed, wind direction, and turbulence intensity. These homogeneous inputs were derived by evaluating the average of the five heterogeneous input values at each time step. The resulting power output of all simulations was recorded with the inclusion of the turbulence correction and without. All cases were simulated using data averaged at time steps of 10 min over a range of 2 months.
In the following discussion, the results from all FLORIS simulations are presented and analyzed to determine the accuracy of power predictions from each test case. Figure 16 includes two horizontal planes showing a partial section of heterogeneous flow calculations during these simulations. This figure demonstrates the visual capabilities of the heterogeneous model and how the effects of the new wake calculations can be translated into visual information for further analysis of wake interactions within a wind farm.
Although these visualizations do not give direct estimates of power prediction, they are helpful in translating the input measurements into a form that characterizes the general behavior of wind farm dynamics for the interpretation of the observer. The cut plane visualization is helpful in performing qualitative analysis of turbine wake interactions and is more useful when displaying the estimated weather conditions characteristic of each location in the flow field, which is improved in the heterogeneous model.
When comparing the performance of the simulations, the calculated power output was tabulated and compared for accuracy. In Fig. 17, the sum of wind farm power output from each FLORIS simulation is normalized with respect to the rated power output for the wind farm and plotted along with the recorded SCADA output. This approach highlights any Figure 17. Power output calculated by FLORIS for homogeneous (red), heterogeneous without the turbulence correction (blue), and heterogeneous with the turbulence correction (green), compared with SCADA data shown in black. Each shaded region represents the difference between predictions of power output and the measured power output from SCADA data. weaknesses in each model, relative to the overall performance of the others. A 24 h period was chosen to demonstrate how the models performed under average diurnal conditions. Figure 17 shows a day with relatively variant weather conditions and many rapid shifts in power output.
In Fig. 17, it is evident that the heterogeneous models are predicting the power output more accurately than the homogeneous model. The trend line of the heterogeneous simulations consistently follows closer to the line representing the power output recorded from SCADA data. Additionally, the heterogeneous simulation that included turbulence-intensity corrections showed an extra advantage in estimating turbine performance, following closely to the trend line of the heterogeneous simulation, and also reliably contributing errorreducing improvements to the heterogeneous model. While this juxtaposition is effective in ranking each model's ability to estimate total farm power output, it should be noted this comparison only indicates the accuracy of a calculations for the entire wind farm power output collectively, without considering the accuracy at each turbine individually.
It is possible for wake models to overpredict the power output of some turbines, and underpredict others, in a way that produces a total wind farm power estimate that seems accurate but is not using reliable and precise methods of calculation. To verify that the recent additions to FLORIS have improved the power-predicting capabilities, it must be confirmed that the new model produces a consistently accurate estimate with respect to each iteration in the time series and each turbine within the wind farm individually. To prove this model's consistency in accuracy, the normalized absolute error was calculated at each turbine at each iteration of the time series for this same day. The sum of the absolute error at all turbines within the wind farm is calculated for each simulation model at each time iteration. To calculate the sum of absolute error (SAE) for all turbines, the following formula was applied to each time iteration of the simulation.
where n is the number of turbines in the wind farm, P actual,i is the measured power output of turbine i, and P model,i is the predicted power output of turbine i from a given FLORIS model. The results of each FLORIS model were calculated and plotted on the same set of axes in Fig. 18. The trends observed in Fig. 18 exhibit similar characteristics that indicate the accuracy of the model at each turbine is increasing with the application of the heterogeneous model and turbulence-intensity correction parameter. The heterogeneous model reliably produces less error when calculating the power at each turbine over the time series, which ensures that the power predictions of the entire farm are not selfcompensating because of simultaneous overpredictions and underpredictions of individual turbine outputs. Furthermore, if Fig. 18 is analyzed with respect to the trends of normalized power in Fig. 17, it is evident that the addition of heterogeneity and turbulence-intensity corrections contributes to improving the accuracy of FLORIS power predictions in instances of overprediction and underprediction, as well as transitions between the two with relative consistency.
To ensure these same trends of accuracy persist over the entire 2-month period, the percent error of the total wind farm power output was calculated at each time-step iteration using the following equation.
where P actual is the measured power output of the wind farm, and P model is the power output of the wind farm predicted by a given FLORIS model. The results of these calculations were grouped into three separate domains: wind speeds of less than 5 m s −1 , wind speeds in the range of 5 to 11 m s −1 , and wind speeds greater than 11 m s −1 . Time iterations when wind speed was less than the cut-in wind speed (2.5 m s −1 ) were considered negligent in regards to power production and therefore omitted from the data set. A histogram of the percent error in each wind speed domain was computed over the entire time series to display the distribution of error with respect to each simulation (Fig. 19).
Although the plots for the wind speed domains vary slightly in distribution, it is clear that each histogram exemplifies a trend toward accuracy in simulations that incorporate heterogeneity and turbulence-correction calculations. It is important to note that only the data points shown in the percent-error range of each histogram were used to calculate the respective binned averages. The outliers were omitted because they tend to skew the presentation of the data set in a way that obscures the actual trend of data.
The mean absolute percent error (MAPE) values of all time-step iterations are also reported in Table 1. The data for this table were calculated by evaluating the percent error of FLORIS power predictions for the full wind farm at each time step and then solving for the mean over the entire time series. This calculation is expressed as where n is the number of time steps in the total simulation, P actual,i is the recorded power output of the wind farm at time step i, and P model,i denotes the predicted power output from the FLORIS model at time step i. When comparing the MAPE values in Table 1 with the histograms of Fig. 19, an increase in MAPE is observed in Table 1 for lower wind speeds of simulations that implemented heterogeneous and turbulence correction models. This is a trend that is not characteristic of the histograms depicted in Fig. 19b. In reference to this observation, it is important to note that the metric of MAPE penalizes overpredictions with more weight than underpredictions. Furthermore, MAPE calculates mean with equal weight for all time steps in the data set, which is not always ideal for an indication of overall farm power output accuracy. It is possible that the reported increase in MAPE with lower wind speeds may be an indication that the heterogeneous and turbulence-intensity correction models tend to produce more frequent overpredictions Figure 19. Percent error of all three FLORIS models, plotted for comparison within varying ranges of wind speeds. of power output in conditions where wind speeds are near the cut-in speed. If this is true, further investigations may be conducted in future work to determine why this is happening and how it could be circumvented. Although MAPE is an informative metric for analyzing the average percent error relative to a specific power output range, methods that use unweighted averaging are sometimes misleading in the analysis of overall power prediction accuracy. The relative error during time-step iterations with lower power output can seem large, even when the absolute error is insignificant in comparison to the magnitude of total farm output.
A more comprehensive representation of relative model accuracy is presented in the following table, where the mean absolute error (MAE) is evaluated for total wind farm output. This was calculated by evaluating the absolute error at each time step and then taking the mean of these error values. This calculation is performed using Eq. (14), where n is the total number of time steps in the simulation, P actual,i is the recorded power output of the wind farm at time step i, and P model,i is the predicted power output from the FLORIS model at time step i.
By taking an average of absolute errors instead of relative errors, MAE is a more effective metric in representing the overall accuracy of total wind farm power prediction. The resulting MAE values are shown in Table 2, where a clear trend of increased accuracy is observed for models that implement heterogeneity and turbulence-adjustment calculations, including the cases where wind speeds are below 5 m s −1 .
Lastly, values of MAE were also calculated to represent the accuracy of the model at each individual turbine within the wind farm. Using this metric, Table 3 shows that simulations implementing the heterogeneous model and turbulence correction calculations outperformed the homogeneous model in the prediction of individual turbine power output. This should be expected, since the overall farm output in Ta-  ble 2 followed a similar trend. The marked improvement of power predictions at individual turbines suggests that the addition of the proposed heterogeneous and turbulence correction methods enhances the FLORIS wake model by simulating farm-flow interactions with more thorough detail and greater accuracy.
To analyze the influence of wake effects in this study, an identical set of simulations were performed excluding FLORIS wake loss calculations, and the results for MAE at the overall farm and individual turbine levels are reported in Tables B3 and B4 in Appendix B. These simulations seem to indicate that the wake losses at the observed subject wind farm are relatively small, due to the large spacing between turbines in the stream-wise direction. The study presented in Sect. 4.1 analyzes the performance of the proposed heterogeneous model in a wind farm with more influential wake losses to give a more detailed analysis of the effects of wake losses.
As noted in Sect. 3.3, the implementation of methods utilized to simulate spatially variant wind direction causes the heterogeneous model to be marginally less efficient in computation. To quantify this increased computational cost, each simulation was timed in this study. On average, these time records show that the simulations using the heterogeneous model took less than 10 % longer to compute than those using the homogeneous model. The choice to sacrifice computational efficiency in the heterogeneous model was seen as a necessary trade-off to achieve greater detail and accuracy in simulations of more dynamic environments. Future developments to FLORIS will attempt to optimize the efficiency of this model and reduce the time necessary to simulate the effects of changing wind direction.

Conclusions
This article introduces a method to include heterogeneous flow fields into the FLORIS simulation tool, as well as a turbulence correction to the power reported at each turbine. To analyze the developed model's improvements in accuracy, several FLORIS simulations with and without these changes were compared to large eddy simulations and SCADA data from a utility-scale wind farm. The results of the FLORIS simulations indicate that these two modifications improve power predictions of the wind farm at the turbine and wind farm level. The increased accuracy of this model's power prediction capabilities shows that this method is more precise in predicting farm-flow interaction in heterogeneous and turbulent environments, which previous versions of FLORIS were not able to simulate.
Overall, the heterogeneous and turbulence-intensity correction modifications presented in this article showed a positive effect on the accuracy of FLORIS capabilities. This improved model provides a more detailed quantitative and qualitative analysis of wind farm flow, including the demonstration of heterogeneous flow in cut-plane velocity plots and improved accuracy in power prediction at individual turbines as well as total wind farm power output. Comparing FLORIS power predictions to LES, the heterogeneous FLORIS model showed an 18.9 % decrease in mean absolute error (MAE) for total wind farm power output and a 19.5 % decrease in MAE for individual turbine power predictions compared to the homogeneous FLORIS model. In comparisons to SCADA data, FLORIS simulations that implemented the heterogeneous flow model showed a 14.6 % decrease in MAE for wind farm power output predictions compared to homogeneous model simulations. With the use of the proposed turbulence-intensity correction method in ad-dition to the heterogeneous model, the MAE in farm power output predictions showed a 31.42 % MAE decrease compared to the homogeneous model.
These modifications to FLORIS have outlined a framework for a wake model that features atmospheric heterogeneity and turbulence-intensity corrections to the power curve and provides a platform for further developments in this area of research. In agreement with this study, the findings of Fleming et al. (2020a) also indicate that this model shows promise in enhancing the performance of FLORIS's existing wind farm optimization controls, in addition to improving the accuracy of wind farm power predictions.
Further studies relating to the effectiveness of this model when applied to wind farm controls would be very beneficial in determining future developments to these algorithms. Additionally, more extensive investigations should be considered to evaluate the efficacy of the proposed model in a wider variety of operational conditions, particularly those with lower wind speeds and extreme variations in wind direction. Other future work will investigate alternative interpolation methods for the flow field that consider the wind farm terrain map, capabilities for simulating more dynamic changes in wind direction, implementing enforcement of momentum conservation, and optimizing the model's computational efficiency.  Table A1. This table lists several key attributes that characterize the nature of the terrain and turbine layout of the observed wind farm discussed in Sect. 4.2. Distance values are reported relative to the average turbine rotor diameter (D). Spanwise and stream-wise directions are defined to be perpendicular and parallel to the average wind direction during the wind farm, respectively.

Measured quantity
Distance in terms of average rotor diameter (D) Average stream-wise inter-distance 20.0D Average spanwise inter-distance 2.0D Range of elevation variation 2.2D Review statement. This paper was edited by Sandrine Aubrun and reviewed by two anonymous referees.