Model intercomparison of the ABL, turbines, and wakes within the AWAKEN wind farms under neutral stability conditions Open Access

In a previous study,^7,12 a series of numerical simulations of wind farms using different model fidelities and for different atmospheric stability conditions were performed as a part of the American WAKE ExperimeNt (AWAKEN). The simulations included FLOw Redirection and Induction in Steady State (FLORIS) wake models, a number of large-eddy simulation (LES) microscale Adaptive Mesh Refinement (AMR-Wind) and Nalu-Wind runs, as well as idealized and complex terrain Weather Research and Forecasting (WRF) runs. The scale of these simulations ranged from small subsets of the individual wind farms to immense computations using the AMR-Wind LES solver to simulate a 100 × 100 km² domain containing all 541 turbines among the five wind farms (King Plains, Chisholm View, Thunder Ranch, Armadillo Flats, and Breckinridge). Results of these simulations have been previously used to investigate blockage effects around wind farms^6,11,12 and provide uncertainty quantification of the dual-Doppler radar system being leveraged for large-scale flow measurements in AWAKEN.⁹ 

The primary objectives behind the extensive simulation work of the AWAKEN project relate to the investigation of seven broad testable hypotheses that are a priority to the wind research community and the wind energy industry. In the previous work of Cheung et al.,¹² simulations primarily focused on the subset of the following three hypotheses below:^33,34

1. Wind farm wakes propagate on land for tens of kilometers and lower the energy production of neighboring wind farms. Characteristics [magnitude and extent of momentum deficits, magnitude and extent of region of increased turbulent kinetic energy (TKE)] of wind plant wakes depend primarily on the spacing of turbines in a wind farm along the primary wind direction, turbine size, individual turbine power level and hub-height TKE, wind speed at hub height, and atmospheric stability. Wind farm wakes can be steered using coordinated individual turbine yaw control, although topography and yaw-misalignment will also influence wake propagation.

2. Wind turbines in the interior of wind farms tend to have more turbulent inflows resulting in higher damage-equivalent loads than those on the exterior. The turbulence levels in land based wind farms asymptote to a fully developed condition after the first three rows of wind turbines.

3. The maximum energy produced by a large (>100 MW) wind farm is constrained by the momentum flux between the surrounding atmosphere and the wind farm.

All simulations for the AWAKEN project were conducted under stability conditions that are representative of the regional climatology at the AWAKEN site. Previously, idealized stable and unstable conditions were defined using experimental data taken from the Atmospheric Radiation Measurement (ARM) site located adjacent to the AWAKEN wind farms. In the present paper, we define a new idealized stability condition with the AWAKEN lidar and meteorological tower measurements taken during the experimental campaign. The method and filtering process behind the selection of the neutral conditions are described in the Simulation Setup section of the paper.

Previous comparisons between simulations under stable and unstable conditions for the King Plains wind farm using a variety of codes (AMR-Wind, FLORIS, Nalu-Wind, WRF-LES) showed large discrepancies among the results. The largest differences were in the downstream wake profiles as shown in Figs. 1 and 2 for the unstable and stable cases, respectively.¹² The models are shown to differ in respect to both wake recovery distance as well as the maximum velocity deficit especially for both WRF-LES and FLORIS in the unstable case.

One primary motivation for the present study is to advance the understanding of the wakes of modern utility-scale turbines using a different suite of models such as WRF-LES and AMR-Wind, through a more comprehensive investigation of the discrepancies in Ref. 12. The large differences in the modeled predictions in Ref. 12 are also problematic as upcoming simulation benchmark activities will also compare models in reference to the AWAKEN experimental data. The format of these future benchmarks will be similar to the study of Doubrawa et al.,¹⁸ which compared models across a wide range of fidelities and assessed the wakes of a single turbine during neutral, unstable, and stable atmospheric conditions.

The results of Cheung et al.¹² demonstrated the need for identifying sources of model uncertainty and error prior to the benchmark studies as significant discrepancies between modeled results are generally difficult to investigate due to the differences between simulation methods. A flow field can be modeled through LES, Reynolds-Averaged Navier–Stokes (RANS) models, and lower-order wind and wake models. Turbine models also have varying levels of fidelity ranging from actuator line and fully resolved blade-boundary-layer schemes to actuator disk and blade element momentum theory (BEMT). A detailed summary of the different numerical methodologies and simulation setups in the present study using AMR-Wind, Nalu-Wind, WRF-LES, and FLORIS is given in the subsequent section. In addition to the separate ABL models, the individual wind turbine models used in each code are also discussed. These wind field and wind turbine models are the same as used in Cheung et al.,¹² but with different simulation methodologies aimed at reducing potential sources of model differences.

In Ref. 12, it was speculated that one reason for the large discrepancies was the unstable flow condition, as large-scale turbulent structures in the unstable atmospheric boundary layer (ABL) can cause sizable local differences in the flow field. Even if the domain-averaged parameters are similar between simulation runs, the turbines can potentially experience significantly different incoming wind. Also, the results compared in Ref. 12 were averaged over a large number of closely packed wind turbines within the King Plains wind farm. The complexity of multi-turbine-wake interactions, by itself a potentially large source of model discrepancy, is another factor to consider.

To remove the model uncertainty accompanying the unstable ABL condition and multi-turbine–wake interactions, the present study employed smaller-scale simulations using a single turbine in a reduced domain with a neutral stability condition. In addition, specific care was taken to improve simulation methods and enforce similar conditions for the domain, atmospheric inflow, and sampling across all model runs to minimize model discrepancies. It should be noted that the purpose of utilizing field data for the neutral ABL condition in the present study was not to create a precise validation case for different wind farm models; rather, it was to ensure a level of realism for the simulated scenarios. The “true” neutral ABL in the field is generally in a state of nonequilibrium due to the short nature of the day-to-night transition. Hence, the statistics garnered from experiment will be more variable and less certain than the quasi-steady-state condition of simulations results.

The literature on idealized LES of neutral boundary layers also reveals many difficulties associated with the region near the surface and the consistency with Monin–Obukhov similarity theory (MOST).^10,31 Often, simulations are evaluated on the basis of the non-dimensional vertical gradient of mean streamwise velocity, $Φ=κzu*∂U∂z$⁠, where $κ$ is the Von Kármán constant, $z$ is the height above ground level, $u*$ the friction velocity, and $∂U∂z$ the vertical gradient of mean streamwise velocity. A consistent treatment of the near-wall region should produce a constant $Φ$ in the surface layer. In Andren et al.,³ four LES codes were used to simulate a neutral boundary layer with a geostrophic wind speed of 10 m/s. It was noted that the models displayed an overshoot in $Φ$ near the surface, which deviated from the MOST prediction that the value should be constant in the surface layer. The subgrid-scale (SGS) model with backscatter produced a better agreement with MOST. Next, Kosović²³ developed a nonlinear SGS model that represents both the backscatter and anisotropy (NBA) and was successful in mitigating the overshoot. Its subsequent implementation in the WRF model was able to mitigate, but not completely remove, the overshoot in comparison with linear SGS models when using aspect ratios in the range between 3 and 4.²⁸ In Brasseur and Wei,⁸ a theoretical framework was developed to address the physical and numerical aspects of this overshoot and its relationship with the SGS model constants, grid resolution, and aspect ratio. This framework was designed for and successful in making the LES prediction for Law-of-the-Wall (LOTW)-normalized gradient constant over the surface layer. In a more recent study, intercomparison for a single wind turbine wake that included neutral conditions, the authors pointed out that obtaining consistent inflow conditions across models was one of the main challenges, with pronounced differences in the TKE profiles.¹⁸ Despite its apparent simplicity, the LES of the neutral boundary layer over flat terrain can be challenging from the modeling comparison point of view.

It must also be noted that a complete discussion regarding the impact of numerical dissipation from the LES solvers is outside of the scope of the current work and should be comprehensively addressed in a separate study. For the sake of comparison, the current study seeks to compare simulations with closely aligned mesh and model settings to achieve similar ABL characteristics. Numerical dissipation and dispersion, along with the mesh sizes employed, characterizes the effective resolution of the LES, which determines which turbulence features in the ABL and wakes that can be captured. For example, the overshoot of the normalized mean shear in the surface layer of the simulations in the current study is a known phenomenon but has not yet been quantified. Jayaraman and Brasseur²¹ demonstrated the dependence of grid resolution and aspect ratio for this overshoot. Due to innate code differences and discrepancies in the applied grid resolutions between AMR-Wind, Nalu-Wind, and WRF-LES, there are expected deviations in numerical dissipation in the LES results. In future studies, we hope to thoroughly document the effects of grid resolution for these codes on turbulence features of the ABL and wake.

The inflow conditions studied in this article were selected from measurements taken at the AWAKEN site and available from the Data Archive Portal (DAP), also known as the Wind Data Hub (WDH).¹⁵ Site A1 is located 300 m from the instrumented turbine of interest (i.e., H05 in the King Plains wind farm) that was selected for study, and the instruments at A1 used for characterization of the inflow velocity profile were a ground-based WindCube v2 profiling lidar (DAP code: sa1.lidar.z02) and a ground-based Halo XR scanning lidar (sa1.lidar.z03.c1). Measurements from these instruments were complemented by a nacelle-mounted, forward-facing WindCube (rt1.lidar.z01.00) on turbine H05. A surface meteorology station located at A1 (sa1.sonic.z01.c0) provided information on the stability state of the atmosphere.

The targeted atmospheric flow condition for this study was a neutral boundary layer as defined by Krishnamurthy et al.,²⁴ with wind speed in Region II where there is maximized rotor efficiency. The 10 min window selected for study began on September 20, 2023, at 18:30:00 local time [i.e., Greenwich Mean Time-5 (GMT)]. This window, which occurred 50–60 min before sunset during fair conditions, represents an unstable-to-stable atmospheric transition as suggested by the local stability parameter, z/L, from surface meteorological data in Fig. 3. No measure of the ABL height was available for this time period, and thus a measurement of the global stability, z_i/L, is not available. However, feasible values for the unstable and stable ABL heights at the AWAKEN site are 2000 and 500 m, respectively,²³ from which estimates of z_i/L at 17:30 and 19:30 local time are −38 and 13, respectively. This change in z_i/L was a large shift in the global ABL over a short, 2 h period so that the ABL was likely out of equilibrium during the entire event. This state of nonequilibrium should be borne in mind when comparing the data in this validation study where the simulated conditions are in equilibrium for validation purpose.

The mean wind speed and wind direction profiles for the selected 10 min period are shown in Fig. 4. Over the height of the rotor, the wind speed profile follows an approximately power-law relationship according to the profiling lidar data, which was deemed more reliable than that of the scanning lidar for the lower altitudes around rotor height as the scanning lidar's coverage did not span the full rotor. About 70 m above rotor height, the wind speed profile peak suggests the presence of a low-level jet. The wind direction profile is approximately linear with only mild veering.

The target atmospheric conditions derived from 10 min measurements are given in Table I. The wind velocity data all derive from the profiling lidar except for the value of turbulence intensity (TI), which stems from the nacelle-mounted lidar data at a range of 400 m ahead of the turbine. This latter measurement, which uses the WindCube definition of TI,²⁵ was considered preferable to the turbulence intensity profile from the profiling lidar because of unexpectedly low turbulence readings from the profiling lidar. In figures below, the TI profile shown is that of the profiling lidar with a uniform scaling of 5.4 applied so that the hub-height TI matches that from the nacelle-mounted lidar, while data from the scanning lidar is not shown due to poor data availability at heights of interest.

The turbulence level represents the greatest source of uncertainty in the observed atmospheric condition because of drastically different values of TI measured by the three lidar instruments, and these differences between lidars might be explained by the lidar-related phenomena of underestimation of TI due to probe-volume averaging and overestimation due to cross-contamination as explained in Sathe et al.,⁴⁰ as well as by the nonequilibrium state of the lower atmospheric during transition from unstable to stable that saw the TI level drop by ∼0.4% over these 10 min of interest according to the nacelle-mounted lidar.

In the simulation precursor development to follow, we ignore the possible low-level jet profile and target an idealized power-law relationship with the 0.116 exponent as calculated with a reference height of 90 m and a least squares fit to the wind speed from the profiling lidar over the rotor height. The veer target is 3.28° over the rotor height as calculated from a linear least squares fit of the wind direction from the profiling lidar over the rotor height. Wind direction data from the profiling lidar were also combined with the turbine yaw heading for this window to calculate a mean yaw misalignment of 4.33°.

The LES code bases AMR-Wind^41,43 and Nalu-Wind,^17,43 both part of the ExaWind code suite, were used to perform high-fidelity simulations of the ABL. AMR-Wind solves the incompressible Navier–Stokes equations using an approximate projection method used in Incompressible Adaptive Mesh Refinement (IAMR)² and incflo.⁴⁴ It is a semi-staggered scheme with a semi-implicit Crank–Nicolson approach. An implicit scheme is used for the viscous terms and the advection term is handled explicitly using an upwind finite-volume method using the WENO-Z (Z-Type Weighted Essentially Non-Oscillatory) algorithm. A single-equation turbulent kinetic energy model is employed for the subgrid scale stresses following the formulation of Moeng.³² 

The multiphysics, massively parallel code, Nalu-Wind, performs LES of the ABL using a node-centered finite volume discretization on an unstructured grid to solve the acoustically incompressible Navier–Stokes equations. An implicit second-order backward difference scheme is used for time integration. Similar to AMR-Wind, a one-equation, constant coefficient, TKE model is used to model the subgrid-scale stresses.⁴⁶ 

Both the AMR-Wind and Nalu-Wind codes have similar approaches to modeling the ABL. An ABL wall function based on the Monin–Obukhov theory is used for the bottom boundary, where a prescribed surface roughness length is used to alter the surface layer profile. An inflow/outflow boundary condition is used at the top boundary.⁴⁵ The inflow/outflow condition uses a potential flow solution in wavenumber space, which dampens high-frequency disturbances at the upper boundary. The boundary condition on the domain sides is periodic. In place of specifying a mean horizontal pressure gradient, ABL forcing source terms are provided to the momentum equation to drive the flow to a predetermined velocity at a specific height.

Both AMR-Wind and Nalu-Wind include a Boussinesq buoyancy model and Coriolis source term to capture the effects of atmospheric stratification and wind veer, the Coriolis term using a specified latitude of 36.572° to match the approximate location of the AWAKEN wind farms. Both codes solve the potential temperature form of the acoustically incompressible Navier–Stokes equations and utilize a potential temperature profile in the initial conditions, where the capping inversion height was set at 650 m for the present study. At the lower boundary condition, surface temperature or surface heat flux may be prescribed and for the present simulations, the surface heat flux was fixed at 0 K/m for both AMR-Wind and Nalu-Wind. Table II displays the associated mesh sizes and computational costs of AMR-Wind and Nalu-Wind compared to the third LES code, WRF-LES. AMR-Wind is shown to be the most scalable and computationally efficient of the three LES codes.

Both AMR-Wind and Nalu-Wind have similar processes for setting up the neutral boundary layer case. This is the same method that is described in more detail in Ref. 20. A precursor using periodic boundary conditions with a domain size of 5120  $×$ 5120  $×$ 960 m³ tall with an isotropic grid size of 10 m in each Cartesian direction was run for 20 000 s to establish a flow with a forced hub-height wind speed of 6.27 m/s, air density of 1.09 kg/m³, and a near-ground temperature of 305 K. The prescribed surface roughness for AMR-Wind and Nalu-Wind was 0.03 m. After the flow was established in the first 20 000 s, the simulations were run for another 5000 s to provide a sufficient search window to find the 10 min statistics that best matched the average hub-height velocity, shear, veer, and TI listed in Table I. East–west corresponded to the x-direction and north–south to the y-direction in the simulation domains. To account for the mean yaw misalignment of 4.33°, the wind direction was offset that amount to the south from being purely in the x direction while keeping the turbine aligned with the grid.

The 10-min window that was chosen for each code was then rerun on a domain with refinement boxes that decreased the grid size from 10 to 5 m, and then to 2.5 m closest to the turbine as shown in Fig. 5. The mesh refinements also retain the cubical nature of the grid cells. The inflow planes of the 10-min windows were used as inflow boundary conditions, and the flow for the entire domain at the beginning of the period was used as the initial condition. The OpenFAST model used for the NREL 2.8-127 was placed at x = 2000 m and y = 2560 m to be slightly upstream of center to allow more room for wake flow. A case using the same mesh was also run without the turbine so that differences in the flow would be purely caused by the turbine and not by mesh resolution.

To compare the results from the different codes, a set of sample planes were agreed upon to be used in each. These included a horizontal x–y slice of the entire domain at the hub-height (z = 90 m), the top tips of the blades (z = 153.5 m), and the bottom tips of the blades (z = 26.5 m). Smaller x–y slices that were 30D  $×$ 8D (where D is 127 m, the diameter of the rotor) were at heights of 90 m (hub height) and 280.5 m. These were centered around the turbine in the y direction and started 10D upstream of the turbine. Several vertical slices in the y–z direction were used to capture flow before and after the turbine at distances of 10D upstream to 20D downstream at intervals of 2D. These extended 4D to either side of the rotor center and were 280 m (2.2D) tall. There were also vertical slices in the x–z direction that were 30D in the x direction and were also 280 m tall. These were placed at 0D as well as 4D and 2D on either side of the turbine. The placement of these planes can be seen in Fig. 6.

The simulations used the Advanced-Research WRF model version 4.4,⁴² which solves the compressible Euler equations for the three spatial dimensions and time. The time integration employs a third-order Runge–Kutta scheme with a smaller time step for acoustic effects. The spatial discretization of advective terms employs a fifth-order scheme for the horizontal direction and a third-order scheme for the vertical direction. An Arakawa C-grid staggering is used horizontally, whereas the vertical direction employs a hydrostatic pressure-based coordinate.

The WRF model allows idealized simulations by neglecting physical processes related to cloud microphysics, radiation, moisture, and surface heterogeneity. These are standard simplifications for studies focused on canonical boundary layers over flat terrain.^{4,28–30,37,39} WRF is a multiscale framework that models turbulence differently depending on the grid resolution. In the mesoscale range, where the grid spacing (dx) is greater than 1 km, turbulence is fully modeled by planetary boundary layer schemes. In the microscale range (dx < 100 m), the LES resolves turbulent motions larger than the spatial filter (grid) size, whereas the effect of the interaction between resolved and unresolved motions in the resolved-scale velocity field is represented by a subgrid-scale model.¹⁶ Our approach involves microscale LES using the nonlinear backscatter and anisotropy (NBA)¹⁸ subgrid-scale model implemented in WRF by Mirocha et al.²⁸ 

A setup consisting of two computational domains one-way coupled by nesting is used (Fig. 7) for the WRF-LES simulations. The parent and the nest domains have dimensions of 5000  $×$ 5100 m and 3500  $×$ 2000 m, respectively, with a model top at 1000 m. The main role of the parent domain is to produce average boundary layer characteristics that are temporally steady, which are then fed into the nest via one-way coupling. Both domains use a constant horizontal grid resolution of 10 m and a vertical grid with a constant resolution of 2.54 m in the first 300 m above ground and stretched out up to 34 m at the model top. The finer grid in the vertical direction with an aspect ratio in the range between 3 and 4 was shown beneficial for WRF-LES simulations.³⁰ The maximum vertical grid size in the first 500 m above ground is limited to 10 m. The single wind turbine, located at x = 500 m and y = 1000 m in the nest domain, is represented in the simulations using the generalized actuator disk (GAD) that was originally proposed and implemented in the WRF model by Mirocha et al.²⁹ and later adapted in Refs. 1 and 4. The GAD method uses turbine control and blade-section aerodynamic properties tables to compute lift and drag forces, which are then projected onto the rotor-swept area. The turbine faces west and the yaw is fixed.

The simulations initialize with nonturbulent, user-specified profiles for wind speed components, potential temperature, and water vapor. Wind direction components were assumed to be constant with height. The potential temperature profile was the same as prescribed for AMR-Wind and Nalu-Wind, with the capping inversion height set at 650 m. The simulations considered dry air, and the water vapor was set to zero at all heights.

The top boundary condition imposes zero vertical velocity and free slip for the horizontal velocity. The upper 300 m of the domain contains a Rayleigh damping layer (damping coefficient of 0.003 s⁻¹) to avoid the downward reflection of gravity waves. A modified surface layer model³ computes the drag force at the first grid cell above the surface using the local horizontal wind components and a drag coefficient based on universal functions for the MOST.³¹ This lower boundary condition imposes a shear stress compliant with the MOST assuming a flat terrain with a homogeneous roughness length of 0.3 m. A constant sensible heat flux of 0 W/m² is also imposed. Periodic lateral boundary conditions are used for the parent domain to allow the atmospheric boundary layer winds and turbulence to develop from the nonturbulent initial condition.

In addition to the adjustment of the initial condition profiles owing to friction and turbulence, the Coriolis effect also produces an inertial oscillation that causes temporal variations in the boundary layer characteristics. Therefore, the approach adopted to produce the inflow considers enough time to fully develop turbulence and for inertial oscillation effects in the mean wind to be small such as in Ref. 37. Considering these constraints, the initial conditions' values were iteratively fine-tuned to reach the desired hub-height wind speed, direction, shear, veer, and TI.

Horizontal plane averages of selected variables are shown in Fig. 8. The wind speed (⁠ $WS90$⁠) and direction (⁠ $WD90$⁠) at 90 m height are computed based on the plane-averaged horizontal wind velocity components. The wind shear exponent across the rotor (⁠ $αr$⁠) is obtained with a nonlinear least squares curve fit of the wind speed profile to a power law profile. The veer across the rotor is the subtraction between the wind directions at the heights of the top and bottom tips (⁠ $ψr=WD153.5−WD26.5$⁠). The vertical velocity variance (⁠ $w′w′90_$⁠) is computed based on the vertical velocity averaged across the horizontal plane (⁠ $w_90$⁠) and the local departures from the average (⁠ $w90′=w90−w90_$⁠) at 90 m height. These computations were consistent across the different simulations and field data.

Before the nest domain is activated, to reduce the computational cost of the parent domain's spin-up phase, a smaller domain (2500  $×$ 1700 m) with periodic lateral boundary conditions is run for at least 8 h to fine-tune the initial conditions to match the target wind profile characteristics. Three iterations were necessary for the initial conditions for wind speed and direction, which were forcing 1 (7.5 m/s and 295°), forcing 2 (8.5 m/s and 287°), and forcing 3 (8.5 m/s and 291°). Once this is achieved, the tiling approach³¹ reconstructs the initial condition of the original parent domain (5000  $×$ 5100 m) using 2  $×$ 3 tiles with the end state of the smaller domain (2500  $×$ 1700 m). Starting at 05:00, the parent domain is run for 4 h to break the periodicity of the turbulence. The differences in the plane-averaged variables between the small and large domain with forcing 3 are very small and reveal the effectiveness of the tiling approach. The wind speed and direction vary slowly after 07:00 and reach the targets between 07:42 and 07:52, whereas the shear stabilizes much earlier and slightly overestimates the target. The wind veer stabilizes slightly below the 3.28° target. The nest domain was activated after 07:30 for 30 min from which the analysis refers to the period between 07:42 and 07:52 (Fig. 8, gray shaded area). Because the parent and nest domains have the same spatial resolution, no turbulence spin-up technique is required in the nest domain.

The FLORIS model is a low-computational-cost, python-based, engineering tool used to simulate steady-state wind turbine and wake behavior.¹⁹ FLORIS contains an array of turbulence, wake, and user-defined models that are calibrated to match the characteristics of the relevant turbine model and atmospheric conditions. In FLORIS, the velocity fields are represented by analytical solutions to the simplified Navier–Stokes equations. The wake itself is composed of four separate mathematical models representing the wake deficit, wake deflection, wake turbulence, and wake combination. The turbine models are defined using the design power and thrust curves of the relevant turbine, and are represented by actuator disks in the simulation.

A number of wake turbulence, deficit, combination, and deflection models are available in FLORIS. In Ref. 12, the selected wake turbulence and deficit models used were the TurbOPark model and Gauss–Curl Hybrid (GCH) model, respectively. However, these FLORIS model selections were specifically made for large cluster arrays of turbines in the previous study configuration. For the one-turbine simulations in the present study, the primary wake deficit model instead was varied between the cumulative curl (CC) and the empirical Gaussian (EMGAUSS) models. The sum of squares freestream superposition (SOSFS) was applied as the wake combination model, a Gaussian wake deflection model was used, and the wake turbulence model chosen was Crespo–Hernandez.¹⁴ 

As a wind farm simulation software designed primarily for optimization and controls-based studies, FLORIS has a significant number of limitations compared to the aforementioned three LES codes. Temperature profiles such as at the surface or capping inversion are not considered in FLORIS, along with more complex atmospheric phenomena such as the impact of Coriolis forcing on veer and estimation of geostrophic wind. Hence, FLORIS is unable to truly model a “neutral” stability state, merely the specified hub-height wind velocity, wind direction, rotor shear, and other inputs. However, in contrast to the high computational cost requirements shown for the LES codes in Table II, FLORIS is very low cost and efficient, completing simulations on the order of minutes using a single computer processor.

Using FLORIS, a single NREL 2.8-127 turbine was simulated using the defined inflow wind conditions listed in Table I. The turbine was placed in the center of the domain. Due to the specified mean yaw misalignment of 4.33°, the generated FLORIS mesh for the simulation was angled. So for post-processing, the results were interpolated onto a 2000  $×$ 1000 m² rectangular grid in x–y with 4 m spacing. For normalization of hub-height statistics, the input hub-height velocity of 6.27 m/s was used. Sampling plane data were extracted at locations corresponding to the LES results. Results using both the cumulative curl and empirical Gaussian wake deficit models are shown due to the difference in wake structures between the FLORIS simulations.

The turbine model chosen for the neutral simulations was the NREL version of the GE 2.8-127 as this model was used in the previous study of Cheung et al.¹² The NREL 2.8-127 OpenFAST³⁶ model was developed by scaling publicly available reference turbines to match the correct hub height, rotor size, and power rating of the actual GE turbines at the King Plains wind farm site.³⁸ No proprietary wind farm or turbine data were used. Similarly, actuator disk models (ADM) coupled to the OpenFAST turbine solver are used to represent the turbine dynamics for the LES cases using AMR-Wind and Nalu-Wind. Note that the Nalu-Wind implementation used an actuator disk with a filtered lifting-line correction (FLLC) model²⁶ and included only the streamwise force component on the actuator disk, while the AMR-Wind implementation did not use the FLLC model, but included both the streamwise and azimuthal forces on the actuator disk. For the WRF model, the wind turbine is represented as a generalized actuator disk (GAD)²⁹ using the NREL 2.8-127 model. FLORIS also represents turbines using actuator disks and its turbine model was specified using the design power and thrust curves of the NREL 2.8-127 model, which were generated using uniform, steady OpenFAST simulations.

To simulate the hub and nacelle, the three LES codes applied a Gaussian drag body force model,²⁷ which used a specified drag coefficient and a reference area corresponding to the hub area of the NREL 2.8-127 model. FLORIS does not have a nacelle or hub model.

In this section, modeled inflow results along with visualizations of the background flow are assessed. Figure 9 shows a comparison of the selected 10 min averaged precursors from the LES codes (AMR-Wind, Nalu-Wind, WRF-LES) against the profiling lidar data, which has been scaled in the case of the TI data, from DAP. As the FLORIS results did not utilize a precursor ABL simulation, FLORIS profiles are not included in this section. The modeled profiles were taken at 2.0–2.5D upstream of the turbine. A summary of relevant rotor and hub-height values from the field data and modeled inflow results is also shown in Table III.

Figure 9 shows good agreement for the horizontal wind velocity profiles between the profiling lidar data with the AMR-Wind, Nalu-Wind, and WRF-LES simulations. The WRF-LES wind profile has slower winds in the bottom region of the rotor, which produces a shear exponent that is slightly higher than that derived from the observations (Table III). For Nalu-Wind and AMR-Wind, the surface roughness was varied to match TI, veer, and the shear exponent. However, for WRF-LES, these parameters were only weakly sensitive to roughness length, preventing much user control. At hub height, veer is similar between all three models and the field data but none of the LES models are able to capture the larger variations in wind direction observed at higher altitudes in the lidar data, indicating a limitation of the present LES models. Both AMR-Wind and Nalu-Wind have similar TI profiles with slightly higher values than the field data. The WRF-LES profile has a lower rotor-averaged TI than the observations, and higher TI below the rotor in comparison with the other LES codes.

The neutral ABL in all simulation methods had an inversion layer specified between z = 600 m and z = 750 m, and used a constant potential temperature profile of $θ$ = 305 K below the inversion layer height (see Fig. 10). This ensured that the buoyancy effects would be negligible near the lower surface and limited the growth of the boundary layer consistently between AMR-Wind, Nalu-Wind, WRF-LES.

All three simulation methods also used Monin–Obukhov similarity theory as a part of the wall model boundary condition applied to the lower surface. The resulting values for the computed friction velocity, $uτ$⁠, over a 10 min duration of the wind turbine simulation are very similar between AMR-Wind and Nalu-Wind (Fig. 11). The observed averaged value of $uτ$ = 0.34 m/s for WRF-LES is significantly higher than the horizontally averaged $uτ$ = 0.29 m/s for both AMR-Wind and Nalu-Wind. A potential reason for this discrepancy is that both AMR-Wind and Nalu-Wind use a method from Basu et al.⁵ to compute the friction velocity from the horizontally averaged velocity at a specified height from the lower surface. In contrast, WRF-LES computes the friction velocity from the local surface shear stress at each grid cell, which is then averaged in space. This difference in friction velocity calculation also results in a more fluctuating friction velocity profile for WRF-LES. The friction velocity time series for AMR-Wind and Nalu-Wind are variable over time, but on a much smaller scale.

Figure 14 compares the 10 min averaged velocity contours among the simulations. The FLORIS EMGAUSS model has a higher wake velocity deficit than the other models, while the three LES codes have similar wake shapes and deficits. The streak of high velocity in the middle of the near-wake region in the three LES models is due to the nacelle model, which allows some flow to pass through the rotor center while FLORIS does not. The FLORIS CC model is shown to have a more similar wake structure compared to the LES than the FLORIS EMGAUSS model.

Figure 15 compares instantaneous velocity flow fields between the three LES models. The FLORIS model is steady state and is not included. The WRF-LES results show some significant differences in the flow field compared to the others. The coarser resolution of the WRF model (10 m compared to 2.5 m in Nalu-Wind and AMR-Wind) causes the smaller-scale structures present in both the Nalu-Wind and AMR-Wind results to be smeared out. Also, the wake persists in the Nalu-Wind and AMR-Wind contours for up to 20D. For WRF-LES, the wake breaks up into smaller structures around 10D.

The simulated turbine quantities were also assessed among the four models and the design curve values of the NREL 2.8-127 turbine, which is used only as a reference for comparison. Table IV summarizes the power, thrust, and torque coefficients along with tip-speed ratio (TSR) among the simulations. The FLORIS CC and EMGAUSS runs produced the same turbine outputs and are summarized as FLORIS below. Neither FLORIS nor WRF-LES used OpenFAST and their standard outputs did not include any turbine quantities besides generator power.

In Table IV, AMR-Wind and WRF-LES are shown to have similar predictions to the turbine design power coefficient, while Nalu-Wind and FLORIS under-predict the power coefficient by ∼10%. In Fig. 17, the generator power time series of the four models is shown and the WRF-LES distribution is shown to have significantly higher-frequency content in its power profile compared to AMR-Wind and Nalu-Wind. This is surmised to be a result of AMR-Wind and Nalu-Wind's coupling to OpenFAST which has complex and integrated aero-servo-elastic models that act as a filter for the temporal variations in turbine power, whereas the simpler GAD model in WRF-LES is significantly more reactive to the incoming wind and hence shows a higher frequency content in the power distribution. AMR-Wind is also shown to have higher predictions for normalized thrust than Nalu-Wind in Table IV, although both are slight over-predictions compared to the turbine design thrust coefficient.

Figure 18 shows the normalized centerline velocity. Due to differences in the wake solutions among the codes, the data were taken at respective hub-height lines taken from the turbine to the location of maximum wake velocity deficit at 10D distance downstream of the turbine. As this case considers the neutral ABL, we assume the maximum velocity deficit is located close to the hub-height plane due to the lack of updraft or downdraft convective structures. This post-processing centerline calculation method differs from Ref. 12, which considered both the stable and unstable ABL, and so the results in Figs. 1 and 2 may not include the true maximum velocity deficit due to the presence of updraft and downdraft structures.

There is an observable spike in velocity immediately downstream of the turbine that is caused by wind acceleration around the nacelle. The spike is higher for WRF-LES in comparison with the other codes, which is possibly related to the differences in spatial resolution. Neither FLORIS model results capture the impact of upstream turbine induction like the three LES models, which agree well in the induction zone. The lack of a nacelle model in FLORIS results in the over-prediction of the wake deficit immediately downstream of the turbine. The cumulative curl model in FLORIS captures the profile shape of the wake much better than the empirical Gaussian model with the latter predicting maximum velocity deficit at the location of the turbine.

Overall, there is good agreement between AMR-Wind, Nalu-Wind, WRF-LES, and FLORIS CC for the downstream wake in contrast to the corresponding plots in Cheung et al.¹² shown in Figs. 1 and 2. Similar to the FLORIS results in Ref. 12, the FLORIS EMGAUSS profile predicts a higher wake velocity deficit compared with the LES results in the near wake at $x/D$ < 2. The FLORIS CC, Nalu-Wind, and WRF-LES profiles have similar maximum velocity deficits and wake shapes. The models tend to agree in the far-wake region (⁠ $x/D$ > 5), but the wake recovery was faster for WRF-LES. Also, WRF-LES and AMR-Wind have similar profiles in the near-wake region (1 < $x/D$ < 4).

Potential causes for the better agreement among models relative to the previous work¹² are worth discussing especially regarding the most pronounced discrepancy in Fig. 1 for the much shorter wake for WRF-LES in convective conditions. In Ref. 12, the convective WRF-LES employed a coarse parent domain with periodic boundary conditions that fed the nest domain. To accelerate turbulence spin-up in the nest domain, the Cell Perturbation Method (CPM)³⁵ was employed. The approach adopted in the present study employs parent and nest domains of the same grid resolution and the CPM is unnecessary in the nest. The latter approach is similar to that used by AMR-Wind and Nalu-Wind, which would explain, in part, the better agreement in contrast with the previous work.¹² The much faster wake recovery in the previous study could also be related to high levels of TKE because of the relatively short spin-up time used for the nest and the somewhat stronger winds (closer to 10 m/s) in comparison with other models (close to 9 m/s). Furthermore, with a hub height wind speed in the range of 9.5–10 m/s, the power curve is close to the knee region after which the thrust force decreases.

Figure 19 shows comparisons of both normalized hub-height velocity and wake-added TKE profiles among the models at different downstream distances from the turbine: 1D, 2D, 4D, 6D, and 8D. Due to being a steady state model, no FLORIS results are included in the wake-added TKE profiles. Wake-added TKE is the difference between the average TKE of the turbine simulation vs the average TKE of the precursor simulation. This quantity represents the turbulence generated from the tip vortex decay and shear inside the wake profiles and this can be masked in the overall TKE profiles.

In the near-wake at 1D and 2D, both FLORIS model profiles show a different Gaussian-like wake shape compared to the LES results, which is due to the lack of a nacelle model in FLORIS and an incorrect assumption of the near-wake behavior. However, the FLORIS velocity profiles do converge onto the LES profiles at farther downstream distances. AMR-Wind predicts higher wake-added TKE in the near-wake region compared to Nalu-Wind and WRF-LES, which potentially is a result of AMR-Wind modeling the swirl component of the flow while Nalu-Wind does not. There is good agreement between all three LES codes for the wake velocity and wake-added TKE profiles in the far-wake (>5D) region.

Figure 20 shows comparisons of the y–z velocity contours among the LES models at different downstream distances from the turbine: 0D, 2D, 4D, 8D and 10D. AMR-Wind and Nalu-Wind show higher wake velocity deficits than WRF-LES, which is most prominent at x/D = 2. Both displacement and momentum areas, along with the y–z contours, capture the three-dimensional distribution of the wake as opposed to a single hub-height plane. In both the near-wake and far-wake regions, AMR-Wind and Nalu-Wind consistently display a higher vertical and wider horizontal wake extent than WRF-LES. This is especially apparent in the 10D contours where the WRF-LES wake has nearly dissipated into the surrounding flow while the wake structure is still present in Nalu-Wind and AMR-Wind. These results correspond with the profiles in Fig. 16, where AMR-Wind and Nalu-Wind had the more displacement and momentum area than WRF-LES, which indicates more momentum is being extracted from the flow for AMR-Wind and Nalu-Wind.

A previous model comparison study in Cheung et al.¹² using large-scale simulations of the AWAKEN site under stable and convective stability states demonstrated large discrepancies among the LES models of AMR-Wind, Nalu-Wind, WRF-LES, and the engineering model FLORIS. These discrepancies were most pronounced in the modeled turbine wakes and the need to have consistent simulation setups between code bases for future AWAKEN benchmark studies prompted the present study.

The present study aimed to align these simulation setups compared to Ref. 12 to see whether the model discrepancies could be reduced or more easily understood. The new study used the same four codes as the previous study and similarly modeled the AWAKEN site characteristics including its turbines and atmospheric conditions. As there were concerns about the large-scale convective structures in the previous study causing local inflow variations among the models, a new neutral stability condition was defined using AWAKEN lidar and meteorological instrument data obtained from the DAP to avoid updraft and downdraft effects. In addition, a smaller-scale one-turbine simulation was chosen for the simulation setup to eliminate the impact of upstream turbine waking and ensure consistency among the modeled inflows.

Some improvements were also made to the simulation methodologies. In FLORIS, it was found that the applied Gauss–Curl Hybrid wake deficit model in Ref. 12 was not ideal for capturing the downstream turbine wake structure and a new cumulative curl wake deficit model was found to better match LES results. In WRF-LES, a new mesh structure with greater resolution was employed compared to Ref. 12, and the new wake results showed much better correspondence with the other LES codes. During the simulation setups, strong emphasis was placed on matching the modeled inflows with the experimental data for a wide range of atmospheric parameters including wind speed, TI, shear, and veer.

These efforts were successful as greater similarities were observed in the modeled inflow, turbine, and wake results in the present study compared to Ref. 12 and the far-wake profiles showed strong agreement among all simulations. The results indicate some combination of picking a more appropriate FLORIS model, increasing the WRF mesh resolution, and being more cognizant of matching inflow conditions (e.g., nonwaked neutral inflow) among models significantly improved model agreement of the far-wake profiles. This study can potentially inform future AWAKEN benchmark efforts on best practices.

It should be noted that these models, especially complex LES codes, are designed with different purpose in mind and it will be difficult to replicate the exact setup between runs. For example, WRF-LES includes some atmospheric physics not present in AMR-Wind and Nalu-Wind, while the latter codes have more detailed turbine representations. FLORIS is much more computationally inexpensive but has significant modeling limitations such as in the near-wake. In the present study, there were still noticeable inflow differences observed among the models such as the lower shear exponent and resolved TKE in the WRF-LES results. These potentially account for discrepancies in the turbine and wake predictions, such as the relatively high power predictions and momentum/displacement areas from WRF-LES. A deeper examination of the wake statistics also showed that, despite the hub-height wake profiles matching among the LES codes, there were significant differences in the downstream wake structures that could be observed in y–z velocity contours. Hence, it is important to be aware of these sensitivities while performing model benchmarking studies and to investigate a wider range of wake metrics.

The authors would like the acknowledge the assistance of others in the AWAKEN simulation group who contributed to previous simulation work including researchers from NREL (Nicholas Hamilton and Ryan Scott) and PNNL (Raj Rai and Colleen Kaul), as well researchers from NREL (Stefano Letizia) and LLNL (Sonia Wharton) who provided insight on instrument characteristics.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA0003525. The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This work was accomplished through funding from the U.S. Department of Energy Wind Energy Technologies Office.

This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725, which was provided through the ASCR Leadership Computing Challenge (ALCC) program.

W. Radünz would like to acknowledge the São Paulo Research Foundation (FAPESP), Grant Nos. 2022/04474-6 and 2023/12599-6, for financial support, and the National Renewable Energy Laboratory (NREL) for the Eagle supercomputer resources, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Also, to acknowledge Miguel Sanchez Gomez for discussions on the tiling approach, Branko Kosovic for discussions on subgrid-scale modeling, and Julie Lundquist for hosting him at the University of Colorado Boulder under the FAPESP Grant No. 2023/12599-6. B. S. Carmo gratefully acknowledges partial support of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brazil (CAPES)—Finance Code 001, as well as support of the RCGI—Research Centre for Gas Innovation, hosted by the University of São Paulo (USP) and sponsored by FAPESP–São Paulo Research Foundation (2020/15230-5) and TotalEnergies, and the strategic importance of the support given by ANP (Brazil's National Oil, Natural Gas and Biofuels Agency) through the R&D levy regulation. B. S. Carmo thanks the Brazilian National Council for Scientific and Technological Development (CNPq) for financial support in the form of a productivity Grant No. 314221/2021-2.

The authors have no conflicts to disclose.

Alan S. Hsieh: Conceptualization (lead); Data curation (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Lawrence C. Cheung: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal). Myra L. Blaylock: Data curation (equal); Methodology (equal); Writing – review & editing (equal). Kenneth A. Brown: Data curation (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Daniel R. Houck: Data curation (equal); Formal analysis (equal); Writing – review & editing (equal). Thomas G. Herges: Conceptualization (equal); Funding acquisition (equal); Project administration (equal). Nathaniel B. deVelder: Resources (equal); Software (equal). David C. Maniaci: Conceptualization (equal); Formal analysis (equal). Gopal R. Yalla: Formal analysis (equal); Software (equal). Philip J. Sakievich: Resources (equal); Software (equal). William C. Radünz: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal). Bruno S. Carmo: Formal analysis (equal); Investigation (equal).

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