*2.2. WindSim*

WindSim is a commercial CFD package that simulates flow over wind farms in complex terrain. The program solves the steady Reynolds-Averaged Navier–Stokes (RANS) equations using a two-equation turbulence closure model. Steady RANS simulates time-averaged flow fields assuming a statistically-stationary condition. This study focuses on the evaluation of this approach and its associated models by comparing simulation results with time-averaged field measurements. When a time average is taken, transient phenomena are smoothed out and become invisible. Hence, unsteady effects such as time-varying large-scale atmospheric forcing, topography-induced vortex shedding, and turbine wake meandering cannot be captured by steady RANS. More advanced methods (e.g., LES models) capable of capturing the unsteady effects are still too computationally expensive for commercial use in wind energy applications.

The Standard *k* − *ε* model (STD), the modified *k* − *ε* model of the Re-Normalization Group (RNG) version [12], and the *k* − *ω* turbulence model of Wilcox [13] were considered in this study. For flow over complex terrain, some studies showed that the RNG *k* − *ε* model produced promising results, a finding corroborated by Peralta et al. [14], while some other studies showed the superiority of the Wilcox *k* − *ω* model over the STD and RNG models [15]. Capturing the effects of forest is essential for this case study. In WindSim, forest can be modeled by the indirect approach mentioned before or by including porous cells with momentum sinks and turbulence sources [16] in the computational grid for areas that include forest. The latter is called the forest model. Our testing results (not shown) indicated that the forest model yields more realistic results than the indirect approach does. Turbine wake effects can be simulated directly by the use of an Actuator Disc (AD) model [17]. However, in WindSim, the AD model cannot be activated together with the forest model. Hence, in this study, turbine wake effects were modeled through the analytical approach. WindSim has implemented the analytical wake models from Jensen, Larsen, and Ishihara [18]. The accuracy of the three models to predict the observed power production was evaluated.

WindSim can optionally account for atmospheric stability by additionally solving the temperature equation. However, this feature requires several inputs that were not available from the measured data. Instead, we validated the assumption of a neutral boundary layer by examining the measurements. Within the selected range of wind direction and speed, we further filtered by time of day, keeping only wind events from dusk and dawn, when the atmosphere was assumed to be neutral. Comparison with the full dataset showed no significant change in time-averaged wind behavior at any of the turbines. We therefore concluded that the assumption of a neutral boundary layer for the WindSim simulations was valid for the model evaluation study here.

#### *2.3. Boundary Conditions and Numerical Settings*

The computational domain, surface elevation data, and turbine locations are shown in Figure 4. For each simulation case, the domain was rotated to make the *x*-axis along the prevailing wind direction, so there was only one inlet (at *x* = 0) and one outlet (at *x* = *Lx*). At the inlet, boundary conditions are given as fully-developed flow profiles taking into account the given roughness at the border and the boundary-layer height *LB* [19]. For the wind speed, the well-known logarithmic profile is defined from the ground up to *LB*, and above this height, the profile is constant. Here, *LB* was set to 1000 m above the mean surface elevation, and the constant speed above *LB* was set to 15 m/s so that the simulated wind speed at Turbine 2 was around 8.5 m/s, which is the median of the wind speed range applied to filter the data. At the outlet, zero gradient boundary conditions are imposed, meaning that a zero diffusion flux for all flow variables is assumed. On the lateral sides, symmetric conditions are applied. The upper boundary condition is specified as fixed pressure. The bottom boundary condition is no penetration together with the equilibrium log-law wall functions.

WindSim uses a Cartesian grid in the horizontal plane and terrain-following grid points in the vertical direction with tighter spacing closer to the ground level. The number of vertical grid points was set to the maximum (60). Test simulations with four different numerical settings as detailed in Table 1 were performed. Since the purpose of those simulations was to check the convergence of numerical results with regard to grid resolution and domain size, wake effects were not considered, and forest was modeled by the less expansive indirect approach.

For the evaluation simulations, the forest model was used. The height of the forest was set to 20 m, which is the mean height of the trees in the region according to a survey [20]. The number of grid cells in the vertical direction for modeling the forest was set to five, corresponding to *dmin z* = 4 m. According to the table in WindSim, the forest resistive force constant *C*2 was set to 0.01, twice that of the default value, because the forest at Mont Crosin was sparse, but dominated by *Picea abies* and *Abies*

*alba*, which are evergreen coniferous trees with a higher leaf area index. The influence of *C*2 on the results is presented in the next section.


**Table 1.** Information of the numerical settings.

**Figure 4.** The computational domain for the WindSim simulations. Here, *z* is the elevation from sea level in m and *zs* is the surface elevation.

#### **3. Results and Discussion**

Figure 5 shows the normalized wind speeds at the hub height of turbines predicted by WindSim using the standard *k* − *ε* model with the four different numerical settings defined in Table 1. It can be seen that the differences between the results of these numerical tests (except for S4, which had the coarsest grid resolution) were rather small. This indicates that the results presented in the following with the numerical setting S1 did not depend on the domain size and grid resolution. S1 (medium grid resolution) was ultimately chosen because it produced results much faster than using S3 (fine grid resolution). Furthermore, the nesting technique (using the results from a larger outer model with coarser resolution as boundary conditions for flow simulation over a smaller domain with higher resolution) was also tested for S1. The nested simulations did not further change the results (not shown). Therefore, the influence of inaccuracies in the assumed boundary conditions on the results of S1 can be regarded as negligible.

Normalized turbine power outputs predicted by WindSim using three different turbulence models are compared with the wind farm SCADA data in Figure 6. Here, the analytical wake model of Ishihara was used, and the effect of multiple wakes was modeled by the linear superposition of the wake deficits. The predicted power outputs were obtained for three incoming wind directions (237◦, 240◦, and 243◦) and different wind speeds (around 8.5 m/s) at the reference turbine, then averaged to yield the mean values and standard deviations (error bars). It is shown that the results of Wilcox were all within the error bars of the data, while the results of STD and RNG largely under-predicted the power outputs of Turbines 1, 11, and 12. To have a quantitative measure of the model performance, the Root Mean Squared Error (RMSE) and Mean Bias (MB) were calculated as follows:

$$\text{RMSE} = \sqrt{\frac{1}{N} \sum\_{n=1}^{N} (P\_s^n - P\_o^n)^2} \tag{1}$$

$$\text{MB} = \frac{1}{N} \sum\_{n=1}^{N} \left( P\_s^{\text{II}} - P\_o^{\text{II}} \right) \tag{2}$$

where *Pns* is the simulated mean power output of turbine *n*, *Pno* is the observed mean power output of turbine *n*, and *N* is the number of turbines used for comparison. RMSE was 0.09 for the Wilcox model, 0.20 for the STD model, and 0.25 for the RNG model. MB was almost zero for the Wilcox model, −0.15 for the STD model, and −0.20 for the RNG model. Overall, it can be concluded that the Wilcox model outperformed the other two models in terms of predicting turbine power outputs in complex terrain.

To have a closer look at the different behaviors of the turbulence models, we plot the fields of predicted wind speed at the height of 95 m for the predominant incoming wind direction in Figure 7. It turns out that at the leeward side of the first hill (marked by the black triangle), where Turbine 2 is located, wind speeds predicted by the Wilcox model were lower than those predicted by the STD model and the RNG model. Since the power of Turbine 2 was used to normalize the results, this explains why the STD and RNG models tended to underestimate the normalized powers at the other turbines. This finding is consistent with other studies showing that the Wilcox model is able to predict mean velocity and turbulent kinetic energy that are closer to the measurements than the other models [15]. The Wilcox model involves the solution of transport equations for the turbulent kinetic energy *k* and the specific dissipation rate *ω* = *ε*/*k* where *ε* is the dissipation rate of *k* [13,21]. Compared to the *k* − *ε* models, the *k* − *ω* model has several advantages, namely that: (1) the model is reported to perform better in mildly-separated flows; (2) the model is numerically very stable; (3) the low-Reynolds-number version is more economical and elegant in that it does not require the calculation of wall distances, additional source terms, and/or damping functions based on the friction velocity. It can be inferred from the results that, among those advantages, the first one is mainly responsible for the best performance of the Wilcox model found here.

**Figure 5.** Normalized wind speeds at the hub height of turbines predicted with the four different numerical settings defined in Table 1.

**Figure 6.** Normalized turbine power outputs observed by field measurement and predicted by WindSim with three different turbulence models: 1. Wilcox; 2. STD; 3. RNG.

**Figure 7.** Wind speeds at the height of 95 m predicted by the three turbulence models for the predominant incoming wind direction (from top to bottom: STD, RNG, Wilcox). The first hill is marked by the black triangle.

With the Wilcox turbulence model, the performances of the three analytical wake models implemented in WindSim were evaluated. Again, the linear superposition of velocity deficits was adopted to handle the multiple wakes. As shown in Figure 8, among the three wake models, the Ishihara model yielded the best overall results with an RMSE being 0.09 and an almost zero mean bias, while the Jensen model had an RMSE of 0.15 and an MB of −0.08, and the Larsen model had an RMSE of 0.21 and an MB of 0.10. The better performance of the Ishihara model may be due to the fact that it introduces a turbulence-dependent rate of wake expansion and adopts the Gaussian shape for the velocity deficit. Nevertheless, it is important to note that none of these analytical wake models considers the change of wake growth with topography due to the pressure gradient, which could be significant according to a recent study [22].

**Figure 8.** Normalized turbine power outputs observed by field measurement and predicted by WindSim with different analytical wake models: 1. Ishihara; 2. Jensen; 3. Larsen.

Figure 9 compares the results obtained by using two different approaches to calculate the superposition of multiple turbine wakes. It turns out that the linear superposition approach led to stronger multiple wake deficits for the last two downstream turbines (13 and 14) and predicted normalized powers that were in better agreemen<sup>t</sup> with the measurements, compared with the other approach that uses the square root of the sum of the squares of the velocity deficits. It is worth mentioning that similar behaviors of the two approaches were found in a study of the Horns Rev offshore wind farm [23].

Table 2 summarizes the prediction errors of the various combinations of modeling options. It is evident that the *k* − *ω* turbulence model of Wilcox together with the analytical wake model of Ishihara and the linear superposition of multiple wake deficits yielded the best performance. Some other combinations of turbulence and wake models were also tested (results not shown), and none of them outperformed the one recommended above. Nevertheless, for this case study, the forest modeling played a key role, and the results were sensitive to the choice of the forest resistive force constant *C*2, as shown in Figure 10.

**Figure 9.** Normalized turbine power outputs observed by field measurement and predicted by WindSim with the multiple wake effects modeled by (1) the linear superposition of velocity deficits and (2) the square root of the sum of the squares of velocity deficits.

**Table 2.** Data of the RMSE and Mean Bias (MB) of the predicted normalized powers for the model combinations.


**Figure 10.** Normalized turbine power outputs observed by field measurement and predicted by WindSim using the forest model with different *C*2 values: 1. 0.005; 2. 0.01; 3. 0.02.
