A multi-fidelity model intercomparison for wake steering of a large turbine in a conventionally neutral atmospheric boundary layer

Wake steering is a promising control strategy for wind farm optimization, yet its effectiveness depends on the accuracy of underlying aerodynamic and structural models. In this study, we evaluate the predictive capabilities of models with varying fidelity for the IEA 22 MW reference turbine, considering both a single turbine and a two-turbine row with 5D spacing under conventionally neutral atmospheric boundary layer conditions. Results are benchmarked against large-eddy simulations (LESs). All models reproduced qualitative trends in power and, where applicable, loads as a function of yaw angle and downstream position, but there was a large spread in quantitative agreement. The dynamic wake meandering (DWM) model implemented in Dynamiks gave very good predictions for mean power, acceptable results for blade and yaw bending damage equivalent loads (DELs) but heavily underpredicted the tower bottom DELs compared to LESs. RANS results from EllipSys3D resolved asymmetric wake features but with reduced magnitude, leading to increasing errors for power prediction with increasing wake deflection. Steady-state engineering models (PyWake and Fuga) performed reasonably well for power prediction in the aligned cases but showed increasing errors under yaw misalignment. None of the engineering models reproduced secondary steering. These findings highlight the limitations of the tested engineering and mid-fidelity models, and emphasize the need for improved treatment of wake asymmetry, veer effects, and meandering physics to enhance reliability in practical optimization applications.

1 Introduction

Wind turbines extract kinetic energy from the incoming wind flow, which creates wakes downstream of each turbine. These wakes are characterized by a velocity deficit and increased turbulence intensity (TI) (Porté-Agel et al., 2020). If a turbine is placed within the wake of another turbine, it produces less power than it would under undisturbed inflow conditions. In addition, the elevated turbulence within the wake leads to higher relative loading on the downstream turbine compared to the case of a simple reduction in wind speed without added turbulence. Furthermore, if only part of the downstream turbine is exposed to the wake, the structural loading can increase even further compared to a fully waked turbine (Herges et al., 2018).

Therefore, the relative placement of turbines is a key factor in wind farm design. More recently, different wind farm control strategies have been proposed that can be used during the farm's operation. These strategies aim to increase the total power output of the farm and/or reduce the structural loading on individual turbines or the farm as a whole (Kheirabadi and Nagamune, 2019).

Static yaw control is one of the more promising control strategies that interests both industry and academia. The basic principle of static yaw control is the deflection of the upstream wake away from the downstream turbine using a predefined yaw misalignment angle to reduce the amount of waked inflow that the second turbine experiences. This can potentially increase the power production of the downstream turbine, but it can also negatively affect the loads on both turbines. Hence, the yaw misalignment angle must be carefully chosen to ensure that the increase in power from the downstream turbine compensates for the loss in power of the upstream turbine and that the changes in loading are acceptable (Meyers et al., 2022). Indeed, how exactly the turbine loading is affected by flow control methods such as static wake steering, how large the yaw offset should be, and for which flow conditions does static flow control outperform other flow control methods are still active areas of research. To answer those questions, aerodynamic models that can represent the effect of yawed inflow on both the rotor and the wake are necessary.

From an aerodynamic perspective, the flow behind a yawed wind turbine is rather complex. The misaligned inflow leads to a reduced thrust coefficient and a deflection of the wake, and the wake shape itself changes from circular to kidney shaped for large yaw angles (Schottler et al., 2018). Close analysis of the cross-flow field in the wake of a yawed wind turbine with uniform inflow reveals two counter-rotating vortices at the top and bottom of the wake. These lead to the deformation from a circular to a curled wake, as visualized experimentally and also numerically by Mikkelsen (2004) and Howland et al. (2016). However, Vollmer et al. (2016) showed that atmospheric stability and the associated veer strongly influence the curling strength. These researchers analyzed a single NREL 5 MW turbine under different stability conditions and yaw angles. For a yaw misalignment of |γ|=30° and at six rotor diameters downstream: (i) in a neutral boundary layer with around 2° of veer over the rotor disk, the kidney shape is still visible and the magnitude of the vertical deflection relative to the tower position is similar for ±30° yaw; (ii) in a stable boundary layer with about 8° of veer, the kidney shape is also visible, but the vertical deflection differs strongly between −30° and +30° yaw due to the strong veer; (iii) in a convective boundary layer, the wake appears largely unaffected by the yaw misalignment.

The asymmetry of the wake with a streamwise and lateral deficit also affects the downstream turbine. Fleming et al. (2018) demonstrated that the wake of a turbine aligned with the incoming wind and downstream of a yawed turbine also mildly deflects, and coined the term “secondary wake steering”. A likely cause of this is that, on average, the lateral components of the deflected wake of the upstream turbine caused the downstream turbine to experience a de facto yawed inflow. Fleming et al. (2018) further noted that for a farm with multiple steered turbines, lateral deficits can also merge and further complicate the wake evolution.

At the early design stage of a wind farm, high-fidelity computational fluid dynamics (CFD) simulations are generally not computationally feasible for comparing flow control strategies or optimizing yaw setpoints. This limitation is particularly restrictive for load predictions, where time-series data from large-eddy simulations (LESs) would be required. Time-averaged Reynolds-averaged Navier–Stokes (RANS) simulations are more tractable for power predictions, yet even these remain inefficient for early-stage design studies. Nevertheless, for specific farm layouts, surrogate or reduced-order models trained on LES data have been successfully applied to wake steering and axial induction control (Hulsman et al., 2020; Debusscher et al., 2022).

For non-specific layouts, however, more computationally efficient approaches are necessary. To this end, a range of engineering wake models have been developed in the literature. These can be broadly divided into two categories: (i) models that rely on a symmetric circular wake description combined with a deflection model, sometimes with a projection into deflected coordinates (Jiménez et al., 2010; Bastankhah and Porté-Agel, 2016; Shapiro et al., 2018; Liew et al., 2023); and (ii) models that propose new deficit formulations explicitly accounting for wake curling (Martínez-Tossas et al., 2019; King et al., 2021; Branlard et al., 2022). In addition, Zong and Porté-Agel (2020) introduced an iterative momentum-conserving wake superposition method for streamwise and lateral deficits. Using the deflection model of Shapiro et al. (2018), they accurately reproduced both primary and secondary steering effects for wind-tunnel test turbines. As an alternative to traditional engineering models, linearized RANS approaches have also been applied to wake modeling and, more recently, adapted to model wake deflection (Larsen et al., 2020).

To date, several publications have compared varying-fidelity models against high-fidelity CFD, LiDAR, and SCADA data. For example, Larsen et al. (2020) compared RANS, linearized RANS, DWM, and engineering model flow-field results with LiDAR measurements for a standalone V52 turbine. In the FarmConners benchmark, Göçmen et al. (2022) conducted a series of blind benchmark tests for wake steering, where lower-fidelity models were compared against both mean and time series LESs and SCADA data for single- and multi-turbine wakes with yaw offsets. These comparisons were quantified in terms of power predictions. The novelty of the present work lies in the level of detail of the comparison – covering flow, power, and loads – as well as in the breadth of the models considered. Furthermore, to date, in the literature, there are very few examples of a detailed load validation of the DWM against LESs or measurement data for very large turbines. For load comparison against measurements, to the author's knowledge, the largest turbine used is a 6 MW (Bernard et al., 2024). For load comparison against LES, there are a few on the IEA 15 MW turbines (Branlard et al., 2024; Rivera-Arreba et al., 2024) but none on the IEA 22 MW turbine. This is important because larger turbines are more sensitive to the influence of shear and veer.

This paper aims to benchmark the whole toolchain of aerodynamic models developed at DTU against high-fidelity large-eddy simulations (LESs) for the large IEA 22 MW reference turbine for primary as well as secondary wake steering at different angles and also veered inflow for a conventionally neutral boundary layer (CNBL). The performance of the models is evaluated for power and, for the dynamic models, also loads. The array of tested models, ranging from lower to higher fidelity, is as follows: (i) the engineering models implemented in PyWake (Pedersen et al., 2023) and Fuga (Ott et al., 2011), (ii) a novel implementation of the dynamic wake meandering (DWM) model in the open-source Dynamiks framework (Pedersen et al., 2026), and finally (iii) the Reynolds-averaged Navier–Stokes (RANS) solver implemented in EllipSys3D (van der Laan et al., 2024).

Among the range of models considered in this study, the LESs resolve the highest level of physical detail relevant to atmospheric turbulence and turbine wake dynamics. For this reason, the LES results are used as a reference for comparison with the lower-fidelity models. In addition, the LES precursor simulations provide a consistent set of inflow planes that can be imposed across all models. This reduces uncertainty associated with inflow specification and allows differences between model results to be primarily attributed to differences in modeling assumptions rather than differences in inflow conditions, thereby enabling a controlled comparison between the models. An alternative reference could have been field measurements; however, such data introduce additional uncertainties related to inflow characterization and atmospheric variability. Although it is acknowledged that the LES does not capture all aspects of atmospheric physics present in measurements, it does resolve the range of physical processes that the lower-fidelity models examined here are intended to represent. Consequently, the use of measurements or meso–micro coupled simulations would not substantially alter the conclusions of the present comparison (Doubrawa et al., 2020).

Section 2 explains the setup and details the methodology used for all the models. Section 3 contains a detailed comparison between the models in terms of power and loads, where applicable. Finally, the conclusions are given in Sect. 5.

2 Methodology

This section describes the methodology behind the setup in Sect. 2.1 and all the models in the remaining section.

2.1 Setup

Two setups are considered: first, a single turbine with different prescribed yaw offset angles; second, two turbines in a row with two different yaw setpoints for the first turbine and no yaw offset angle for the second turbine. The main parameters of the IEA22 MW turbine are specified in Table 1, and additional details are described by Zahle et al. (2024). The coordinate system is defined by the x axis being aligned with the undisturbed flow direction at hub height, while the z axis points upward. With this convention, the sign of the yaw angle is opposite to the right-handedness of the coordinate system about the z axis.

Table 1 Turbine and setup properties.

2.2 LES

The LES results are obtained using EllipSys3D, a finite-volume incompressible Navier–Stokes solver initially developed by Michelsen (1992), Michelsen (1994), and Sørensen (1994), which solves the governing equations (including the potential temperature equation, as a conventionally neutral boundary layer (CNBL) is used as inflow) with a finite volume method in general curvilinear coordinates using a collocated grid. The pressure correction equation is solved with an extended SIMPLEC approach (Shen et al., 2003; Kobayashi and Pereira, 1991) using Rhie/Chow interpolation. Time advancement uses a second-order three-level implicit method with sub-iterations in each time step. Convective terms are discretized with a fourth-order central differencing scheme with a fourth-order dissipation term to reduce numerical instabilities (Wit and van Rhee, 2013). Sub-grid scale modeling is achieved using the AMD model (Abkar et al., 2016).

The inflow for all cases is based on a CNBL generated using an LES precursor simulation, conducted in a domain with size $\mathrm{16}\mathrm{640}\times \mathrm{16}\mathrm{640}\times \mathrm{3000}$ m, with grid resolution $\mathrm{d}x=\mathrm{d}y=\mathrm{4}\mathrm{d}z=\mathrm{20}$ m. The applied initial conditions are ground temperature T_wall=288.15 K; an initial temperature profile constant until 600 m height with perturbations up to 100 m height to encourage breakdown into turbulence, then a linear gradient of 10⁻³ K m⁻¹ above; geostrophic wind values U=9.48 m s⁻¹ and $V=-\mathrm{0.59}$ m s⁻¹; velocity profiles initialized with a log-law based profile with 600 m boundary layer height (Allaerts and Meyers, 2015); Coriolis parameter $f_{\mathrm{c}}=\mathrm{1.1947}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹; and wall roughness z₀=0.001 m. A wind direction controller (Sescu and Meneveau, 2014) is active in order to achieve aligned flow at turbine hub height (170 m). Rayleigh damping is applied above a height of 2000 m (Klemp and Lilly, 1978). The precursor is run for 26 h while the boundary layer develops, followed by a further period of 4500 s over which the three velocity components and the temperature are extracted over a cross-stream plane in the domain center, which are used as inflow to the successor simulation. The inflow profile is visualized in Fig. 1.

In the successor simulations, wind turbines are modeled using the actuator disk (AD) method (Mikkelsen, 2004), which is fully coupled to the aero-servo-elastic solver HAWC2 (Larsen and Hansen, 2025) (which has a Timoshenko beam element-based multi-body formulation and hence can account for non-linear effects) through the Dynamiks coupling framework (Pedersen et al., 2026). The velocity components are extracted along the positions of the three rotating blades in EllipSys3D and passed to HAWC2, which calculates loads and deflections. These are transferred back to EllipSys3D, and used to set the magnitude and position of the body forces that make up the AD. The AD body forces are applied in three overlapping 240° sections, where forcing decreases from a maximum at the blade location to zero at the neighboring blades, which means that non-uniform loading can be captured. The two-way nature of the coupling means that the interaction between loading, deflections, and flow is also represented (Hodgson et al., 2021). HAWC2 also uses the DTU Wind Energy Controller, so the turbine varies rotational speed and pitch in response to its loading.

Successor simulations are conducted in a domain with size $\mathrm{8520}\times \mathrm{4970}\times \mathrm{3000}$ m ($\mathrm{30}D\times \mathrm{17.5}D\times \mathrm{10.6}D$), with a refined region stretching from the inlet to 20D downstream, laterally extending 3D in each direction from the domain centerline, and from the ground to 4D in height. Within the refined region, the grid resolution is $\mathrm{d}x=\mathrm{d}y=\mathrm{d}z=D/\mathrm{32}$, and outside the cell size is stretched to the domain boundaries. The first turbine is placed 3.5D from the inlet, and for the two-turbine setup, a second turbine is placed 5D behind the first. Each simulation has 900 s, which are used as spin-up, with the final 3600 s used for analysis.

The modeling of atmospheric boundary layers, wind turbines, and wakes in EllipSys3D has been extensively verified and validated through comparisons with measurement data (Mikkelsen, 2004; Doubrawa et al., 2020), fully resolved blade simulations (Hodgson et al., 2022), and other LES codes (Hodgson et al., 2023). Like all numerical models, LES includes uncertainties, in particular related to grid resolution and sub-grid scale model. Hence, the chosen grid resolutions for both precursor and successor simulations are informed by verification studies and cross-code comparisons (Hodgson et al., 2023) in order to ensure the grid convergence of important statistics (atmospheric boundary layer (ABL) and wake profiles, turbine thrust coefficient). The constant geostrophic wind used in the precursor may lead to lower wind direction variability – at least in long-term statistics – than observed in reality, where the driving atmospheric conditions are more variable. However, the setup captures much of the physics that dictate wind turbine wake behavior and recovery: the turbulent spectra, interaction with the Coriolis force and veer, atmospheric stability, dynamics of the near and far wake, and realistic turbine loading and response through the aeroelastic coupling.

2.3 RANS

The RANS results are computed by PyWakeEllipSys v5.4 (DTU Wind and Energy Systems, 2025), which uses a steady-state version of EllipSys3D. The k-ε-f_P turbulence model (van der Laan et al., 2015b) is employed, coupled with an atmospheric boundary layer (ABL) inflow model including Coriolis forces and a simplified turbulent buoyancy source term that sets an ABL height implicitly (van der Laan et al., 2024). The latter is used to represent the stratification of the ABL, where lowering the ABL height results in a more stable ABL characterized by a larger wind shear and veer with respect to neutral conditions. The simplified turbulent buoyancy source term is only dependent on a prescribed Brunt-Väisälä frequency 𝒩, leading to an inflow model that does not require a temperature equation. The inflow is generated by a precursor simulation using EllipSys1D (van der Laan and Sørensen, 2017); the results are compared with LES in Fig. 1. The hub-height wind speed and turbulence intensity based on k from LES (9.2 m s⁻¹ and 3.8 %) are obtained using a geostrophic wind speed G=9.47 m s⁻¹ and $\mathcal{N}=\mathrm{1.22}\times {\mathrm{10}}^{-\mathrm{3}}$ s⁻¹. The roughness length is set lower than LES, $z_{\mathrm{0}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ m, to obtain a similar ABL height.

The turbine is modeled as rigid using an actuator disk (Réthoré et al., 2014) with an analytic force distribution model of Sørensen et al. (2020) that includes tangential forces and local effects of shear and veer. A precomputed lookup table controls the AD aerodynamic coefficients as a function of the disk-averaged wind speed using single AD simulations without yaw and tilt (van der Laan et al., 2015a). The RANS simulations with the controller include a tilt angle of 6°and optional yaw angles.

The numerical domain of the single-turbine simulations is a Cartesian box with dimensions ($L_x,L_y,L_z)=(\mathrm{1025}D,\mathrm{1006}D,\mathrm{25}D)$. The turbine is placed near the horizontal center, $(x,y,z)=(\mathrm{0},\mathrm{0},z_{\mathrm{H}})$, and the flow around the turbine is resolved using a refined inner domain with dimensions $-\mathrm{5}D<x<\mathrm{20}D,\left|y\right|<\mathrm{3}D,z<\mathrm{3}D$ using a minimum spacing of 16D and a first cell height of $D/\mathrm{200}$. A total of 12.6 million cells are used. The refined and outer domains of the two-turbine case is 5D longer in the streamwise and lateral directions, leading to a total of 22.4 million cells. The bottom boundary is a rough wall boundary condition (Sørensen et al., 2007). The inflow and top boundaries are inlet conditions at which the inflow profiles are prescribed (hence the flow is lid driven). The lateral boundaries are periodic, and the outlet is a boundary condition at which all gradients in the streamwise direction are assumed to be zero.

2.4 DWM

The new implementation of the DWM in the open-source Dynamiks framework is accessible under Pedersen et al. (2026) and is heavily based on HAWC2Farm as benchmarked by Liew et al. (2023), but it is much more modular. The DWM implementation in Dynamiks should not be confused with the older version implemented in HAWC2 itself (Larsen and Hansen, 2025).

Given a background flow U_ambient, a turbulence box covering the wind farm U_turb and an array of n_WT wind turbines with wake deficits (including the wake-added turbulence) ΔU_i, the velocity at any point in the domain is given by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{U}(\mathit{x},t)& ={\mathit{U}}_{\text{ambient}}(\mathit{x},t)+{\mathit{U}}_{\mathrm{turb}}(\mathit{x},t)\\ & +f_{\text{Superposition}}\left(\mathrm{\Delta }{\mathit{U}}_{\mathrm{1}}(\mathit{x},t),\mathrm{\dots },\mathrm{\Delta }{\mathit{U}}_{n_{\mathrm{WT}}}(\mathit{x},t)\right).\end{array}\end{array}$$

Different wake superposition methods are available in Dynamiks. This publication uses the root-of-summed-squares superposition method with the rotor-averaged inflow velocity as a reference for the deficit calculation (Zong and Porté-Agel, 2020).

The base assumption behind the dynamic wake meandering model is that a wind turbine wake behaves like a passive tracer advected by the large-scale turbulence. The cutoff frequency for the turbulent structures is defined as structures larger than twice the local wake width, which was conjectured in Larsen et al. (2008) and subsequently full-scale validated in Larsen and Lio (2025). At approximately each time step, the turbine emits a tracer particle that is then advected based on the filtered background flow. In this publication, we use a spatial averaging filter with the local wake diameter D_w(x) on the background flow field to get the advection velocity U_advection for a particle k emitted from turbine j with the location x_j,k:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{U}}_{\text{advection},j,k}& =\sum _{i=\mathrm{0}}^{n_{\mathrm{avg}}}{w}_i\stackrel{\mathrm{̃}}{\mathit{U}}({\mathit{x}}_{j,k}+D_{\mathrm{w}}(x)\cdot \mathrm{\Delta }{\mathit{x}}_i)\\ & \text{for}k\in [\mathrm{1},n_{\mathrm{particles}}]\text{and}j\in [\mathrm{1},n_{\mathrm{WT}}],\end{array}\end{array}$$

where n_avg, w_i, and Δx_i are the number, the weight, and the location of the stencil points for the spatial averaging filter. Each particle has a wake deficit attached to it that is initialized with the deficit at the rotor when the particle was emitted, and the wake evolves as the particles move downstream. Different settings are available for the background flow $\stackrel{\mathrm{̃}}{\mathit{U}}$, but the settings that make the most sense are the local flow, including either all deficits or just the deficits upstream of the turbine from which the particle was emitted. Note that since the veer is included in the background flow, this also leads to a deflection of the particles due to the veer.

The advection velocity is further processed if the turbine is yawed or tilted. At the moment, two deflection models are available in Dynamiks, namely (i) the Hill Vortex model from Larsen et al. (2020) and (ii) the Jiménez model from Jiménez et al. (2010). In this publication, we use the Hill Vortex model because it gave better results than the Jiménez model, even when tuning the calibration parameter of the model for specific yaw angles. For the wake deflection with a yaw misalignment γ and a tilt angle θ, the model writes as follows:

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathit{U}}}_{\mathrm{advection},j,k}& {\displaystyle }={\mathit{U}}_{\mathrm{advection},j,k}-\mathrm{0.4}\\ \text{(3a)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{\Delta }{U}_{\mathrm{w},j,k}\left[\begin{array}{c}-\cos \left(\mathit{\gamma }\right)\cdot \cos \left(\mathit{\theta }\right)\\ \sin \left(\mathit{\gamma }\right)\cdot \cos \left(\mathit{\theta }\right)\\ \sin \left(\mathit{\theta }\right)\end{array}\right]\text{(3b)}& {\displaystyle }\mathrm{\Delta }{U}_{\mathrm{w},j,k}& {\displaystyle }={\displaystyle \frac{\mathrm{2}}{R_{\mathrm{w}}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{R_{\mathrm{w}}}{\int }}\mathrm{\Delta }{U}_j\left({\mathit{x}}_k\right)r\mathrm{d}r,\end{array}$$

where k indexes the particles belonging to turbine j with coordinates x_j,k, and R_w is the wake expansion at the particle's location. For evaluating the advection velocity, the wake of the turbine from which the particle is emitted is excluded.

The most commonly used wake deficit model in the DWM context is the Ainslie model. It is derived from the steady-state RANS equations in cylindrical coordinates, assuming (i) symmetry with respect to rotation such that $\frac{\partial }{\partial \mathit{\theta }}=\mathrm{0}$ and ${\stackrel{\mathrm{̃}}{U}}_{\mathit{\theta }}=\mathrm{0}$, (ii) that the shear layer between the wake and freestream is thin such that the gradients of the mean flow fields are much bigger in the radial than in the axial direction, (iii) that the pressure gradient in the axial direction is zero which only holds in the far wake, and (iv) that the Reynolds stresses can be modeled using a simple eddy viscosity model. The momentum equation in the axial direction and the mass conservation equation for a rotor-averaged inflow wind speed U_rotor will then look as follows:

$$\begin{array}{}\text{(4a)}& {\displaystyle }\stackrel{\mathrm{̃}}{U}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}}{\partial x}}+{\stackrel{\mathrm{̃}}{V}}_r{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}}{\partial r}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left({\mathit{\nu }}_{T}r{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}}{\partial r}}\right)\text{(4b)}& {\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\stackrel{\mathrm{̃}}{V}}_r\right)+{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}}{\partial x}}& {\displaystyle }=\mathrm{0},\end{array}$$

with the boundary conditions

$$\begin{array}{}\text{(5a)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}(x,r)}{\partial r}}{|}_{r=\mathrm{0}}& {\displaystyle }=\mathrm{0}\text{(5b)}& {\displaystyle }\underset{r\to \mathrm{\infty }}{lim}\stackrel{\mathrm{̃}}{U}(x,r)& {\displaystyle }=U_{\mathrm{rotor}}\text{(5c)}& {\displaystyle }\stackrel{\mathrm{̃}}{U}(\mathrm{0},r)& {\displaystyle }={\stackrel{\mathrm{̃}}{U}}_{\mathrm{0}}\left(r\right).\end{array}$$

The initial velocity deficit at the rotor ${\stackrel{\mathrm{̃}}{U}}_{\mathrm{0}}\left(r\right)$ can be calculated from a BEM model or a simple actuator disk model. Dynamiks has an interface to the BEM model implemented in HAWC2 and also the simple actuator disk turbines from PyWake. This publication uses the DWM model in combination with HAWC2 only. Since the pressure term in the momentum equation is dropped, the wake expansion and accompanying velocity deceleration in the near wake are not captured. To counteract this, the output of the actuator model is scaled before using it as an inlet boundary condition for the Ainslie model:

$$\begin{array}{}\text{(6a)}& {\displaystyle }{\stackrel{\mathrm{̃}}{U}}_{\mathrm{0}}\left(r\right)& {\displaystyle }=U_{\text{rotor}}\left(\mathrm{1}-(\mathrm{1}+f_U)a\left(\stackrel{\mathrm{̃}}{r}\right)\right)\text{(6b)}& {\displaystyle }r& {\displaystyle }=\stackrel{\mathrm{̃}}{r}\sqrt{{\displaystyle \frac{\mathrm{1}-a}{\mathrm{1}-(\mathrm{1}+f_R)a}}},\end{array}$$

where U_rotor is the rotor-averaged inflow wind speed, $a\left(\stackrel{\mathrm{̃}}{r}\right)$ is the radially varying induction obtained from the actuator model, and f_R=1.1 and f_U=0.98 are the tuning parameters for the wake expansion and velocity deceleration, respectively. Keck et al. (2012) claim that this makes the model accurate from 3D downstream.

The velocity deficit ΔU_j of the turbine j that is used in the superposition is then given by $\mathrm{\Delta }{U}_j=U_{\text{amb},j}-{\stackrel{\mathrm{̃}}{U}}_j$.

A simple mixing length model is used for the eddy viscosity as developed by Keck et al. (2012). It consists of two terms accounting for ambient and wake-induced turbulence, respectively.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\nu }}_T& =F_{\mathrm{amb}}{k}_{\mathrm{amb}}{\text{TI}}_{\mathrm{amb}}+F_{\mathrm{wake}}{k}_{\mathrm{wake}}\\ & \cdot max\left({\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{D}_{\mathrm{w}}\left(x{)}^{\mathrm{2}}\right|{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}}{\partial r}}|,{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{D}_{\mathrm{w}}(x\left)\right(\mathrm{1}-\underset{r}{min}\left(\stackrel{\mathrm{̃}}{U}\right(x,r\left)\right))\right)\end{array}\end{array}$$

Here, D_w(x) is the local wake diameter; TI_amb is the turbulence intensity at the rotor, the parameters k_amb=0.0914 and k_wake=0.0216; and finally

$$\begin{array}{}\text{(8a)}& {\displaystyle }{F}_{\text{amb}}\left(x\right)=\left\{\begin{array}{l}\frac{x}{\mathrm{4}},\text{for }x<\mathrm{4}\\ \mathrm{1},\text{for }x\ge \mathrm{4}\end{array}\right.\text{and}\text{(8b)}& {\displaystyle }{F}_{\mathrm{wake}}\left(x\right)=\left\{\begin{array}{l}\mathrm{0.035},\text{for }x<\mathrm{4}\\ \mathrm{1}-\mathrm{0.965}\exp \left[-\mathrm{0.35}\left(\frac{x}{\mathrm{2}}-\mathrm{2}\right)\right],\text{for }x\ge \mathrm{4}\end{array}\right.\end{array}$$

are filter functions that model the delay in turbulent wake diffusion due to entrainment and compensate for the non-equilibrium turbulence close to the rotor. The influence of turbulence build-up and atmospheric shear, as described by Keck et al. (2015), are not implemented yet.

For modeling the wake-added turbulence (WAT), a homogeneous and isotropic Mann box is used. The box size is given by $L_x\times L_y\times L_z=[\mathrm{15}D,D,D]$, and the spatial discretization is $N_x\times N_y\times N_z=[\mathrm{1024},\mathrm{64},\mathrm{64}]$. The box uses periodic boundary conditions and the same advection speed as the background turbulence. The Mann box is fully defined by the anisotropy parameter Γ=0, the length scale $L=D/\mathrm{16}$, and unit variance σ²=1. To obtain a unit variance, the box is rescaled to give $\mathit{\alpha }{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}$ according to Mann (1994):

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}^{\mathrm{2}}\approx \mathrm{0.688}\mathit{\alpha }{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{L}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where α is an empirical parameter and ϵ is the dissipation of TKE.

The local velocity of the box is then multiplied by a scaling parameter $k_{m_t}$ as first proposed by Madsen et al. (2010):

$$\begin{array}{}\text{(10)}& k_{m_t}\left(r\right)=\left|\mathrm{1}-\stackrel{\mathrm{̃}}{U}\left(r\right)\right|k_{m\mathrm{1}}+\left|{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{U}\left(r\right)}{\partial r}}\right|k_{m\mathrm{2}},\end{array}$$

where the tuning parameters k_m1 and k_m2 are adapted from Branlard et al. (2024) as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }f(\stackrel{\mathrm{̃}}{x},f_{min},D_{min},D_{max},e)=f_{min}\\ \text{(11a)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}-f_{min}){\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{x}-D_{min}}{D_{max}-D_{min}}}\right)}^e\text{(11b)}& {\displaystyle }& {\displaystyle }{k}_{\text{def}}\left(x\right)=\mathrm{0.6}f(\stackrel{\mathrm{̃}}{x},\mathrm{0},\mathrm{0},\mathrm{2},\mathrm{1})\text{(11c)}& {\displaystyle }& {\displaystyle }{k}_{\text{grad}}\left(x\right)=\mathrm{3.0}f(\stackrel{\mathrm{̃}}{x},\mathrm{0},\mathrm{0},\mathrm{12},\mathrm{0.65}).\end{array}$$

For the simulations presented in this paper, the LES precursor simulations are used to generate the background turbulence and mean flow for the most direct comparison. The precursor inflow planes are stacked together, and the mean velocity at hub height is the advection speed of the box. The wake-added turbulence box is computed as described above using one seed. As in the LESs, the same HAWC2 setup file and controller are used for the structural turbine model.

2.5 PyWake

Steady-state flow conditions for various prescribed yaw angles were simulated using PyWake (Pedersen et al., 2023). The inflow conditions were selected to match the hub-height wind speed and turbulence intensity obtained from LES results. The simulation does not include blockage effects and propagates downstream a wake deficit modeled using the super-Gaussian model, with the recalibration according to Blondel (2023). Wake-generated turbulence is added on top of the prescribed background turbulence intensity using the model from Crespo and Hernández (1996). A root-of-summed-squares superposition method is used to combine multiple wakes. The wake deflection is modeled using the Jiménez model with β=0.1 (Jiménez et al., 2010), an engineering wake deflection model empirically derived from LES data that accounts for lateral displacement of wakes due to turbine yaw.

2.6 Fuga

The turbine wake evolution is modeled using a fast linearized CFD solver based on the Reynolds-averaged Navier–Stokes (RANS) equations with an eddy viscosity closure. Fuga (Ott et al., 2011; Ott and Nielsen, 2014) assumes that the forcing from the turbine is small relative to the background shear flow and applies a first-order perturbation expansion. The governing equations are solved efficiently using a mixed spectral formulation, where a Fourier transform in the horizontal directions reduces the system to a set of decoupled ordinary differential equations in the vertical coordinate. Lookup tables are precomputed and reused to account for different turbine positions and forcing without solving the full equations for each new case.

Non-linear effects of yaw-induced wake deflection are captured by adapting a linearization technique (Ott et al., 2019; Larsen et al., 2020). This involves a curvilinear coordinate transformation in which the transverse direction is shifted by a small displacement $\mathit{\lambda }(x,y,z)$ such that the new transverse coordinate $\stackrel{\mathrm{´}}{y}=y-\mathit{\lambda }$ remains constant along streamlines. This transformation effectively removes the wake deflection in the new coordinates, only allowing it to be recovered upon returning to the original space. The displacement λ is computed from the lateral velocity field and the background streamwise velocity, and stored as an additional lookup table.

The inflow conditions were designed to match the mean wind profile obtained from LESs, using a target wind speed at a hub height of U_hub=9.2  m s⁻¹ and a prescribed surface roughness length of z₀=0.001 m. These conditions are used to define the background velocity profile U⁰(z) and the eddy viscosity K(z) via Monin–Obukhov similarity theory, with the stability function for momentum ${\mathit{\varphi }}_m(z/L)$ set for neutral conditions.

Since Fuga is linear by definition, wake superposition is also linear. This means the results of superimposed deflected wakes might be slightly affected by breaking the assumption of null deflection at the wind turbine rotor. In our case, the downstream turbine is not yawed, so the effect of this simplification on the resulting wind fields is negligible.

2.7 Definition of different turbine representations

Along with an array of flow models, a corresponding array of turbine representations are presented in this paper. For better readability, in this section, we define different terms for the different turbine representations that then correspond with the captions in the plot:

• HAWC2 turbine: A power-extracting turbine that is modeled with rotor forces and induces a wake. Its loads and power are computed from the modeled rotor response (e.g., via a two-way coupling between HAWC2 and LES or a one-way coupling between HAWC2 and the DWM model).

• HAWC2 ghost turbine: A virtual turbine whose loading is evaluated by running a standalone HAWC2 simulation using a series of inflow planes sampled from a regular simulation as input. There is no coupling between the inflow and the turbine, the induction model in HAWC2 is used to approximate the induction. For the DWM, ghost turbines should give the same result as simulations where turbines are included directly, because the DWM is using the BEM-based induction model in HAWC2 anyway. For the LESs, ghost turbines give an approximation of the loading of a real turbine.

• RANS LUT model: For the RANS simulations, an actuator disk is used to represent the turbine. The thrust of the turbine is evaluated from a lookup table and the oncoming wind speed in a two-way coupled convergence process.

• Velocity + power-curve: A turbine whose power is evaluated on the basis of rotor-averaged inflow velocity rather than a HAWC2 simulation.

The virtual turbines were introduced to get more value out of the dataset since they are much more computationally efficient and several scenarios can be evaluated. Further, ghost turbines are also used to train load surrogate models (Hulsman et al., 2020). A comparison between the results of the virtual and the “real” turbine is also presented.

3 Results

The sections contain a comparison of the inflow profiles in Sect. 3.1, the results for the single-turbine setup in Sect.3.2, and finally, the results for the two-turbine setup in Sect. 3.3.

3.1 Inflow

Care is taken to match the inflow profiles for mean velocity, inflow direction, and turbulence intensity to the LES precursor. The specific matching procedures are described in the methodology sections for each model. For steady-state models (RANS, Fuga, PyWake), the LES inflow is averaged over the 1 h simulation window. Section 1 shows the results at the inlet of the domain.

The wind direction controller in the LES precursor setup did not fully converge, so there is a slight offset in the wind direction at hub height of $\mathrm{\Delta }\mathit{\varphi }=-\mathrm{0.5}\mathit{°}$. Averaged over the rotor disk, this becomes $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathit{\varphi }}=\mathrm{0.4}\mathit{°}$. This offset in wind direction was applied in the remaining model simulations to ensure consistency between the setups. Initial calculations neglected this offset because it appeared small. However, rotating the reference coordinate system by 0.5° leads to an offset of the wake center as compared to the unrotated coordinate system of $\mathrm{\Delta }{y}_{\mathrm{WC}}/D(x/D=\mathrm{15})=\tan \left(\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathit{\varphi }}\right)\cdot \mathrm{15}\approx \mathrm{0.1}$, and as visible in Fig. 5 the order of magnitude is similar due to the one for wake steering. Since Fuga and PyWake do not model veer, the wind direction remains constant with height.

For the comparison of the turbulence intensity (TI) to LES, total and streamwise formulations are used for RANS and the DWM, respectively:

$$\begin{array}{}\text{(12a)}& {\displaystyle }& {\displaystyle }\text{TI}={\displaystyle \frac{\sqrt{\mathrm{2}/\mathrm{3}k}}{U_{\mathrm{ref}}}}={\displaystyle \frac{\sqrt{\mathrm{1}/\mathrm{3}({\mathit{\sigma }}_U^{\mathrm{2}}+{\mathit{\sigma }}_V^{\mathrm{2}}+{\mathit{\sigma }}_W^{\mathrm{2}})}}{U_{\mathrm{hub}}}}\text{(12b)}& {\displaystyle }& {\displaystyle }{\text{TI}}_u={\displaystyle \frac{{\mathit{\sigma }}_U}{U_{\mathrm{hub}}}}.\end{array}$$

For the RANS model, the inflow turbulence is isotropic, which does not correspond to reality. This means that the RANS model underpredicts the streamwise TI, but the total TI is close to the LES profile, since the TKE is fitted to match the LES. Hence, the total TI is used for comparison with LES for the RANS model. For the DWM model, only the streamwise part of the TI is fully approximated; hence, the streamwise TI is used for comparison to LES.

Figure 1 Comparison of inflow profiles from different models matched to the LES precursor.

Table 2 Inflow properties.

3.2 Setup I: a single wind turbine

First, a visual comparison of the model predictions through horizontal slices of the streamwise velocity at hub height for the fully aligned case and the misaligned case with $\mathit{\gamma }=-\mathrm{30}\mathit{°}$ is shown in Fig. 2. For both the aligned and the misaligned cases, the prediction of the far wake deficit and wake recovery rate seems reasonably well captured by all the models. In addition, the lateral wake expansion rate appears similar across all models for both cases (see Appendix C for deficit plots for the aligned case), indicating that differences in wake width are not the primary source of the discrepancies discussed below. For the near wake, there is more variation. The RANS model is close to the LES. The DWM model overpredicts the deficit close to the rotor per the boundary conditions for the deficit model; this neglects the gradual pressure recovery after the rotor but is close to the LES mean from about 2D onwards. The linearized RANS solution from Fuga underpredicts the deficit in the near wake but starts to agree well with LES from about 5D onwards. At first glance, the velocity field at the rotor location in Fuga suggests that the turbine yaws in the wrong direction. However, this is simply an artifact of how the turbine forcing is implemented in Fuga: it uses a non-yawed disk that applies forcing in both the longitudinal and transversal directions. When looking further upstream from the turbine, the blockage pattern shows that the rotation occurs in the correct direction. Finally, the engineering model implemented in PyWake predicts the flow deceleration behind the rotor due to wake expansion accurately in terms of magnitude, but the transition from near to far wake comes too late, and the deficit starts to agree with the LES from about 6D onwards. Since the mean wind direction misalignment over the rotor disk points downward and the wake moves upwards due to tilt in an area where the mean veer also points downwards, the wake should deflect a bit even in the aligned case. Fuga and PyWake do not model the veer, so they can only model the part of that effect related to the mean wind direction misalignment.

Figure 2 Time-averaged streamwise velocity at hub height for all models for (left) the aligned and (right) the yawed case with $\mathit{\gamma }=-\mathrm{30}\mathit{°}$.

Moving on to the TI, the RANS, PyWake, and DWM fields are compared to the LES fields; the results from Fuga are not shown since it does not model turbulence added by the wake. As explained in Sect. 3.1, differing TI formulations are used for the DWM and RANS/PyWake. The top plots of Fig. 3 depict the total TI for the mean LES, RANS, and PyWake predictions. As expected, RANS overpredicts TI as compared to LES, but surprisingly, the overprediction is less pronounced for the yawed case, even though the deficit magnitude predictions are similar between the two. The WAT model in PyWake from Crespo and Hernández heavily overpredicts the TI as compared to LES and RANS. The bottom plots of Fig. 3 show the streamwise turbulence intensity TI_u for the mean LES and DWM predictions. With the updated calibration of the WAT model in the DWM from Branlard et al. (2024), at hub height, the DWM predictions are very close to the ones from LES, both in the near and far wake, even for the deflected case. The high TI region due to the tip and root vortices in the near wake persists for a bit longer in the LES than in the DWM. However, since the WAT formulation does not consider shear, below and above hub height, the TI is under- and overpredicted, respectively. In Appendix A in Fig. A1, streamwise slices of the TI and rotor-averaged TI in the wake of the turbine are shown for the DWM and the LES results.

Figure 3 Turbulence intensity (TI) at hub height across five columns: (I)–(III) total TI from LES, RANS, and PyWake; and (IV)–(V) streamwise TI from LES and the DWM model. On the left is the aligned case, and on the right is the yawed case with $\mathit{\gamma }=-\mathrm{30}\mathit{°}$.

To gauge more accurately how well the wake deflection is captured for different yaw misalignment angles, a wake-tracking algorithm was run on the (time-averaged) flow fields for all models. The Constant Area wake center tracking algorithm from the Quon (2025) was run on the streamwise velocity fields. Other simpler wake-tracking algorithms were tested, but they performed unreliably due to the asymmetry of the wake due to veer and curling. The waketracking algorithm identifies velocity deficit isocontours and looks for the contour line that roughly corresponds to the rotor area. Then a weighted average over the area is used to calculate the wake center (y_WC,z_WC) as

$$\begin{array}{}\text{(13a)}& {\displaystyle }{y}_{\text{WC}}\left(x\right)& {\displaystyle }={\displaystyle \frac{\iint y\mathrm{\Delta }U(x,y,z)\mathrm{d}y\mathrm{d}z}{\iint \mathrm{\Delta }U(x,y,z)\mathrm{d}y\mathrm{d}z}}\text{and}\text{(13b)}& {\displaystyle }{z}_{\text{WC}}\left(x\right)& {\displaystyle }={\displaystyle \frac{\iint z\mathrm{\Delta }U(x,y,z)\mathrm{d}y\mathrm{d}z}{\iint \mathrm{\Delta }U(x,y,z)\mathrm{d}y\mathrm{d}z}}.\end{array}$$

Figure 4 shows the streamwise velocity of the full wake at 5 diameters downstream, along with the results from the wake-tracking algorithm for all models. The variation in wake deficit shape shows that the wake center is a rather diffuse concept for asymmetric time-averaged wakes and cannot fully capture all effects at play in the complex flow field. For example, at that location, as will be shown in Fig. 5, while the RANS model underpredicts the wake deflection according to the wake-tracking algorithm, it captures all the wake asymmetry due to yaw and veer. In contrast, at that specific location, the DWM model and PyWake give the closest estimate of the wake center as compared to LES, but the wake shape looks completely different. As will be shown later, for the inflow conditions considered here, this proves adequate for power predictions, as the rotor-averaged velocity in the wake remains well matched. However, for the load predictions, matching the actual shape of the wake deficit becomes more important.

Figure 4 Streamwise time-averaged velocity at 5D downstream for all models for (top) zero, and (bottom) −30° yaw misalignment. The white circles outline the rotor of the upstream turbine.

Figure 5 shows the location of the wake center determined by the wake-tracking algorithm for all models and yaw angles.

Figure 5 Wake center position in (left) horizontal and (right) vertical as obtained from the wake-tracking algorithm of the time-averaged flow field for all models and yaw angles (γ) at veer of ϕ=6°.

For the lateral deflection due to yaw misalignment, all models capture at least the right sign for the deflection, but otherwise, there is a large spread in the results. RANS underpredicts the lateral deflection for all simulated angles except for the positive deflection angle γ=10°, where it mildly overpredicts the deflection as compared to LES in the far wake. The authors speculate that this might be due to shortcomings of the turbulence model in the near wake, which are more amplified in yawed versus non-yawed flow. Averaged over all deflection angles, the DWM model using the Hill Vortex deflection model and a spatial averaging filter for the particle advection gives the closest results to LES. Fuga predicts lateral deflection values similar to those of the RANS model, with a stronger underprediction of the deflection for negative yaw angles. Finally, the PyWake model using the Jiménez model gives more mixed results with no clear trend. For the largest yaw angle, both for the DWM and the PyWake results, the shape of the lateral deflection line is more curved in the near than the far wake, as compared to LES. The difference in vertical deflection shows an even larger spread between models than the lateral deflection. LES and RANS show similar trends with an upward deflection in the near wake and an up- or downward deflection in the far wake, depending on how substantial the yaw misalignment is. The vertical deflection in the DWM model is only driven by the tilt, and hence it shows a steady upward deflection, either under- or overpredicting as compared to LES in the far wake, with decent agreement in the near wake. Fuga seems to predict a downward deflection even for the non-yawed case for two reasons: (i) Fuga does not include the tilt deflection yet, and (ii) its closure model varies eddy viscosity linearly with height, showing over-diffusion of the wake in the vertical direction. A Gaussian filter is applied to the Fuga predictions to get some turbulence mixing and to avoid spurious results from non-sampled Fourier modes in the high frequency side. Since the turbulence closure in Fuga only considers vertical shear effects, and the Gaussian filter is potentially too strong, it can lead to a wake which diffuses too quickly (Ott and Nielsen, 2014). Improvements in these aspects are being implemented but are not ready for this publication. PyWake, similar to the DWM model, also predicts a steady upward projection of the wake, which is not visible in the LES reference data. Finally, there is potentially some coupling between the lateral and vertical deflection of the wake since the wake-averaged veer increases with height. Only the models that model veer can capture this.

To finish the discussion of the mean results for the single wake case, Fig. 6 shows the normalized power error of a hypothetical waked turbine for different downstream positions compared to LES. The aggregated errors are also shown on the side and at the bottom of the plot. As a general trend for all models, the error decreases with downstream distance for the aligned case, higher yaw angles lead to larger error, and for the yawed cases, the power tends to be underpredicted in the far wake. For the aligned case, the RANS model shows small far-wake errors of ∼10 %, whereas underpredicted wake deflection in the misaligned cases increases the far-wake relative error to ∼25 %, yielding an overall error across all cases of about 20 %. Of all the models, the DWM shows the smallest total error, just shy of 10 %, and the maximum error hovers around 15 %. Fuga shows errors of up to 300 % in the near wake (note the y-axis scaling; see the figure caption for more explanation), where non-linear effects dominate and are not possible to be captured by the linear nature of the model; for the aligned case, the far-wake prediction error is below 10 %, but for the non-aligned cases, the error in the far wake maxes out at about 30 %. Finally, the results from PyWake are off by up to 150 % in the near wake, but from 4D onwards, the error decreases to a maximum of around 25 % for all cases, even without modeling the veer.

Figure 6 Relative power error $\mathit{ϵ}=(P-P_{\text{LES}})/P_{\text{LES}}$ compared to LES for the first turbine as well as for several ghost turbines placed at different downstream locations in the wake. For the first turbine, results for a real turbine with coupling to the inflow are shown using untextured bars (indicted by an empty square in the legend) for the HAWC2 turbine in the LES/DWM simulations and a plus symbol for the corresponding RANS LUT result; these data are only available for the first row, as the remaining turbines are ghost turbines. Consequently, for the first turbine, the LES reference is a HAWC2 turbine, whereas HAWC2 ghost turbines are used as the LES reference for the waked turbines. Results for all yaw angles and models are shown. The aggregated error is given by ${\mathit{ϵ}}_{\text{RMS}}=\sum \sqrt{{\mathit{ϵ}}^{\mathrm{2}}}$. A symlog scale is used on the y axis such that values with $\left|\mathit{ϵ}\right|<\mathrm{0.5}$ are shown on a linear scale and larger values on a logarithmic scale.

The instantaneous wake and turbine loading properties will be presented for the LES and the DWM predictions since these are the only two models that predict dynamic results and not just mean properties.

A wake-tracking algorithm was applied to instantaneous flow slices at various downstream locations to dynamically compare the wake behavior of the turbine. Figure 7 presents the power spectral density (PSD) and time series of the wake centerline displacement for the aligned case at 5 rotor diameters downstream. As indicated in the figure labels, a rolling averaging filter with width T_avg is also applied to reduce the sensitivity of the wake-tracking algorithm to small-scale turbulence. For the dynamic wake meandering (DWM) model, two results are shown: the output from the wake-tracking algorithm applied to the flow field and the actual positions of the wake-defining particles, representing the “true” wake center as defined by the model. The particle-based wake center exhibits less motion than the one inferred from the tracking algorithm. This highlights the ambiguity and sensitivity of defining a wake center, especially when working with instantaneous flow slices.

No scientific consensus on the root cause of wake meandering has been reached in the literature, according to Yang and Sotiropoulos (2019). Several plausible mechanisms have been proposed: (i) advection by large-scale inflow structures – the core assumption behind the DWM model; (ii) shear layer instability, akin to vortex shedding seen in flow past bluff bodies; and (iii) wake-to-wake interactions, when multiple turbines are present. In the PSD plot, the cutoff frequency f_C identified by Larsen et al. (2008) and Larsen and Lio (2025) – above which vortices are considered too small to influence wake meandering – is indicated. Also shown is the Strouhal number range $\text{St}=\frac{fD}{U_{\mathrm{hub}}}\in [\mathrm{0.2},\mathrm{0.4}]$, where shear layer instabilities and wake interactions have been reported in the literature.

Compared to the LES reference, the DWM model shows good agreement in the lateral direction, where the dominant frequencies of wake center meandering align well between the two models. In contrast, in the vertical direction, the correspondence in dominant frequencies is less clear. Regarding the magnitude of the spectral peaks, in the lateral direction the DWM overpredicts the energy below the outlined Strouhal numbers and underpredicts it at higher frequencies. In the vertical direction, however, the intensity of the wake center motion is underpredicted across the entire spectrum. The discrepancies in the vertical direction can largely be attributed to the absence of a ground model in the DWM as well as differences in wake deflection. In the lateral direction, the deviations are likely also influenced by the choice of filter used to extract the meandering motion. Other sources have noted that the uniform spatial averaging filter tends to overpredict the energy at certain frequencies while underpredicting it at others (Jonkman and Shaler, 2021). One additional feature that stands out in the frequency plots is the higher energy at very low frequencies in the LES compared to the DWM. The procedure used to compute the spectra from the time series consists of three steps: (i) subtraction of the mean, (ii) application of windowing to account for the non-periodicity of the signal, and (iii) calculation of the Fourier coefficients. These steps ensure that the observed low-frequency energy is not a numerical artifact. Examining the time series shown in Fig. 7, the dashed lines – representing rolling averages over a 5 min window – indicate that although the two signals are correlated, there exists a non-constant offset between them. This offset corresponds to the increased energy at very low frequencies observed in the spectra. Such drift is likely a consequence of differences in wake shape and wake deflection (and its interaction with veer) between the two models. The correlation coefficients between the wake center location time series of the two models are $r_{\text{Y,max}}(\mathrm{\Delta }T=-\mathrm{39}\text{s})=\mathrm{0.80}$ and $r_{\text{Z,max}}(\mathrm{\Delta }T=-\mathrm{35}\text{s})=\mathrm{0.69}$, where ΔT denotes the time offset that maximizes the correlation between the two signals.

Figure 7 Time series and power spectral density (PSD) of the time series of the wake centers of the LES and the DWM predictions at 5D downstream for the aligned case. The top-left side shows the lateral, and the top right side shows the vertical wake displacement. The time series for the (middle plot) lateral and the (bottom plot) vertical deflection are plotted. For the DWM, the PSD of the wake-tracking algorithm applied to the output flow slices and the actual wake particle position are shown. A rolling filter with a width of 20 time steps is applied for the averaged results.

Finally, Fig. 8 shows the relative error, with respect to the LES results, of the damage equivalent loads (DELs) computed for the selected moments over the entire simulation period (i) for a real turbine fully integrated into the flow simulation at $x/D=\mathrm{0}$, and (ii) for ghost turbines placed in the wake of the freestanding turbine. Again, the results for the ghost turbines are standalone HAWC2 simulations supplied with a time series of inflow planes from either LES or DWM simulations that use the induction model in HAWC2 to account for the effect of the turbine on the flow. The DELs were computed for the blade root flapwise moments, the tower bottom foreaft moments, and the tower-top yaw bending moments. Those load channels were selected based on comparison with the literature (Hanssen-Bauer et al., 2023, 2025; Branlard et al., 2024). The errors for the DELs for the blade root bending moment reach up to around 20 %, with larger errors for smaller yaw angles or the fully aligned case. For the tower bottom foreaft DELs, the error also increases with smaller yaw angles and is largest for the aligned case, with a maximum relative error of around 60 %. The predictions for yaw bending moment at the top of the tower are very accurate and increase with increasing yaw angles, but the errors remain below around 10 %.

Figure 8 Relative prediction error of the turbine DELs evaluated using the DWM as compared to the LES. The results are shown for (left) flapwise blade loads at the root, (middle) foreaft tower loads at the bottom of the tower, and (right) yaw bending loads at the top of the tower. The errors are also evaluated for all yaw angles from top to bottom; see the plot labeling. At $x/D=\mathrm{0}$, a real DWM turbine is compared to a real LES turbine, and at $x/D>\mathrm{0}$, DWM ghost turbines are compared to LES ghost turbines.

While the accuracy for the DWM for the prediction of the blade root moments and the yaw bending moment is satisfactory, the tower bottom moment predictions are very far off and warrant further investigation. In literature, similar comparisons between LES and DWM simulation results for the tower bottom foreaft DEL show the same trend. Branlard et al. (2024) plot the ratio of the tower DELs between the second and first row turbine for two IEA 15 MW turbines placed 7D apart in stable and neutral conditions. For the stable case, the ratio between the DELs for the LES is similar to the ones presented for the LES here. In the paper, the DWM underpredicts the ratio by about 20 %, whereas in the case presented here, the DWM underpredicts the ratio by about 40 %. Likewise, Hanssen-Bauer et al. (2023) and Hanssen-Bauer et al. (2025) analyze a row of four NREL 5 MW turbines spaced 7.5D apart, both below- and above-rated conditions. In below-rated conditions with a turbulence intensity (TI) of around 4.6 %, all DWM models underpredict the tower bottom DELs of the waked turbines by roughly 50 %. In above-rated conditions with a TI of 5 %, the underprediction is about 30 % for the second-row turbine and increases up to nearly 100 % for the third and fourth rows. For the below-rated case, additional turbulence intensities of 9 % and 12 % were tested: at 9 % TI, the tower DELs were underpredicted by roughly 10 %, while at 12 % TI, they were predicted with reasonable accuracy. In this paper, the tower bottom bending DEL for the second-row turbine at 7.5D is underpredicted by about 50 %. Although the prediction error for the tower bottom foreaft DEL is larger than the reference cases in literature, it is plausible that the increase in error is due to the increase in the size of the turbine. Underprediction of the DEL has also been observed by Bernard et al. (2024) when comparing different DWM implementations against measurements for a large 6MW offshore wind turbine: a 7 %–35 % underprediction for tower tilt moments and up to 23 % underprediction for tower torsion, which are different but related load channels. For a more detailed causal analysis of the prediction errors, see the discussion in Sect. 4.

3.3 Setup II: two aligned turbines

The results for the second setup with two turbines spaced 5D apart are shown in this section. This setup has been added to see how the ghost turbines compare with an actual turbine for the LES and to see if the models can capture the secondary steering.

Figure 9 shows the mean flow fields for the aligned and steered cases with a yaw offset of $\mathit{\gamma }=-\mathrm{30}\mathit{°}$. Qualitatively, all models capture the increased wake deficit for the second turbine and the wake asymmetry caused by partial overlap with the upstream wake. However, the DWM and PyWake predictions show no additional deflection of the second-turbine wake from secondary steering, because they do not model the lateral velocity deficit that drives the steering. While Fuga models the lateral velocity deficit due to the yaw misalignment, since it already underpredicts the wake deflection for the first turbine, it also underpredicts the deflection for the second one. In the DWM model, the Ainslie velocity deficit was initially projected using the rotor yaw and tilt misalignment, which produced some secondary steering but strongly underpredicted the streamwise deficits. Projecting the deficit only according to the local flow deflection angle was also tested; this yields a minimal lateral velocity deficit and therefore insufficient steering of the downstream wake. Appendix B presents the deflection of the time-averaged wake from the wake-tracking algorithm (Fig. B1) and compares two modified PyWake setups against LES data to further explore secondary steering (Fig. B2). The latter includes two additional cases: (i) shifting the upstream turbine downward to overlap with the wake that would occur under yaw, and (ii) yawing the second turbine according to the rotor-averaged veer observed in the LES. Finally, for wake properties unrelated to deflection, larger inter-model variations appear in both the depth of the velocity deficit and the wake expansion of the second turbine. For more details, see Appendix C where explicit velocity deficits are plotted. Overall, the models tend to overpredict the deficit magnitude to varying degrees, and the predicted wake expansion shows a wider spread than in the single-turbine case. A fuller analysis of these differences is beyond the scope of the present study, as it would require addressing wake superposition, turbulence build-up, and their combined influence on wake expansion.

Figure 9 Time-averaged streamwise velocity at hub height for all models for (left) the aligned and (right) the yawed case with $\mathit{\gamma }=-\mathrm{30}\mathit{°}$.

Figure 10 shows the relative error in the predicted mean power for three cases: the first turbine, the waked second turbine at x=5D, and a ghost turbine at x=10D. In this setup, results are presented for both real and ghost turbines.

For the first-row ghost turbine with LES inflow planes, the errors are about 10 % for the aligned case and 3 % for the yawed case. This indicates that the induction model in HAWC2 does not fully capture the inflow complexity. For the second-row ghost turbine, the aligned case error decreases to 5 %, while the yawed case error increases to 20 %. This increase suggests that the induction model becomes less accurate for more complex inflows, which is expected. When comparing the real second-row turbine instead of the ghost turbine, the differences are smaller for the aligned case. Hence, the conclusions from the previous section remain valid for the other models as well. In the yawed case, all models previously overpredicted the power compared to LES. However, the errors are smaller relative to the real turbine, meaning that the apparent overprediction against the ghost turbine was partly an artifact of the ghost representation. For the third-row turbine, compared again against an LES ghost turbine, all models underpredict the power, likely partially due to not capturing the secondary steering. Yet, as seen for the second row, the lower-fidelity models may in fact be closer to the “true” LES behavior than suggested by the ghost turbine comparison alone. Finally, the total root mean square (RMS) error (shown in the right and bottom panels of Fig. 10) is smaller than in the ghost-only case. This is mainly because the power underprediction for the second turbine is less severe when compared to the real turbine, and because the dataset now contains fewer yawed cases, where all models generally perform worse.

Figure 10 Relative power error compared to LES for the first two turbines of the setup and a ghost turbine placed at 10D. Results for two yaw angles and all models are presented. The dotted bars are HAWC2 ghost turbines, and the bars with circles are obtained from rotor-averaged values plus the power curve.

In the previous section, three DEL channels in the DWM were compared to LESs, where the tower bottom flapwise bending moments deviated by up to 50 % at 5 diameters downstream. Figure 11 now shows the DEL prediction errors for the second setup with two turbines, again including both real and ghost turbine results from the LESs. In contrast to the power predictions, where ghost turbines tended to overpredict relative to the real turbine, the loads show the opposite trend: ghost turbines underpredict the DELs compared to the real turbine. For the blade root flapwise DEL, the prediction errors remain in the same order of magnitude, whether the LES reference is the ghost or the real turbine. For the tower bottom moments, however, the DWM error increases to roughly 75 %, since the LES ghost turbine itself deviates by about 25 % from the LES real turbine for this channel. For the tower-top yaw DEL, the sign of the error in the DWM remains unchanged regardless of whether the ghost or the real turbine is used as a reference. Overall, because loads are more complex than power predictions, there is no simple one-to-one relationship between the DWM accuracy relative to ghost and real turbines in the LES.

Figure 11 Relative prediction error of the DELs for both ghost turbines with LES inflow planes and the DWM as compared to real turbines in the LES. The results are shown for (left) flapwise blade loads at the root, (middle) foreaft tower loads at the bottom of the tower, and (right) yaw bending loads at the top of the tower. The errors are evaluated for two yaw angles; see the labeling of the plot.

Finally, Fig. 12 shows the load spectra of the second turbine for the aligned case. Qualitatively, the spectra look similar between LES and DWM, with the highest peak at 1P and 3P for blade/yaw and tower loads, respectively. However, the peak magnitude is underpredicted in the DWM as compared to LES. The underprediction is greater for higher frequencies, i.e., 2P, 3P, and 6P. Since the tower bottom loads have the highest energy at 3P, the underprediction of the signal energy at higher frequencies in the DWM is consistent with the underprediction of the DEL for this channel.

Figure 12 Load spectra for the waked turbine at 5D downstream of the upstream turbine for different yaw misalignment angles. The DELs are for (left) flapwise blade loads at the root, (middle) foreaft tower loads at the bottom of the tower, and (right) yaw bending loads at the top of the tower.

Table 3 Approximate computation time for a single flow case with a total simulation period of 1.5 h including a 0.5 h spin-up period for the unsteady models, and just one simulation for the steady-state models. “n/a” means not applicable for this specific model. Where appropriate, the number of nodes and cores is also indicated as $n_{\text{nodes}}\times n_{\text{cores}}\times \mathrm{\Delta }{t}_{\text{time per core}}$. Fuga 1 and 2 turbine simulations of the full three-dimensional field include the one-time LUT loading time of 11 s. PyWake simulation time includes only the execution time to obtain power at the turbine locations.

3.4 Comparison of computation time

In addition to the previous comparison between the models results, in the following, the computational cost for each model and setup is described. Table 3 lists the CPU time of all the different models. Please note that these simulations were all run on the Sophia cluster (DTU Computing Center, 2026) which uses first-generation AMD EPY2 7351 processors released in 2017; a more modern cluster setup would provide further speedup. Hence, it is recommended using the table for the intercomparison of models only. Further, the DWM implementation has not been fully optimized for speed and adjustments to the particle advection schemes are expected to provide a decrease in computation time. Finally, for the HAWC2 model of the IEA 22 MW turbine, no optimization has been carried out on the discretization of the structure and the BEM model; the authors expect that with some optimization the model should be able to run close to real time. Despite all of this, the order of magnitude difference between the different fidelities is apparent. Further, the duration of one ghost turbine simulation is much lower than the duration of one LES, and several of them can be run in parallel with different turbine position. This clarifies why using ghost turbine simulations is a useful tool, especially for generating large amounts of data for training load surrogates as for example in Hulsman et al. (2020). Even for the DWM, it is useful as several layouts can be tested without running several longer full simulations.

4 Discussion

The comparison between models of different fidelity revealed consistent trends but also systematic discrepancies that point to model limitations.

The RANS model qualitatively captured asymmetric wake effects due to veer and yawed inflow, lateral velocity deficits, and secondary steering, but the magnitude of these effects was underpredicted. This led to growing errors in power prediction with increasing yaw angles, and an underprediction of wake deflection due to veer, even for the aligned case.

The DWM gave the most accurate power predictions overall, mainly due to the good performance of the wake deflection model, but it did not capture wake asymmetry. Since the inflow conditions tested here only produced moderate veer and curling, the accuracy is expected to decrease under more strongly sheared or curled conditions. For turbine loading, blade and yaw bending moments were predicted with reasonable accuracy for both waked and unwaked turbines. Tower bottom loads were satisfactory for unwaked turbines, but for waked turbines they showed large errors (up to about 75 % in DEL for the second setup). These discrepancies are likely multifactorial, including but not limited to the underprediction of the rotor-averaged TI, underprediction of the wake meandering, neglection of non-axisymmetric wake effects in the Ainslie model, and other simplifying assumptions in the DWM turbulence treatment.

Several papers have investigated the sensitivity of the tower bottom bending DELs to properties of the flow (Dimitrov et al., 2018; Shaler et al., 2023; Doubrawa et al., 2023). Dimitrov et al. (2018) calculated site-specific Sobol indices for unwaked turbines using Mann boxes, finding the strongest sensitivity of the tower bottom bending DEL to mean wind speed and streamwise turbulence along with less-pronounced sensitivity to wind shear and integral length scale of the turbulence box. Shaler et al. (2023) repeated a similar analysis using a DWM model with a Kaimal turbulence box for both waked and unwaked turbulence. They found for ambient turbulence in primary wind direction and shear as the most load-driving parameters, followed by integral turbulence scale and other parameters in the coherence model of the Kaimal generator. They found that for both waked and unwaked turbines, streamwise ambient turbulence and shear were the most load-driving parameters, followed by integral turbulence scale and other parameters in the coherence model of the Kaimal generator. For waked turbines of tertiary importance, accounting for up to 5 % of the load significant events were related to non-streamwise turbulence parameters, wind direction and wake calibration parameters. Rinker et al. (2021) compared the prediction of different DEL channels for a waked turbine using both LES data and constrained turbulence boxes on the basis of the LES data. They showed that the shaft torsion and tower bottom DELs were much more sensitive to the non-stationarity of the inflow than the blade-related load channels. Bernard et al. (2024) showed that introducing some asymmetry to the mean wake deficit in the DWM can lead to some increase of the loads at 3P frequencies for related load channels, which aligns with the underprediction of the tower bottom loads at 3P in Fig. 12. To summarize, even with help from literature, pinpointing a specific cause as the main driver of the discrepancy is not possible. Hence, in the following, suggestions for each of the potential contributions are given, and it is recommended to look into this issue further in future work.

Using a uniform spatial averaging filter for meandering may contribute to the underprediction of the wake meandering unrelated to shortcomings in the model hypothesis. Comparison of the PSD of the wake center line (with and without the filter enabled) reveal that this filter tends to overdampen for certain frequencies. Hence, a more advanced filter could improve the performance. For example, in FastFarm, a non-uniform spatial averaging filter based on the Jinc function is implemented (Jonkman and Shaler, 2021).

For the wake asymmetry, alternative modeling approaches have been proposed in the literature. Branlard et al. (2022) introduced the curled wake model in the DWM framework implemented in FastFarm, which can predict not only lateral velocity deficits but also reproduce some of the curling in the wake. However, for this specific case, a preliminary in-house implementation of the curled wake model reproduced the curling of the wake from about 5D onwards, whereas in the LES dataset, the curling is only visible up to about 5D for the single-turbine case and a yaw misalignment of $\mathit{\gamma }=-\mathrm{30}\mathit{°}$. Further, the curled wake model does not capture the effect of veer on the wake shape, which is more important, since if there is no-wake deflection, accurately predicting that the wake shape becomes more important. More broadly, fast mid-fidelity wake models that capture the combined effects of veer, shear, and yaw on both power and loads remain an important development area. Potential candidates include fast 3D parabolic solvers (Mittal et al., 2017), simple wake deformation models (Abkar et al., 2018), and data-driven surrogates trained on high-fidelity data (Andersen and Murcia Leon, 2023; Schøler et al., 2024).

Additionally, even though the DWM implementation in Dynamiks uses the updated WAT model from Branlard et al. (2024), the TI in the wake is missing key features: (i) the effect of shear on the TI is not modeled, hence the TI is correct at hub height but under- or overpredicted at the top and bottom of the wake, respectively; and (ii) the accumulation of the TI in the farm is not correctly captured. Keck et al. (2015) have proposed improvements for both these problems by connecting the wake TI to the Reynolds stress in the wake and shear in the inflow, but they are not implemented in Dynamiks yet. If a non-axisymmetric deficit model is implemented as suggested in the previous paragraph, then the first problem would also be much easier to solve.

Among the steady-state engineering models, the super-Gaussian wake model in PyWake using the Jiménez model for deflection gave similar accuracy to RANS in the far wake. The inclusion of veer-induced deflection similar to the DWM and case-specific retuning of the near-wake region could potentially improve the results further. Fuga performed similarly to PyWake in the far wake but showed errors up to 300 % in the near wake and overall slightly larger errors than PyWake for all cases, despite modeling more of the physics. Including tilt deflection and potentially modifying the Gaussian smoothing applied to the wake by including a larger number of Fourier modes in the Fourier space solution could improve the results.

Finally, the DWM and PyWake do not model lateral velocity deficits and thus cannot capture secondary steering. This leads to flawed results for yaw setpoint optimization. A straightforward improvement could be to include a lateral velocity deficit term following Zong and Porté-Agel (2020).

5 Conclusions

This study compared power and load predictions for wake steering on a solitary IEA 22 MW turbine and a two-turbine row with 5D spacing, using models of varying fidelity under neutral ABL conditions. Large-eddy simulations were used as the reference since they were the most high-fidelity model out of the tested models and provided a consistent inflow for model intercomparison. Overall, all models reproduced the qualitative power and load variation trends with yaw angle and downstream position. Quantitatively, however, substantial differences in accuracy with respect to large-eddy simulations were observed. The DWM provided the most accurate power predictions but heavily underestimated tower bottom loads for various reasons. RANS captured asymmetric wake effects but underpredicted their magnitude, leading to errors in power prediction. Steady-state engineering models (PyWake and Fuga) produced reasonable predictions in the far wake for all yaw angles but could benefit from small model-specific improvements. PyWake could not reproduce secondary steering since the lateral velocity deficits were not modeled as opposed to Fuga where lateral and vertical components of the deficit are modeled. Several suggestions for further improvements to the models have also been discussed. The improvements focus on capturing wake asymmetry due to non-uniform inflow, mean wake deflection, and lateral velocity deficits, as well as better capturing the main physics behind wake meandering.

The main takeaway is that while the presented models can capture broad trends for the neutral boundary layer case presented here, their quantitative accuracy is limited by the aforementioned reasons. For practical applications, this means that optimization strategies relying on the tested models should be applied cautiously, especially for yaw offsets larger than 10°. Further, the results may become less reliable in cases with stronger veer.

Appendix A: More detailed TI plots for the DWM results

Figure A1 (Left) streamwise turbulence intensity for the DWM and the LESs at $x/D=\mathrm{5}$ and (right) streamwise evolution of the rotor-averaged TI for both models.

Appendix B: Deflection plots for the steered setup

Figure B1 Wake center position for the two-turbine case in (left) horizontal and (right) vertical direction as obtained from the wake-tracking algorithm of the time-averaged flow field for all models and two yaw angles.

Figure B2 Relative wake deflection for the steered two-turbine setup using PyWake and LES. The plotted deflection is relative to the deflection for the aligned case $y=y_{\text{steered}}-y_{\text{aligned}}$ to account for the influence of veer, which is only present in the LES. For PyWake, also second option is shown where the second turbine is yawed in accordance with the rotor-averaged yaw angles obtained from the LES data.

Appendix C: Deficit profiles

Figure C1 Normalized velocity deficit $\mathrm{\Delta }U=\frac{U_{\text{hub}}-U}{U_{\text{hub}}}$ at hub height for the aligned case for (left) a single turbine and (right) two turbines in a row with spacing Δx=5D.