Large-eddy simulations (LES) are performed on the flow over a wind farm sited behind an abrupt rough-to-smooth surface roughness jump. The change in surface roughness affects both the first-order and second-order turbulent statistics. The usual deficit, i.e., the difference between the velocities upstream of the entire wind farm and downstream of a turbine, attains negative values close to the ground, which makes it difficult for modeling within the usual Gaussian radial-shape framework. A different definition, i.e., the difference in velocity at the same location with and without a turbine on a heterogeneous surface, is always positive and is amenable to Gaussian shape-based modeling. For the setup considered here, wind farms sited downstream of a surface roughness jump produce more power than a wind farm sited on a homogeneously rough surface. This increase is primarily because of the larger power generated by the downstream turbines and only slightly due to the increased power of the first-row turbine. The farm performance is affected by the distance between the abrupt change in surface roughness and the position of the first row of turbines. The wind farm performance is also dependent on the aerodynamic roughness upstream of the surface roughness jump. Two single-turbine analytical models and three wake-merging strategies are evaluated for their ability to predict the velocity deficits. A corrected form of the standard Gaussian model with a recently proposed wake-merging methodology, applicable for a varying background field, is found to be insensitive to the tunable model parameter and is consistently in line with the LES results.

## I. INTRODUCTION

The flow within the atmospheric boundary layer (ABL) directly impacts the energy extracted by a wind farm. Modern wind farms often span vast horizontal distances, over which the land is likely to undergo a change in its type and use. This leads to spatial variability of the momentum fluxes imposed by the ground surface on the ABL and wind farm flow. The momentum flux imposed by the ground is parameterized by the aerodynamic roughness height (*z*_{0}). Therefore, it is crucial to investigate the influence of surface roughness heterogeneities on ABL flow and the consequences for wind farms located on such heterogeneously rough terrain.^{6,27}

The impact of surface roughness heterogeneities on ABL flow has been investigated through several experiments or field observations^{7,9} and numerical simulations,^{1,27,33,37} with the most commonly studied configuration being that of an abrupt transition between a rougher and smoother surface in the normal direction to that of the mean wind. These investigations have demonstrated the formation of an internal boundary layer (IBL) at the roughness transition point, where the change in surface roughness modulates the mean streamwise velocity and turbulent fields. Several analytical models have been formulated for the mean streamwise wind speed^{1,9,13,14,28,30} and for the turbulence intensity (TI)^{27} over surfaces with heterogeneous roughnesses.

Several experimental and numerical studies have been reported of wind turbines/farms and are summarized in the review articles by Stevens and Meneveau^{38} and Porté-Agel *et al.*^{32} A number of simplified analytical models have been developed to predict the evolution of the wake of an isolated wind turbine. Among these, the model by Bastankhah and Porté-Agel^{2} particularly has been widely adopted because it conserves mass as well as momentum and assumes a Gaussian shape in the radial direction starting at moderate distances downstream of the turbine, which approximates the actual shape very closely. Engineering wind farm wake models (“bottom–up” models following the terminology of Stevens *et al.*^{40}) rely on strategies to merge the wakes of multiple turbines together and come up with a composite wake. The oldest examples of these wake models that account for interactions between multiple turbines are the linear merging^{26,29} and quadratic merging.^{23,47} A new wake merging strategy was proposed recently by Lanzilao and Meyers (LM),^{25} which, unlike previous models, allows for the background velocity field on which the wind turbine wake-induced deficits are merged, to be non-uniform. Almost all of the above-mentioned numerical simulations or analytical modeling studies have considered wind turbines/farms sited on homogeneously rough surfaces.

Studies involving surface roughness heterogeneities and wind turbines have been limited in number. Kethavath *et al.*^{24} studied the evolution of an isolated wind turbine wake with the turbine sited downstream of an abrupt transition in *z*_{0}. The surface roughness heterogeneity was shown to significantly alter the turbulent statistics. Furthermore, a few corrections to the Gaussian model were proposed, which account for a non-uniform incoming flow and incorporate the effect of a few terms in the momentum deficit transport equation that are usually neglected by other single-turbine models. The model by Kethavath *et al.*^{24} was shown to provide better predictions of the large-eddy simulation (LES) results of a single wind turbine wake sited behind an abrupt surface roughness jump than the Gaussian model.

Several studies have investigated coastal effects on offshore wind farms as reviewed in Schulz-Stellenfleth *et al.*^{36} For example, the Anholt wind farm in Denmark has been studied using a combination of mesoscale simulations, supervisory control and data acquisition (SCADA) data, and Reynolds-averaged Navier–Stokes (RANS) simulations for understanding the impact of horizontal wind-speed gradients, the effect of the distance between the coast and turbines, and assessing the accuracy of wake models in predicting power losses and annual energy production.^{31,45} The work of Van Der Laan *et al.*^{45} showed that the distance of the leading turbine in a wind farm array to the shoreline affects the wind speed seen and the power generated by the turbine. Recently, Goit and Önder^{18} investigated the effect of the land–sea transition on the performance of an offshore wind farm. The first row of turbines was placed 25 turbine diameters downstream of the land–sea transition. Incorporating the effect of this surface roughness transition into the simulations led to an enhanced wake recovery, greater inflow turbulence levels, and increased the kinetic energy entrainment. For the specific wind farm layout and surface roughnesses considered, the wind farm placed near the land–sea transition was shown to produce 17% higher power than a wind farm placed further offshore, far from the surface roughness transition.

A key input to wake models for wind farms is the wake width behind each turbine, usually denoted by *σ*. While analytical modeling of the wake width is a topic of active research,^{11,12,44} the simplest and most widely used approximation is a linear fit, i.e., $ \sigma = \sigma 0 + k * ( x \u2212 x t )$, where *x* and *x _{t}* are the streamwise coordinate and the location of the turbine, respectively, and the fitting parameters are $ k *$ (wake growth rate) and

*σ*

_{0}(the initial wake width). Several modeling studies

^{4,21,29,42}have focused on understanding the impact of and developing empirical relations for the wake growth rate, $ k *$, for example, as a function of the turbulence intensity upstream of a turbine. The effect of

*σ*

_{0}has been explored to a limited extent, with values proposed by Bastankhah and Porté-Agel

^{3}being used for most studies. In this paper, we explore the impact of

*σ*

_{0}systematically for wind farms sited on a homogeneously rough surface as well as for wind farms sited downstream of an abrupt change in surface roughness.

In this paper, we study the flow in wind farms sited downstream of a rough-to-smooth abrupt surface roughness transition using large-eddy simulations (LES). Numerical simulations are conducted to study the effect of three parameters, namely, the distance between the surface roughness transition and the first turbine row, the surface roughness values (*z*_{01}, *z*_{02}), and their ratio, $ m = z 01 / z 02$. Analytical models for the mean streamwise velocity deficit inside the wind farm are investigated. Specifically, it is of interest to investigate the corrected-Gaussian^{24} single-turbine model and the wake merging methodology of Lanzilao and Meyers,^{25} as they are designed to be applicable to the spatially varying background velocity fields that arise on account of surface roughness heterogeneities. We build upon our previous work, which focused on the effect of an abrupt rough-to-smooth surface roughness jump on the ABL^{14,27} and on the wake of an isolated turbine.^{24} The present study has similar scope as that of Kethavath *et al.*^{24} but for a thirty-turbine wind farm. Additionally, this paper builds upon the work of Goit and Önder,^{18} which did not consider a systematic variation of parameter values and the evaluation of engineering wake models. Furthermore, we also focus on the effect of *σ*_{0} as a key tunable parameter in wind farm wake models and make recommendations for its appropriate values.

The LES methodology, cases simulated, and details of the analytical models evaluated in this study are outlined in Sec. II. Results of the LES are presented and discussed in Sec. III. Evaluation of analytical models to predict the mean velocity deficit in wind farms placed downstream of an abrupt change in surface roughness is discussed in Sec. IV. Summary and conclusions of the work are provided in Sec. V.

## II. NUMERICAL METHODOLOGY

### A. LES framework

Large-eddy simulations (LES) are employed to simulate the flow in wind farms situated over heterogeneously rough surfaces. The incompressible Navier–Stokes equations are solved under the assumption of truly neutral conditions for the atmospheric boundary layer. The influence of the Coriolis force in the atmosphere is not considered in these simulations. A standard concurrent-precursor simulation methodology is utilized, wherein two computational domains are employed.^{15} The precursor domain has a uniform surface roughness *z*_{01} without any turbine sited in it, and the flow here is driven by a constant pressure gradient ( $ \u2212 u * 2 / L z$). Here, $ u *$ is the friction velocity, and *L _{z}* is the domain size in the vertical direction. The main domain includes the features of interest, i.e., surface heterogeneity and/or a wind farm. Each domain has the same size $ ( L x \xd7 L y \xd7 L z )$ along the streamwise (

*x*), spanwise (

*y*), and vertical (

*z*) directions, respectively. The aerodynamic roughness in the main domain is

*z*

_{01}up to $ x s = 10 d$ (where

*d*is the diameter of each turbine) and

*z*

_{02}beyond it in cases involving a surface roughness heterogeneity.

The last quarter of the main domain is referred to as the “fringe” region, where an additional forcing term nudges the flow toward the same conditions as in the precursor domain. The aerodynamic surface roughness in the fringe region is the same as in the upstream and precursor domains, i.e., *z*_{01}. An array of turbines (shown in Fig. 1) are placed in the main domain with three columns and ten rows. The first row of turbines is at a streamwise distance of *x _{t}* from the inlet plane, and the columns are located at 2

*d*, 6

*d*, and 10

*d*in the spanwise direction. The spacings between the wind turbines along the streamwise and spanwise are

*s*and

_{x}*s*, respectively. Turbine row

_{y}*n*is at $ x t , n = x t + ( n \u2212 1 ) s x$, so that $ x t , 1$ is identical to

*x*. We use spacings $ s x = 5 d$ and $ s y = 4 d$, a fixed value of the local thrust coefficient $ C T \u2032 = 0.9$ (nominal thrust coefficient of $ C T = 0.6$), and hub-height $ z t = 0.6 d$. All these values are reasonable approximations of actual wind farms and have been used in previous idealized studies of wind farm wakes.

_{t}^{15,24,39}While not pursued in this paper, sensitivity of the results to these parameters should be studied in future work.

Both the horizontal (*x*, *y*) directions are assumed to be periodic. The top boundary (*z* = *L _{z}*) has free-slip conditions, while the bottom boundary (

*z*= 0) has a shear stress specified via a wall model. No-penetration vertical velocity conditions are applied to both walls.

The simulations are carried out using the in-house “PadeOps-igrid” code,^{41} which employs Fourier-spectral discretization in the horizontal (*x*, *y*) directions and a sixth-order staggered compact finite-difference scheme in the vertical direction. A third-order Runge–Kutta method is employed for the time advancement. All simulations maintain a constant Courant-Friedrichs-Lewy (CFL) value of 0.8, which is within the stability limit of the time stepping scheme. The anisotropic minimum dissipation (AMD) model^{35} is used for subgrid-scale closure. A standard actuator-disk model without rotation (ADM-NR)^{48} and with uniform distribution of the axial thrust force is used to model the turbine forces.^{22} The viscous sub-layer in this flow is too thin to be resolved by moderately sized grids due to the high Reynolds number ( $ R e \u223c 10 8$) based on the boundary layer height, the velocity at the hub-height of a typical wind turbine, and the viscosity of air. Consequently, the region where the direct influence of viscosity is significant is minuscule, and resolving this would require very large grids. This problem is sidestepped by dropping the viscous terms from the Navier–Stokes equations and introducing an additional shear stress at the bottom wall using a wall model. The wall model proposed by Bou-Zeid *et al.*^{5} is employed here. Unlike several other wall models, this wall model does not assume homogeneity in the aerodynamic roughness of the underlying surface, is fully local, and hence, is applicable over heterogeneously rough surfaces.^{5,24,27}

The PadeOps-igrid code employed here has been extensively validated and utilized for various studies previously, including those involving LES of stratified and unstratified turbulent boundary layers, wind turbines, and wind farms.^{15–17,19,20} In our recent previous work,^{27} we have validated the code for a high-Re half channel flow (as a surrogate for an ABL) over an abrupt rough-to-smooth surface roughness transition by comparing with the wind tunnel experimental data of Chamorro and Porté-Agel^{8} and demonstrating good agreement for first-order as well as second-order turbulent statistics. Furthermore, in this paper (in Appendix B), we validate the capability of the code to accurately simulate wind farm wakes by reproducing the experimental data of Chamorro and Porté-Agel^{10} for a miniature-turbine wind farm sited on homogeneously rough ground surface.

### B. Cases simulated

Five sets of LES are conducted, which are presented in Table I. In the first set, illustrated in Fig. 1, case HT1 includes a wind farm with the first turbine at a short distance $ x t \u2212 x s = 4 d$ from the abrupt transition in surface roughness. This case is compared to case HomogT, where the wind farm is retained, but the surface heterogeneity is eliminated. In case HNT (heterogeneous; no turbine), the wind farm is eliminated, but the surface heterogeneity is preserved. The surface roughness values *z*_{01} and *z*_{02} are taken to match the experiments of Chamorro and Porté-Agel,^{8} and the case HNT reproduces this experimental dataset (described in detail in Mondal *et al.*^{27}). The resulting ratio of surface roughnesses, $ z 01 / z 02 = 83.3$, is of the order of magnitude of a transition from smooth land (e.g., grassland) to sea^{8,43} or from rough suburban landscape to bare soil or snowy surfaces. The computational domain size for all the cases mentioned is $ ( L x \xd7 L y \xd7 L z )$ = $ ( 96 d \xd7 12 d \xd7 4.5 d )$. The simulations are carried out on a uniform grid with $ ( n x \xd7 n y \xd7 n z )$ = $ ( 768 \xd7 96 \xd7 108 )$ points in each of the precursor and main domains.

Case . | Description/purpose . | $ ( x s , x t ) / d$ . | $ ( L x \xd7 L y \xd7 L z ) / d$ . | $ ( n x \xd7 n y \xd7 n z )$ . | $ z 01 / d$ . | $ m = z 01 / z 02$ . |
---|---|---|---|---|---|---|

HomogT | Homogeneous; turbine | $ ( \u2212 , 5 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | NA |

HNT | Heterogeneous; no turbine | $ ( 10 , \u2212 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT1 | Heterogeneous; turbine | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT2 | Varying x _{t} | (10, 17) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT3 | (10, 20) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT4 | (10, 30) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-m1 | Varying z_{02} and z_{01} | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 25 |

HT1-m3 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 125 | |

HT4-m1 | (10, 30) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 25 | |

HT4-m3 | (10, 30) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 125 | |

HT1- $ z 01 \u2032$ | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 1.125 \xd7 10 \u2212 3$ | 25 | |

HT1-g1 | Sensitivity to grid and domain size | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 256 \xd7 32 \xd7 30$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT1-g2 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 512 \xd7 64 \xd7 60$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-g3 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 640 \xd7 80 \xd7 90$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-g5 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 1024 \xd7 128 \xd7 120$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-g4 $\u2032$ | (10, 14) | $ 96 \xd7 12 \xd7 6$ | $ 768 \xd7 96 \xd7 144$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HomogE | Validation with experiments^{10} | $ ( \u2212 , 5 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 2.0 \xd7 10 \u2212 4$ | NA |

HomogT-z0 | Upstream TI value matched with HT4 | $ ( \u2212 , 5 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 8.75 \xd7 10 \u2212 4$ | NA |

Case . | Description/purpose . | $ ( x s , x t ) / d$ . | $ ( L x \xd7 L y \xd7 L z ) / d$ . | $ ( n x \xd7 n y \xd7 n z )$ . | $ z 01 / d$ . | $ m = z 01 / z 02$ . |
---|---|---|---|---|---|---|

HomogT | Homogeneous; turbine | $ ( \u2212 , 5 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | NA |

HNT | Heterogeneous; no turbine | $ ( 10 , \u2212 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT1 | Heterogeneous; turbine | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT2 | Varying x _{t} | (10, 17) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT3 | (10, 20) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT4 | (10, 30) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-m1 | Varying z_{02} and z_{01} | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 25 |

HT1-m3 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 125 | |

HT4-m1 | (10, 30) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 25 | |

HT4-m3 | (10, 30) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 3.75 \xd7 10 \u2212 3$ | 125 | |

HT1- $ z 01 \u2032$ | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 1.125 \xd7 10 \u2212 3$ | 25 | |

HT1-g1 | Sensitivity to grid and domain size | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 256 \xd7 32 \xd7 30$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 |

HT1-g2 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 512 \xd7 64 \xd7 60$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-g3 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 640 \xd7 80 \xd7 90$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-g5 | (10, 14) | $ 96 \xd7 12 \xd7 4.5$ | $ 1024 \xd7 128 \xd7 120$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HT1-g4 $\u2032$ | (10, 14) | $ 96 \xd7 12 \xd7 6$ | $ 768 \xd7 96 \xd7 144$ | $ 3.75 \xd7 10 \u2212 3$ | 83.3 | |

HomogE | Validation with experiments^{10} | $ ( \u2212 , 5 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 2.0 \xd7 10 \u2212 4$ | NA |

HomogT-z0 | Upstream TI value matched with HT4 | $ ( \u2212 , 5 )$ | $ 96 \xd7 12 \xd7 4.5$ | $ 768 \xd7 96 \xd7 108$ | $ 8.75 \xd7 10 \u2212 4$ | NA |

In the second set of simulations, the position of the first row in the wind farm ( $ x t = 17 d$, 20*d*, and 30*d*) is changed. These cases enable an assessment of how the position of the turbines relative to the location of the abrupt jump in surface roughness impacts the results.

In the third set, simulations are carried out for different ratios of the surface roughness values, i.e., $ m = z 01 / z 02$. Without changing the upstream surface roughness (*z*_{01}), simulations are conducted for two additional values of *m* = 25 and *m* = 125 (cases HT1-m1 and HT1-m3). In case HT1- $ z 01 \u2032$, we keep the *z*_{02} same as in case HT1, but the ratio is *m* = 25.

In the fourth set, the effect of grid size on the results is evaluated by comparing cases HT1-g1, HT1-g2, HT1-g3, and HT1-g5 to the baseline case HT1. This series of simulations demonstrates that $ 768 \xd7 96 \xd7 108$ points are adequate for obtaining converged results as shown in Appendix A. An additional simulation, case HT1-g4 $\u2032$, is carried out to demonstrate that the vertical domain size is adequately large.

Case HomogE is used to validate our numerical framework for its ability to accurately simulate wake interactions in a wind farm sited on a homogeneously rough ground surface by demonstrating ( Appendix B) good agreement with the experimental results of Chamorro and Porté-Agel.^{10}

Finally, a comparison is made between cases HomogT, HT4, and HomogT-z0, which has a homogeneously rough surface but with *z*_{0} adjusted such that the turbulence intensity upstream of the first turbine matches the value seen in the HT4 case. This comparison is discussed in Appendix C and makes it clear that the effects of an abrupt surface roughness transition upstream of a wind farm are not explained only by a lower upstream turbulence intensity but are qualitatively different.

All quantities are normalized using the velocity scale $ u *$ and the length scale *d*. Typical values for these may be considered to be $ u * = 0.45$ m/s and *d* = 100 m. Simulations are conducted for 180 time units, where one time unit equals $ d / u *$. The flow field is given 60 time units to reach a statistically stationary state, and the remaining 120 time units are utilized to calculate statistical averages. These spin up and averaging durations are consistent with prior studies.^{15,39} The results of the precursor are averaged in both time and the horizontal planes, whereas the statistics for the main domain are averaged in time and along the vertical planes located at turbine mid-span locations ( $ y = 2 d$, 6*d*, and 10*d*) of the three turbine columns.

### C. Analytical models

The analytical model proposed by Bastankhah and Porté-Agel^{2} (called Gaussian model henceforth), which has been widely adopted and successfully used, assumes that the wake velocity deficit decays in the streamwise (*x*) direction and follows a self-similar Gaussian profile in the radial (*y*, *z*) directions. In this paper, we consider two different definitions of velocity deficit. The first is the difference between the incoming upstream velocity and the downstream velocity at any point, i.e., $ D U ( x , y , z ) = U ( x u p , y , z ) \u2212 U ( x , y , z )$. The second is the difference in the velocity at the same downstream location between case HNT and a heterogeneous case (e.g., HTx), i.e., $ \Delta U ( x , y , z ) = U 0 ( x , z ) \u2212 U ( x , y , z )$. The second definition is necessitated because, as shown below, the first definition attains negative values close to the ground for the heterogeneous cases. The second definition of velocity deficit is always positive for all cases and is considered for analytical modeling purpose in this paper.

*x*and

*z*, predicts the velocity deficit due to turbine

*n*at a downstream point (

*x*,

*y*,

*z*) as

*U*

_{0}is a function of

*x*and

*z*, which is not incorporated in the Gaussian model. The varying background velocity field is incorporated into the Gaussian model with slight modifications, similar to what was described in our previous work (Kethavath

*et al.*

^{24}). The velocity deficit due to turbine

*n*is given by

*n*, respectively. $ U 0 , h = U 0 ( x , z t )$ is the velocity at the hub-height in the HNT case, which is obtained from the model described in Ghaisas.

^{14}The function $ F ( x ) \u2248 0.025405 ( x \u2212 x t , n ) / d \u2212 0.060504$ is a correction to the Gaussian model proposed by Kethavath

*et al.*

^{24}Another correction proposed by Kethavath

*et al.*,

^{24}to account for additional terms that are missing from the Gaussian model, gives the non-dimensional constant $ \gamma \u2248 0.01629 \u2212 87.33 \zeta + 7.738 \xd7 10 4 \zeta 2$, where

*ζ*is obtained from the relation

*W*

_{0}is the vertical velocity in the HNT case and is obtained by imposing continuity with the

*U*

_{0}field.

^{24}

^{,}

*δ*is the internal boundary layer (IBL) height, calculated using the empirical relation

_{i}^{13}$ \delta i ( x ) = z 02 [ 0.75 \u2212 0.03 \u2009 ln \u2009 ( z 02 / z 01 ) ] ( ( x \u2212 x s ) / z 02 ) 0.8$. The last term in Eq. (2) is added only up to the top tip-height, i.e., $ z \u2264 z t + d / 2$.

For a homogeneous surface, the last term is always zero. Furthermore, if we consider the corrections *γ* and *F*(*x*) as zero, Eq. (2) reduces to the Gaussian model [i.e., Eq. (1)]. We refer to Eq. (2) as the Gaussian-corrected (GC) model.

^{26,29}

^{23,47}

^{25}is also considered here, specifically because it accounts for a heterogeneous background velocity field and is expected to be applicable to flows involving surface roughness heterogeneities. The composite wake is given by

## III. RESULTS

### A. Effect of heterogeneity on a wind farm

Cases HomogT, HNT, and HT1 are discussed first. This allows a comparison of the performance of a wind farm sited behind a surface roughness jump with that of a wind farm sited on homogeneously rough ground surface. The results of the precursor simulations and the profiles upstream of the first-row turbines from the main domain (i.e., $ x / d = x t , 1 \u2212 1$) for the cases HomogT, HNT, and HT1 are shown in Fig. 2. In the HNT case, where there is no wind farm and hence *x _{t}* is not defined, the profiles are plotted at the same location as in the HT1 case. Profiles from the left to right are the time-averaged streamwise velocity (

*U*), streamwise turbulence intensity (

*TI*), and negative total shear stress (−

*TSS*). The streamwise turbulence intensity is defined as $ T I = U \u2032 U \u2032 \xaf / U h$, where

*U*is the velocity at the hub-height 1

_{h}*d*upstream of the first turbine row. The total shear stress is the sum of resolved, subgrid-scale and wall model components, $ TSS = U \u2032 W \u2032 \xaf + \tau 13 sgs \xaf + \tau 13 w m \xaf$. The profiles of

*U*,

*TI*, and

*TSS*vary by 0.3%, 1.7%, and 5.9%, respectively, between the precursor and the main domain for the HomogT case. This close agreement is expected since the aerodynamic surface roughness is the same $ ( z 01 )$ in the precursor domain and in the main domain in the HomogT case. The small differences can be attributed to induction effects caused by the wind farm. For the cases HT1 and HNT, close to the lower wall, there is an increase in

*U*and a decrease in both

*TI*and

*TSS*compared to the HomogT case. This is because the surface roughness has changed from

*z*

_{01}to

*z*

_{02}at the location 1

*d*upstream of the first turbine row. The profiles of

*U*,

*TI*, and

*TSS*match closely with the precursor and HomogT profiles above a certain height, as is also expected, highlighting the well-controlled nature of the inflow conditions applied in our numerical simulations.

#### 1. Mean velocity

Contour of the instantaneous streamwise velocity at both x–z and x–y planes for the case HT1 are shown in Fig. 3. The surface roughness change is marked by the dashed vertical line, and the location of the wind turbines is marked by the solid white line. As a result, wind turbine wakes are seen and identified by the low-speed regions downstream of each turbine. Additionally, vertical profiles of the time-averaged streamwise velocity (*U*) at mid-span and averaged over the three columns of turbines obtained from the LES are shown in Fig. 4. The results of cases HomogT, HNT, and HT1 are shown at different locations downstream of three rows (i.e., row 1, row 5, and row 9). In the cases with wind turbines, HomogT and HT1, a low-speed region downstream of each turbine is seen, which indicates the turbine wake. In cases HNT and HT1, the surface roughness value changes from *z*_{01} (rough) to *z*_{02} (smooth) at $ x s = 10 d$. After this step change, an increased streamwise velocity near the bottom wall is seen at all downstream locations compared to case HomogT. As seen in Fig. 4, close to the bottom wall, a slight reduction in the streamwise velocity is seen for case HT1 as compared to case HNT.

#### 2. Turbulence intensity

Vertical profiles of the streamwise turbulence intensity (*TI*) at the turbine center plane, averaged over the three turbine columns, are shown in Fig. 5. It is defined as the ratio of the standard deviation of the streamwise velocity fluctuations to the mean velocity at the turbine hub-height (i.e., $ T I = U \u2032 U \u2032 \xaf / U h$). The results are shown at four downwind locations $ ( x \u2212 x t , n ) / d$ = 1–4 behind the same three rows as previously, i.e., row 1, row 5, and row 9. Under both, HomogT and HT1 conditions, a significant increase in turbulence intensity is observed behind each turbine and between $ z = z t \xb1 d / 2$. At each downstream location, the maximum of the *TI* profile occurs at the top tip height because of the strong shear and correspondingly large turbulent kinetic energy (TKE) production at that location.

In cases HNT and HT1, the impact of the surface roughness transition is evident, especially close to the lower wall. The *TI* values are lower in these cases compared to HomogT, where a higher surface roughness exists throughout the domain. The *TI* profiles obtained from cases HomogT and HT1 are identical behind row 1 except very close to the bottom wall. Behind rows 5 and 9, the turbulence profiles are appreciably smaller for case HT1 than for case HomogT over a larger portion of the domain in the vertical direction. Thus, it is clear that the effect of reduced surface roughness is limited to the region close to the ground surface ( $ z / d \u2272 0.8$) at smaller downstream distances from the step change, but it increasingly affects the entire boundary layer at larger downstream distances.

#### 3. Total shear stress

Vertical profiles of the normalized total negative shear stress (*TSS*), averaged over the three turbine column midplanes, obtained from the LES, are shown in Fig. 6. The total shear stress is given by the sum of the resolved, subgrid-scale and wall model contributions, $ TSS = ( U \u2032 W \u2032 \xaf + \tau 13 sgs \xaf + \tau 13 w m \xaf )$. The surface roughness transition from *z*_{01} to a lower *z*_{02} in cases HNT and HT1 results in reduced values of *TSS* near the bottom wall. In cases HomogT and HT1, the effects of turbines are clearly seen, with an increased −*TSS* above the turbine hub-height and reduced −*TSS* below the hub-height in the wake region (roughly $ z t \xb1 d / 2$). Consistent with the profiles of *TI*, the shear stress magnitudes are smaller for the HT1 case than for the HomogT case. These effects are confined to only close to the bottom wall ( $ z / d \u2272 0.8$) behind row 1, to roughly $ z / d \u2272 2.5$ behind row 5 and to about $ z / d \u2272 3.5$ behind row 9.

#### 4. Internal boundary layer height

The profiles of *TI* and *TSS* shown earlier deviate from each other among the cases HomogT, HNT, and HT1 and from the precursor plots at a certain height. These observations indicate the formation of an internal boundary layer that progressively increases in height with increasing streamwise distance. The evolution of the IBL height is shown for the HomogT, HNT, and HT1 cases in Fig. 7. The IBL height at each *x* location is calculated as the highest position where the upstream *TSS* profile and the *TSS* profile at *x* differ by more than 10%. Calculating the IBL height based on other profiles, e.g., *TI*, gives similar values.^{27} The *x* coordinate for the HomogT case is translated such that the position of the first turbine row matches between the HomogT and HT1 cases.

The height of the IBL formed due to the surface roughness jump in the HNT case is zero upstream of the step change and increases steadily with increasing $ x adj = x \u2212 x s$. The height of the IBL formed due to the wind farm in the HomogT case is zero until the turbine location, jumps abruptly to $ \delta i = 1.2 d$ at the turbine location, and then increases gradually. For the HT1 case, the IBL is formed due to the surface heterogeneity at *x _{adj}* = 0 and due to the wind farm starting from $ x adj = 4 d$. The IBL height follows the curve for the HNT case until the first turbine row location and closely follows the HomogT curve after the first row of turbines is encountered. These results indicate that the IBL height due to the wind farm is not affected by the presence of the roughness heterogeneity of the underlying surface.

#### 5. Velocity deficit

In Fig. 8, a comparison of three velocity deficit profiles, evaluated from the time-averaged LES results, is presented. The conventional definition of velocity deficit at any point downstream of a turbine is the difference between the incoming upstream velocity and the velocity at that point, i.e., $ D U ( x , y , z ) = U ( x u p , y , z ) \u2212 U ( x , y , z )$. Here, we take *x _{up}* to be equal to $ x t \u2212 1 d$. In case HomogT, close to the bottom wall,

*DU*has positive values at all streamwise locations. The velocity deficit peaks at the hub-height and decreases gradually to zero above and below this height. In contrast, for case HT1, the

*DU*profile shows significant negative values close to the lower boundary. This is due to the acceleration of the flow downstream of $ x s = 10 d$ due to the surface roughness change, which results in $ U ( x , y , z )$ being larger than $ U ( x u p , y , z )$.

An alternate definition of velocity deficit is considered, defined as the difference between the velocity at the same downstream locations in cases HNT and HT1, i.e., $ \Delta U ( x , y , z ) = U 0 ( x , z ) \u2212 U ( x , y , z )$. Here, *U*_{0} refers to the velocity field obtained from the simulation with heterogeneity but no turbine (case HNT). This velocity field is spatially averaged in the spanwise (*y*) direction due to homogeneity in that direction. The profiles of $ \Delta U$ shown in Fig. 8 exhibit non-negative values near the bottom wall and are similar to the *DU* profiles obtained from the HomogT case. The advantage of the alternate definition of velocity deficit is that, in the region $ 0 \u2264 z \u2264 2 d$ that is directly affected by the turbine wakes, it has positive values, unlike the *DU* profile. As a result, unlike the *DU* profile, the $ \Delta U$ profile can be modeled by a function that is always positive (e.g., a Gaussian).

### B. Effect of wind farm position

The effect of the distance between the surface roughness jump (*x _{s}*) and the first turbine row in a wind farm (

*x*) is studied in this section. The qualitative flow features outlined for $ x t = 14 d$ (Case HT1) in Subsection III A are seen in cases HT2, HT3, and HT4 as well, where the distance

_{t}*x*is progressively increased to 17

_{t}*d*, 20

*d*, and 30

*d*, respectively. These results are not repeated here for brevity.

The effect of increasing *x _{t}* on the IBL height is shown in Fig. 9. The IBL height increases gradually from the step change location ( $ x adj = x \u2212 x s = 0$) until the position of the first turbine row, as seen in Fig. 9(a). Once the first turbine row is encountered, the IBL height abruptly jumps to $ \delta i = 1.2 d$ for all cases. The magnitude of this abrupt jump reduces across the HT1 through HT4 cases because the IBL upstream of

*x*, due to the surface heterogeneity alone, grows to a larger height as

_{t}*x*increases. Figure 9(b) shows that the IBL heights downstream of the first turbine row are almost the same for all cases, which indicates that they are not affected by the underlying surface roughness.

_{t}The effect of increasing *x _{t}* is quantified in terms of the power in Fig. 10. Because all turbines have the same rotor diameter and the same local thrust coefficient, the power of each turbine is calculated as $ P n = U n 3 \xaf$, where

*U*represents streamwise velocity averaged over the rotor disk of turbine

_{n}*n*. The overhead bar denotes averaging over time and turbine columns, and the subscript

*n*represents the turbine row within the wind farm. While this power is normalized by $ u * 3$, we refer to it as “normalized absolute” power to distinguish it from the relative power. The relative power for turbine row

*n*is computed by normalizing by the power of the first-row turbine, as $ P n / P 1$.

The normalized absolute power [Fig. 10(a)] of first-row turbines increases with increasing *x _{t}* because of the acceleration of the flow due to the step change in surface roughness leading to an increased disk-averaged velocity. The first-row turbine powers are 0.7%, 3.5%, 5.5%, and 13.4% higher in the cases HT1, HT2, HT3, and HT4, respectively, than in the HomogT case. In the cases HT1 through HT4, the powers generated by the rest of the turbines (rows 2 through 10) are larger compared to case HomogT. The powers

*P*

_{2}through

*P*

_{10}are slightly different between the HT cases with different

*x*. This indicates a complicated interaction between turbine wakes and the acceleration of the mean flow due to the surface roughness jump.

_{t}The relative powers of the downstream turbines (rows 2 through 10), shown in Fig. 10(b), are smallest for case HomogT. Changing to a smooth surface reduces the wake losses and increases the relative powers of the downstream turbines. As *x _{t}* increases, the relative powers start reducing, which suggests that as the distance from the step change to the first row becomes large, the wake losses again approach the same level as for a homogeneous surface. This is quantitatively shown in Fig. 11, where the average of the relative powers of the downstream turbines is plotted for different cases. Compared to the HomogT case, the relative powers in the HTx cases are larger by about 19.3% (for HT1) to about 10.8% (for HT4). These results show that for a wind farm sited on a heterogeneous surface, the distance between the surface roughness transition and the first turbine row ( $ x t \u2212 x s$) is a crucial parameter, and that an optimal value of $ ( x t \u2212 x s )$ exists, for which the wake losses are minimized.

### C. Effect of surface roughness ratio

The normalized absolute and relative powers for different $ m = z 01 / z 02$ values with the first turbine row at *x _{t}* = 14 (corresponding to case HT1) are shown in Fig. 12. First, we consider a fixed $ z 01 = 3.75 \xd7 10 \u2212 3 d$ and varying

*z*

_{02}to obtain different

*m*ratios. The values of

*z*

_{02}are $ 1.5 \xd7 10 \u2212 4 d , \u2009 4.5 \xd7 10 \u2212 5 d$, and $ 3 \xd7 10 \u2212 5 d$ for the

*m*ratios 25, 83.3, and 125, respectively. The normalized absolute and relative powers of the HT1-m1 and HT1-m3 cases are almost the same as for the case HT1. Specifically, Fig. 12(a) demonstrates that the normalized absolute power of the front-row turbine is nearly equal under all conditions. The same is found to be true between cases HT4-m1, HT4, and HT4-m3, shown in Fig. 13. Combined with Fig. 10(a), it is clear that, within the parameter ranges studied, the power of the front-row turbine and the wake losses in the wind farm are dependent on the distance from the heterogeneity, but not strongly dependent on the strength of the heterogeneous roughness jump, as long as the upstream surface roughness is kept the same. It is possible that this conclusion does not hold for much longer (more than ten rows) extended wind farms, which can be studied in the future.

Figure 14 shows the results of case HT1- $ z 01 \u2032$, where the ratio *m* = 25 is obtained by keeping *z*_{02} unchanged from case HT1 but by changing $ z 01 = 1.125 \xd7 10 \u2212 3 d$. For case HT1- $ z 01 \u2032$, the front-row turbine power is much larger than that in cases HT1 and HT1-m1, which is directly related to the reduced upstream roughness *z*_{01}. Figure 14(b) shows that a significant difference exists in the relative powers of the downstream turbines for case HT1- $ z 01 \u2032$ compared to HT1 and HT1-m1.

To summarize, the surface roughness *m* affects the performance of a wind farm if the upstream roughness (*z*_{01}) is changed keeping the downstream roughness (*z*_{02}) fixed but not if the downstream roughness is changed keeping the upstream roughness fixed.

## IV. EVALUATION OF ANALYTICAL WAKE MODELS

The ability of different analytical models to predict the flow features described earlier is discussed in this section. We focus on describing the HTx cases for brevity. The other cases are predicted with a similar accuracy and are not discussed here.

### A. Front-row turbines

The performance of the front-row turbines is considered separately from that of the subsequent, downstream, rows. In the HTx cases, the power generated by the most upstream turbine increases as the distance from the surface roughness jump is increased, as shown in Fig. 10(a). This increase is directly linked to the increased wind speed at the hub-height due to the abrupt change in surface roughness. The velocity field over a heterogeneous surface in the absence of turbines, $ U 0 ( x , z )$, is predicted using the model of Ghaisas.^{14} This predicted velocity field agrees well with the LES results, as seen in Fig. 15.

### B. Wake growth rate

Modeling the wake effects of upstream turbines is critical to accurately modeling the performance of the downstream turbines. A crucial input to wake models is the wake growth rate behind each turbine row.

The evolution of the wake widths along the streamwise direction behind each turbine row for the HomogT and HT1 cases is shown in Fig. 16. At each *x*, the width is obtained as the geometric mean ( $ \sigma n = \sigma z \sigma y$) of the wake widths along the spanwise (*σ _{y}*) and vertical (

*σ*) directions. The

_{z}*σ*and

_{y}*σ*are computed from the velocity deficit profiles along

_{z}*y*and

*z*, respectively, using the expressions provided in Bastankhah and Porté-Agel.

^{3}

Figure 16 shows that the wake width evolution is similar for the HomogT and the HT1 case. In each case, the first row turbine wake width is the smallest, and the wake widths increase monotonically from row 1 to row 10. This is expected since the turbulence intensity is the smallest upstream of the first row and increases as one proceeds into the wind farm. Behind each turbine row, the wake has an almost linear expansion, similar to previous observations.^{2,3,15} A linear fit, $ \sigma n \u2248 \sigma 0 , n + k n * ( x \u2212 x t , n )$, gives the growth rate, $ k n *$, and the intercept (initial wake width), $ \sigma 0 , n$, for each *n*. The wake growth rates extracted from the LES for the HomogT and the HTx cases are shown in Fig. 17(a). It is clear that the behavior is qualitatively similar for the HomogT and all the HTx cases with small quantitative differences between the HTx cases.

Figure 16 further shows that, for a given turbine row, the *σ _{n}* values are larger for the HT1 case than for the HomogT case. The reason for this difference can be traced to the fact that the deficit $ \Delta U$ for the HT1 case is larger than the deficit

*DU*for the HomogT case, particularly close to the ground below the turbine hub-height, as seen in Fig. 8. The difference in the background velocity, namely, the incoming upstream velocity profile for HomogT and the spatially evolving profile from the HNT case for the HT1 case, lead to the deficits

*DU*and $ \Delta U$ being different, which in turn leads to differences in the

*σ*and $ k *$ values between the HNT and HT1 cases. The same observations hold for the other HTx cases and explain the larger $ k *$ values seen in Fig. 17(a) for the HTx cases over the HomogT case.

_{n}The wake growth rate is usually found to be correlated with the streamwise turbulence intensity.^{15,29} This does not hold directly for HTx cases since the $ k *$ values are larger [Fig. 17(a)], whereas the *TI* values are smaller (Fig. 5) in the HTx cases compared to the HomogT case, particularly for downstream rows. The wake growth rate is seen in Fig. 17(b) to correlate well (correlation coefficient of 0.9) with the added turbulence intensity, $ \Delta T I = T I 2 \u2212 T I 0 2$. Here, *TI*_{0} is the turbulence intensity in the HNT case and, in practical applications, can be obtained as a function of *x* and *z* from the model described by Mondal *et al.*^{27} for cases involving an abrupt surface roughness jump. For wind farms on homogeneously rough surfaces, *TI*_{0} profile is the same at all streamwise locations, equal to the profile upstream of the wind farm, and the $ \Delta T I$ is the same as the added turbulence intensity used in previous studies.^{4,29}

### C. Downstream turbines

Two single-turbine wake models (Gaussian and Gaussian-corrected) coupled with three wake merging methods (linear, quadratic, and Lanzilao and Meyers), described in Sec. II C, are compared with the LES data. The normalized velocity deficits at the turbine hub-height and in the turbine center-plane along the spanwise direction, i.e., $ \Delta U ( x , y t , z t ) / U 0 , h , u p$ for the HomogT and HT1 cases, are shown in Figs. 18 and 19, respectively. To quantitatively assess the accuracy of the different methods, the mean absolute error between the LES and the model predictions is calculated and shown in Figs. 18(d)–18(f) and 19(d)–19(f). The $ k *$ values are those extracted directly from the LES, while different choices are considered for *σ*_{0}. The lines labeled as “Act *σ*_{0}” use the actual values extracted from the LES. The lines labeled as “Fixed *σ*_{0}” use a uniform, fixed, value for each turbine, varying between 0.34 and 0.44 in steps of 0.01. The large thickness of these lines captures the upper and lower bounds of the velocity deficit predictions as *σ*_{0} is varied over this range.

Several observations can be made from Figs. 18 and 19. First, for both (HomogT and HT1) cases, linear wake merging generally overestimates the velocity deficits [Figs. 18(a) and 19(a)]. On the other hand, quadratic wake merging generally underestimates the deficits [Figs. 18(b) and 19(b)]. The Lanzilao and Meyers wake merging methodology gives the closest predictions for both cases, as shown in Figs. 18(c) and 19(c).

The results of linear and quadratic wake merging methodologies are very sensitive to *σ*_{0}, when the same value is used for each turbine wake. In the HomogT case, the errors of the GL and GCL models vary from roughly 7% to 19% and from $ 5.5 %$ to 10%, respectively [Figs. 18(d) and 19(d)]. In the HT1 case, the errors similarly vary from around 5% to about 12%. Considering the GQ and GCQ models, the errors in predictions of the HomogT [Fig. 18(e)] and HT1 [Fig. 19(e)] cases reach as high as 20% and 16%, respectively. The Lanzilao and Meyers merging is not as sensitive to *σ*_{0} as seen quantitatively [Figs. 18(f) and 19(f)], with the errors varying between roughly 5% and 10% for the HomogT and HT1 cases.

The sensitivity of the errors to *σ*_{0} is primarily because of the sensitivity of the predictions of the first few (up to 5) turbines in each wind farm. By the seventh row ( $ x / d = 35$), the model predictions are almost the same regardless of *σ*_{0}. This is seen from the vertical spread in the green and blue solid lines in Figs. 18(a)–18(c) and 19(a)–19(c).

In most cases, the Gaussian-corrected wake model gives a smaller error and closer prediction of the LES results than the corresponding Gaussian wake model. This is particularly true in the HomogT case for linear and Lanzilao and Meyers merging [Figs. 18(d) and 18(f)] and in the HT1 case for the quadratic and Lanzilao and Meyers merging [Figs. 18(e) and 18(f)].

Carefully evaluating all model combinations, it is observed that the GCLM method with a fixed, uniform, *σ*_{0} set to any value between 0.38 and 0.42 consistently gives small errors of around 5%–6% across all cases. Thus, it is possible to make a clear recommendation of using the corrected Gaussian model with the Lanzilao and Meyers merging method and a fixed uniform value of *σ*_{0}. The mean absolute errors for the other three HTx cases are compiled in Table II. It is clear that the recommended method gives accurate predictions, with errors being consistently around 5%–6%.

Fixed $ \sigma 0 = 0.42$ . | HT2 . | HT3 . | HT4 . | |||
---|---|---|---|---|---|---|

G . | GC . | G . | GC . | G . | GC . | |

Linear | 4.8 | 5.8 | 4.9 | 6.5 | 5.1 | 9.9 |

Quadratic | 17.6 | 10 | 17.6 | 10.1 | 17.4 | 7.9 |

LM | 6.4 | 5.1 | 6.3 | 5.4 | 6.2 | 6.3 |

Fixed $ \sigma 0 = 0.42$ . | HT2 . | HT3 . | HT4 . | |||
---|---|---|---|---|---|---|

G . | GC . | G . | GC . | G . | GC . | |

Linear | 4.8 | 5.8 | 4.9 | 6.5 | 5.1 | 9.9 |

Quadratic | 17.6 | 10 | 17.6 | 10.1 | 17.4 | 7.9 |

LM | 6.4 | 5.1 | 6.3 | 5.4 | 6.2 | 6.3 |

## V. CONCLUSIONS

The evolution of the wake of a wind farm sited on a surface undergoing an abrupt change in aerodynamic roughness from a relatively high value to a relatively low value is studied using large-eddy simulations. The flow in the heterogeneous case is compared to the flow over a homogeneous surface with a wind farm (HomogT) on the one hand and to the flow over a heterogeneous surface without wind farm (HNT) on the other. The change in surface roughness leads to significant differences in the different turbulent quantities. Specifically, the mean velocity is larger, while the magnitudes of the turbulence intensity and the vertical shear stress are smaller in wind farms sited downstream of an abrupt rough-to-smooth roughness transition as compared to a wind farm sited on a homogeneously rough surface. The vertical extent over which these differences are significant increases with increasing downstream distance from the roughness jump. An internal boundary layer develops starting from the position of surface roughness jump and affects the turbulent statistics in the wake of the wind farm.

The velocity deficit (denoted *DU*), when defined in the usual manner as the difference between the velocity upstream of and inside the wind farm, attains negative values close to the ground for a wind farm sited downstream of an abrupt change in surface roughness. An alternate definition (denoted $ \Delta U$), namely, the difference between velocities on a heterogeneously rough surface in the absence of and in the presence of the wind farm, gives positive values in the rotor region. This alternate definition of velocity deficit is more beneficial since it can be modeled as, e.g., a Gaussian in the radial direction, which is always positive. The $ \Delta U$ values are found to be larger than the *DU* values due to differences in the “background” velocity field. As a result, the wake width values (*σ*) and the wake growth rates ( $ k *$) calculated from wind farms sited on homogeneously rough surfaces and surfaces with a roughness jump are different, with the values being larger in the latter cases.

The position of the wind farm with respect to the step change in surface roughness has a significant impact on the wake losses and the power generated. The power generated by the front-row turbines changes by 0.7%–13.4%, while the normalized powers of the downstream turbines, which are indicative of the wake losses, change by about 10.8%–19.3% over the cases simulated. Placing the first row of the wind farm 4 diameters downstream of the surface roughness jump is found to be optimal. As the distance to the first turbine row increases, the wake losses of the downstream turbines systematically increase. For very large distances, the wake losses of the downstream turbines approach the same level as in the HomogT case.

Over the range of parameters simulated, changing the aerodynamic roughness height of the downstream (smooth) surface does not affect the wake losses. On the other hand, changing the roughness of the upstream (rough) surface affects the relative power production by the turbines and the wake losses.

Finally, a detailed description of analytical modeling considerations for wind farms sited on homogeneously rough surfaces and on surfaces with an abrupt roughness transition is presented. First, the power of the upstream turbines for different cases is directly related to the velocity seen by the first row turbines. This velocity is shown to be predicted well by our previously developed model (Ghaisas^{14}). Second, the wind turbine wake growth rates in wind farms sited behind a surface roughness jump are shown to be correlated with the added turbulence intensity, $ T I 2 \u2212 T I 0 2$. Here, *TI*_{0} refers to the turbulence intensity of the background flow due to only the step change in surface roughness (the HNT case) and can be prescribed by the simple model described in Mondal *et al.*^{27}

Third, different analytical modeling choices are evaluated for their ability to predict the wake losses in wind farms sited downstream of an abrupt surface roughness jump. Two models for the wake of an isolated turbine are coupled with three distinct wake superposition methods. Results of linear and quadratic wake merging strategies are very sensitive to *σ*_{0}, which can be considered to be a tunable parameter in the model. The Lanzilao and Meyers wake merging strategy is less sensitive to *σ*_{0}. A value between 0.38 and 0.42 is found to be suitable for the two single-turbine models coupled with either the linear or the Lanzilao and Meyers merging method. Overall, the corrected Gaussian wake model with the Lanzilao and Meyers wake merging methodology gives good agreement, within 5%–6% of the LES data for all cases studied here.

Future work can be focused on studying a wider range of parameter values, such as other turbine operating conditions (larger thrust coefficients), streamwise and spanwise spacings, the ratio of the turbine hub-height to the diameter, and alternate arrangements such as partial or full staggering. This study can be extended to smooth-to-rough roughness transitions as well as other canonical configurations of surface roughness heterogeneities, such as stripes of finite widths, oblique stripes, and rectangular patches with varying aerodynamic roughness values.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the SERB Core Research Grant (CRG/2022/006735) for enabling this work as well as the computational resources made available under the National Supercomputing Mission (Grant DST/NSM/R&D_HPC_Applications/2021/28) on Param-Brahma cluster at IISER Pune and Param-Seva cluster at IIT Hyderabad.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Naveen N. Kethavath:** Conceptualization (equal); Data curation (lead); Investigation (lead); Writing – original draft (lead); Writing – review & editing (equal). **Niranjan S. Ghaisas:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

### APPENDIX A: GRID CONVERGENCE AND SENSITIVITY TO DOMAIN SIZE

In this section, we study the sensitivity of the results to the grid size. The HT1 simulation (involving the surface roughness heterogeneity as well as the wind farm) is performed on a series of grids. The domain size is $ ( L x \xd7 L y \xd7 L z )$ = $ ( 96 d \xd7 12 d \xd7 4.5 d )$. The aerodynamic surface roughness jumps from $ z 01 = 3.75 \xd7 10 \u2212 3 d$ to $ z 02 = 4.5 \xd7 10 \u2212 5 d$ at $ x s = 10 d$ and the first turbine row is at $ x t = 14 d$ from the left end of the main domain.

Profiles of the time-averaged streamwise velocity (*U*), the streamwise turbulence intensity ( $ T I = U \u2032 U \u2032 \xaf / U h$), and total negative shear stress (*TSS*) are shown in Figs. 20–22, respectively. The results are also averaged over the three turbine columns. Vertical profiles at four downstream locations behind three rows (row 1, row 5, and row 9) are shown. Except for the coarsest mesh (case HT1-g1), the results are seen to be largely insensitive to the grid resolution with minimal differences between the two finest grids (HT1 and HT1-g5).

To further quantify the sensitivity to the grid, the time-disk and rotor-disk averaged streamwise velocity (*U*) and the streamwise turbulence intensity ( $ T I = U \u2032 U \u2032 \xaf / U h$) as a function of turbine row for the five grid spacings are shown in Figs. 23(a) and 23(b), respectively. Richardson extrapolated (RE) values are determined through a mixed-order analysis,^{34,46} utilizing HT1-g1, HT1-g2, and HT1-g5 cases. Similar to the observations in Figs. 20–22, the profiles corresponding to the cases HT1-g2, HT1-g3, HT1, and HT1-g5 are fairly close to each other and to the profiles obtained by Richardson extrapolation. The relative errors in the mean velocity and *TI*, assuming the Richardson extrapolated profiles to be the exact solution, as a function of turbine row are shown in Figs. 24(a) and 24(b). For the HT1 grid spacing, the mean absolute errors are found to be 1.7% and 6.2%, respectively, for the normalized mean velocity and the turbulence intensity. Further studies for all other cases are carried out using the same resolution as in the HT1 case, i.e., with $ ( n x \xd7 n y \xd7 n z )$ = $ ( 768 \xd7 96 \xd7 108 )$ grid points.

The sensitivity to the domain size is shown in Fig. 25. The HT1-g4 $\u2032$ case employs a domain that is larger in the vertical direction. To keep the same grid spacing as in the HT1 case, the number of points, *nz*, is increased with respect to the HT1 case. Profiles at 3*d* downstream of the three turbine rows (row 1, row 5, and row 9) are seen to be largely insensitive to the vertical domain size. The vertical domain size is also in accordance with that used by previous simulations.^{39}

### APPENDIX B: VALIDATION WITH EXPERIMENT

The *U* and *TI* profiles (Figs. 26 and 27, respectively) obtained from the LES “HomogE” are validated with experimental data of Chamorro and Porté-Agel.^{10} Here, the boundary layer depth is $\delta =4.5d$, the uniform aerodynamic roughness height is $ z 0 / \delta = 4.444 \xd7 10 \u2212 5$, and the local thrust coefficient ( $ C T \u2032$) values used for the different rows in the wind farm are consistent with those reported in Stevens *et al.*^{39} The LES results at different distances downstream of only four rows (row 1, row 4, row 7, and row 10) are shown for brevity, and the results are similar behind other turbine rows. The mean velocities and turbulence intensity values agree well with the experiments. Some differences in the mean velocity (Fig. 26) are seen up to 3*d* downstream of the first-row turbine and at 1*d* downstream of the other rows. The *TI* profiles (Fig. 27) are underpredicted behind the first-row turbine but are well predicted behind the other downstream turbine rows. These differences between the LES and the experiments may be attributed to the turbine support structure utilized in the wind tunnel measurements, which are absent from the LES. A similar level of agreement is observed in previous LES (without the turbine support structures) of this wind farm reported by Stevens *et al.*,^{39} which are also shown in Fig. 26, thus building confidence in our numerical methodology.

### APPENDIX C: COMPARISON WITH HOMOGENEOUSLY ROUGH SURFACE WITH MATCHING *TI*

To ensure that the differences observed in this paper are actually caused by the changed surface roughness, rather than a different ambient turbulence intensity, we compare cases HT4, HomogT, and an additional case HomogT-z0. This last case considers a wind farm with a homogeneously rough surface with the *z*_{0} tuned such that the ambient turbulence intensity upstream of the first turbine matches that in the HT4 case. The ambient rotor-averaged *TI* values (i.e., 1*d* upstream of the first turbine row) are 0.16 for HomogT and 0.12 for HT4 and HomogT-z0. In the HomogT-z0 case, this *TI* is obtained by setting $ z 0 = 8.75 \xd7 10 \u2212 4 d$. Case HT4 is chosen for this exercise rather than the cases HT1, HT2, or HT3 because the *TI* values 1*d* upstream of the first turbine differ the most between cases HomogT and HT4.

Figure 28 shows the rotor-averaged *TI* values 1*d* upstream of each turbine row, the normalized absolute and relative powers, and the wake growth rate ( $ k *$) values at the turbine locations as a function of turbine row. The *TI* values upstream of the first turbine row [Fig. 28(a)] are identical between the HT4 and HomogT-z0 cases by design. However, the *TI* values upstream of the rest of the turbines differ slightly between the HT4 and HomogT-z0 cases. This is expected since, in the absence of the wind farm, the *TI* values are not expected to change with *x* over a homogeneously rough surface but are expected to decrease behind a rough-to-smooth transition. Figure 28(b) shows that decreasing the *z*_{0} leads to an increased absolute power for the HomogT-z0 case compared to the HomogT case. The absolute powers in the HT4 case are also very different compared to the HomogT-z0 case because of the differences in the upstream surface roughness values. The relative powers, shown in Fig. 28(c), are quite different between the HomogT-z0 and the HT4 cases. These differences are consistent with the wake growth rate values shown in Fig. 28(d). For the HomogT-z0 case, the wake recovery is slower, as characterized by the smaller $ k *$ values, which leads to larger wake losses than in the HT4 case.

These results make it clear that the upstream ambient *TI* is not sufficient to characterize the wake losses in a wind farm sited behind a surface roughness heterogeneity. The streamwise variation of the background (due to the surface heterogeneity and in the absence of wind farm) velocity and turbulence intensity fields is crucial for explaining wake losses in wind farms sited behind a rough-to-smooth surface roughness heterogeneity.

## REFERENCES

*European Wind Energy Association Conference and Exhibition*