Wake turbulence modeling in stratified atmospheric flows using a novel $k−ℓ$ model Open Access

Atmospheric boundary layer (ABL) flow environments with wind shear and stratification can pose challenges for fast-running engineering wake models.^1–4 Error and uncertainty in wake modeling can impact predictions of wind farm power production used for both wind farm design⁵ and wind farm flow control.⁶ Many state-of-practice engineering wake models were derived under assumptions of uniform or neutral inflow.^7–9 These assumptions greatly simplify wake physics by neglecting or parameterizing ABL effects such as atmospheric stratification and the impact of wind shear on wake turbulence. Yet wind turbine wakes and wind farm performance strongly depend on atmospheric stability.^10–14 

Recent studies have developed extensions of existing engineering models that account for the effects of wind direction shear (veer) on wind turbine wake advection¹⁵ or through a shape function that depends on azimuthal angle with superimposed veer.^3,16 These extended models show substantially improved performance relative to standard axisymmetric analytical wake models.^3,15,16 However, by relying on an analytical model form as the starting point (i.e., a Gaussian wake model), these models often neglect the nonlinear interactions between the wake momentum and the turbulence added by the wake itself. The turbulence generated by the presence of a wake is often called wake-added turbulence.^17,18 While turbulence can have a significant impact on farm power and farm efficiency in analytical wake models,¹⁹ turbulence is often parameterized in a simple manner that does not allow for feedback between the wake deficit and the wake-added turbulence. In many analytical wake models, the wake-added turbulence is only relevant in the (often linear) wake spreading parameter for waked turbines.²⁰ This is because the wake recovery is generally assumed to only depend on the turbulence intensity incident to the turbine,²¹ rather than the turbulence in the wake. In other words, the wake of a wind turbine operating in freestream ABL flow does not account for the effect of wake-added turbulence on wake recovery. In contrast, in the balance of mean momentum in the wake, the wake recovery is controlled by the spatially varying turbulence, which is the combination of inflow ABL turbulence and wake-added turbulence.⁴ Recognizing this, several recent engineering wake models have incorporated the effect of wake-added turbulence on the wake recovery.^22–24 Yet, a challenge in fast-running wake models is how to incorporate the effect of stratification on the wake-added turbulence, and therefore, the spatially varying wake recovery.

Wake-added turbulence is strongly dependent on atmospheric stratification due to the effects of atmospheric stability on ambient ABL turbulence intensity, wind speed shear, and wind direction shear.⁴ In addition, the buoyant destruction (generation) of turbulence kinetic energy in stable (unstable) ABL flows plays a secondary role.⁴ As discussed earlier, existing analytical wake models often do not directly account for these effects due to their simplified, analytical form. Instead, to improve wake turbulence and momentum modeling in the stratified ABL, we leverage an alternative engineering modeling approach based on the parabolized Reynolds-averaged Navier–Stokes (RANS) equations, such as the curled wake model.^25,26 The curled wake model solves a simplified and parabolized form of the RANS equations to calculate the velocity deficit along with a Boussinesq eddy viscosity model for the turbulence closure. Working within this framework allows for improvement to the physical representation of the flow field through the incorporation of ABL processes and forcings, while still maintaining computational efficiency. The primary benefit of this approach is that the mean wake momentum deficit and wake-added turbulence are naturally coupled in the dynamics. However, a primary gap is in accurately modeling wake-added turbulence in stratified ABL flows.⁴ Previous works that employ the curled wake model have used simple zero-equation mixing length modeling^25,26 and data-driven turbulence modeling.²⁷ In the present work, we utilize the parabolized RANS curled wake model framework, and we develop a novel, physics-based wake-added turbulence model that is applicable to wind turbine wake modeling in a wide range of stably stratified ABL flows based on physical insights from prior work.⁴ 

The remainder of this work is organized as follows: the methods are presented in Sec. II with a description of the extended curled wake model developed in this study (Sec. II A), the development of a $k−ℓ$ model (Sec. II B), details of the baseline wake models that are used for comparison (Sec. II C), numerical details of the engineering wake model developed here (Sec. II D), and (Sec. II E) details of the large eddy simulation (LES) data used for testing and validation. The results are presented in Sec. III, with validation against the LES data and comparison with the baseline, existing wake models. Finally, conclusions are provided in Sec. IV.

Previous studies that solve similar parabolized RANS equations for wind turbine wakes as Eq. (3) have utilized analytical or zero-equation eddy viscosity models for the Reynolds stresses.^25–27,34 While these models have the advantage of simplicity and computational efficiency, they lack turbulence representation that comes with the additional complexity of traditional RANS models.³⁵ The working hypothesis of this study is that this additional detail in turbulence modeling, especially the three-dimensional spatial information and the incorporation of buoyant forcing and wind shear effects, is especially important for modeling wakes in the stratified ABL.

Both existing analytical wake models^7,8 and the curled wake model²⁶ take initial conditions from momentum theory,³⁶ such that $Δu¯=−2auh$ at $x=x0$ downwind of the turbine, where $uh$ is the hub-height wind speed, a is the axial induction factor $a=12(1−1−CT),$ and $CT$ is the thrust coefficient. Specifically, the wake initial condition for $Δu¯$ is applied at the rotor plane $(x=0)$ and the eddy viscosity form in Eq. (6) yields relatively less recovery until $x=x0$⁠, where wake-added turbulence then contributes to the eddy viscosity. The Heaviside function limits turbulence in the near-wake region as discussed earlier. Wake-added TKE becomes active in the far wake (⁠ $x>x0$⁠) and is modeled via a transport equation. We target the modeling of $kwake$⁠, rather than the full TKE, as it has been shown that the wake-added turbulence is highly dependent upon stability and by extension speed and direction shear.⁴ 

In Eqs. (6) and (8), the mixing length $ℓ$ is taken to be the wake width $yD$⁠, which is defined as the distance from the wake centerline at which the wake has recovered to 5% of the ambient. The wake width is a common choice for the mixing length $ℓ$ in wind turbine wake models.^37–40 In RANS modeling of the stratified ABL, the mixing length is often influenced by a variety of factors that limit mixing, including the wall normal distance and stability.⁴¹ Previously, Klemmer and Howland⁴ showed that stability has a strong impact on wake-added TKE. Therefore, the primary goal of this study is to develop an improved model for wake-added TKE and to incorporate the wake-added TKE model into a wake model for mean momentum. Since aspects of stability are encoded in the wake-added TKE, we select the wake width as a simple model of the mixing length. The mixing length approximation using the wake width is quantitatively evaluated in  Appendix B, where the results show that the approximation yields reasonable accuracy for neutral and moderate stabilities, but has increased error at higher stability. The model constants $Ck1$⁠, and $Ck2$ are set to 1, and $Cν$ is set to 0.04. To arrive at these values, we performed a sensitivity analysis across all stability cases with aerodynamic roughness lengths of $z0=1$ cm and $z0=50$ cm (further details on the LES simulation data are in Sec. II E). Details of this parameter sensitivity analysis are given in  Appendix A.

To assess the performance of the proposed $k−ℓ$ model, we compare this model against the prediction from three other models: the analytical Gaussian wake model,^8,9 the analytical Jensen wake model,⁷ and a turbulence model from Scott et al.²⁷ that uses the curled wake model. Importantly, the parameters used for these models are taken from the literature, meaning that unlike the $k−ℓ$ model presented in this work, the comparative models are not tuned for the sweep of ABL cases for $z0=1$ cm and $z0=50$ cm. Despite this, all models are evaluated on unseen data (⁠ $z0=10$ cm). In other words, the new $k−ℓ$ model has not been calibrated to the LES data for $z0=10$ cm, which is the data to which we make comparisons in this study.

We compare the proposed $k−ℓ$ model to two existing, analytical wake models: the Gaussian wake model^8,9 and the Jensen wake model.⁷ For the evaluation of both analytical wake models, we use implementation within the open-source FLOw Redirection and Induction in Steady State (FLORIS) tool version 4.1,⁴² developed by the National Renewable Energy Laboratory (NREL).

For all models, the coefficient of thrust $CT$ is set to 0.75 to match the LES data. For the $k−ℓ$ and Scott models, which used the curled wake model framework, the base flow is taken from LES. At $x/D=0$⁠, both the $k−ℓ$ and Scott models are given an initial condition of $Δu¯=−2auh$ in the rotor area (following the example in Martínez-Tossas et al.²⁵) Additionally, for the $k−ℓ$ model, $kwake$ is given an initial condition of 0, uniform in space, at $x/D=0$⁠. A no slip condition is applied at the bottom boundary, and a stress-free condition at the top wall.

The near-wake length $x0$ is a complex function of the turbine control and the incident wind conditions.^9,33 There is a limited understanding of the near-wake length in sheared and stratified ABL flow conditions that we investigate in this study. To isolate the evaluation of the modeling of the turbulent portion of the wind turbine, we provide both the baseline engineering models and the proposed parabolized RANS $k−ℓ$ model with the near-wake length directly measured from LES as the position in the wake where $Δu¯>0$ for the first time.³³ This is the point at which wake recovery begins. This is defined by the location at which the streamtube-averaged wake deficit is at its minimum value. We leverage these LES measured near-wake lengths $x0$ rather than using an engineering model for $x0$ (e.g., Liew et al.³³) because the existing models do not account for stability. Future work is suggested to improve modeling of $x0$ in stratified ABL flow conditions.

We use large eddy simulation (LES) of modeled wind turbines immersed in stratified ABL flows for validation and performance assessment for the model. All data are generated using the PadéOps LES code (https://github.com/Howland-Lab/PadeOps),^32,46 which is an open-source, pseudo-spectral computational fluid dynamics solver. The solver uses Fourier collocation in the horizontal directions and a sixth-order staggered compact finite difference scheme in the vertical direction.⁴⁷ For the temporal integration, a fourth-order strong stability-preserving variant of a Runge–Kutta scheme is used.⁴⁸ The sigma subfilter-scale model is used⁴⁹ for the subgrid stresses with a turbulent Prandtl number of 0.4 for the scalar diffusivity.

All simulations employ a fringe region in the streamwise and lateral directions to force the inflow to the desired profile.⁵⁰ The inflow is specified via the concurrent precursor method,⁵¹ in which two simulations are run concurrently: a primary simulation with the turbine and a precursor simulation of the same domain as the primary but without the turbine. This ensures that the wind turbine encounters physically consistent stratified ABL turbulence as its inflow condition. Each LES simulation has a single turbine that is modeled with the actuator disk model without rotation⁵² and has a diameter of $D=126$ m, a hub height of 90 m, and a thrust coefficient of $CT=0.75$⁠. The actuator disk model is used here to evaluate engineering wake models that are also built on actuator disk modeling assumptions.^9,26,33 For first validation of the wake model and for computational tractability for the 15 different atmospheric conditions investigated here, we only consider the modeling of a single turbine wake within the ABL. Future work, outside the scope of this paper, should evaluate the modeling of multiple wind turbine wakes within wind farms operating in the stratified ABL.

The model developed in this study is calibrated and tested against 15 ABL cases that correspond to a sweep over five surface cooling rates $∂θ∂tz=0=[0,−0.25,−0.5,−0.75,−1]$ K h⁻¹ and three surface roughness lengths $z0=[1,10,50]$ cm. In total, there are three conventionally neutral boundary layer (CNBL) cases and 12 stable boundary layer (SBL) cases. More specifically, $z0=[1,50]$ cm is used for calibration and $z0=10$ cm is used for testing. Thus, $z0=10$ cm is not seen by any model for calibration. In this study, we focus on neutral and stable ABLs because wake losses are typically less important in unstable ABL conditions,^12,13 but future work may consider unstable ABL conditions. Data for all flows with $z0=10$ cm come from Klemmer and Howland.⁴ The additional ten cases with $z0=[1,50]$ cm are generated for this study to provide stratified flows with varied freestream turbulence intensity for the same surface cooling rates. The 15 cases examined in this study encompass a range of flow conditions, spanning from neutral through to very stable stratification. The inclusion of turbulence intensity variations (through the surface roughness parameter) within each boundary layer enables a more detailed investigation into the role of inflow turbulence on wake-added turbulence in stratified flows. The SBLs are named according to their surface cooling rates (SBL $−0.25$⁠, SBL $−0.5$⁠, SBL $−0.75$⁠, SBL $−1$⁠), and—unless otherwise specified—results are shown for $z0=10$ cm.

All cases are driven by a geostrophic wind speed of 12 m/s. The grid spacing in the streamwise and lateral directions is 12.5 m, with a streamwise domain length $Lx$ of 4800 m and a lateral domain length $Ly$ of 2400 m. The vertical extent is 2400 m in the CNBLs and 1600 m in the SBLs. All cases have vertical grid spacing of $Δz=6.25$ m. The Coriolis frequency is $1.03×10−4$ s⁻¹, which corresponds to a latitude of 45° N. A constant surface heat flux $w′θ′¯w=0$ K m s⁻¹ is prescribed in the CNBL. In the SBLs, the surface boundary condition is specified as a cooling rate.⁵³ In all cases, the initial potential temperature profile is set to 301 K up to 700 m for the CNBL and 50 m for the SBLs. Above this height is an inversion with strength 0.01 K m⁻¹. Further simulation details are provided in Table I.

In this section, the proposed $k−ℓ$ model is validated through a comparison with LES, and performance is assessed through comparison with the three models detailed in Sec. II. For both the validation and performance assessment, we show qualitative visual results generated using the CNBL, SBL $−0.25$⁠, and SBL $−1$ cases, all with $z0=10$ cm. These cases cover the full range of stabilities, and trends remain consistent across roughness lengths as shown in the aggregated metrics, which are shown for all stabilities and roughness lengths.

We assess the models via four metrics: pointwise error in $Δu¯$⁠, the streamtube-averaged wake deficit $〈Δu¯〉s$⁠, the centerline wake deficit $Δu¯c$⁠, and the estimated power of a hypothetical waked turbine P. Together, these four metrics illustrate the predictive capabilities of the proposed model over a varied range of levels of description pertinent to flow field and power prediction.

Starting with the most descriptive qualitative assessment of the flow field predictions, we compare the models and LES pointwise through slices of $Δu¯$ in the $y−z$ plane. Figures 1–3 show the wake deficit of the CNBL, SBL $−0.25$⁠, and SBL $−1$ at three streamwise locations: 5D, 10D, and 15D downstream of the turbine (and initial condition).

Starting with the CNBL (Fig. 1) at $x/D=5$⁠, all the models underpredict the magnitude of the deficit at this location, meaning that the turbulence, which is responsible for wake recovery, is overpredicted. The $k−ℓ$ model most closely approximates $Δu¯$ at this location, while the Jensen and Gaussian models most severely underpredict the deficit. This can be seen moving downstream to $x/D=10$ and $x/D=15$⁠, where the Gaussian and Jensen models differ most from the LES field, and the $k−ℓ$ and Scott models more closely capture the wake velocity deficit and shape.

This pattern is additionally observed in the two SBLs in Figs. 2 and 3, where again we see that the Scott and Gaussian models significantly underpredict the wake deficit magnitude (corresponding to an overprediction of turbulent diffusion). In contrast, the $k−ℓ$ model more closely approximates the correct magnitude. The Scott and $k−ℓ$ models have a further advantage of being able to capture the skewing of the wake in the stable cases, which is not present in the Gaussian model predictions. The wake skewing is a direct consequence of the advective terms within the prognostic equation for the mean wake momentum deficit in Eq. (3), rather than a heuristic parameterization. While both the Scott and $k−ℓ$ models fail to exactly reproduce the wake shape, the $k−ℓ$ model appears to more faithfully capture the wake skew, likely due to the inclusion of a three-dimensional eddy viscosity.

We reduce the granularity of the comparison and performance metrics by introducing the streamtube-averaged wake deficit $〈Δu¯〉s$⁠. In contrast to the pointwise comparison of $Δu¯$ and the associated error, the streamtube average focuses the analysis on the turbine wake and its recovery. However, in contrast to the pointwise error metric, the streamtube-averaged wake deficit does not directly penalize errors in predictions of the wake shape that were shown in Sec. III A. Figure 5 shows $〈Δu¯〉s$ for the CNBL, SBL $−0.25$⁠, and SBL $−1$ (with $z0=10$ cm). The $k−ℓ$ model tends to be relatively well-matched to the LES across the three ABLs. The analytical models tend to underpredict the magnitude of the streamtube-averaged wake deficit (as observed in Sec. III A).

It should be noted that all the models fail to capture the precise streamwise development of the flow, a result that is seen in all error metrics. This is likely due to the difficulty in capturing the near-wake behavior. None of the models have a mechanism for increasing the magnitude of the wake deficit downwind of the initial condition. This would necessitate more accurate near-wake models, particularly for cases with complex shear.

The centerline velocity deficit $Δu¯c$ contains less spatial information than $〈Δu¯〉s$⁠. This reduced granularity allows for a comparison that does not directly encode the effect of stability on the shape of the wake deficit. The centerline corresponds to the location of hub height in each $y−z$ slice of the wake, assuming the wake centerline is not deflecting in a mean sense given inflow-rotor alignment and limited Coriolis effects from the high Rossby number.²⁸  Figure 7 compares the centerline wake deficit for each model. In the CNBL in Fig. 7(a) and SBL $−0.25$ in Fig. 7(b), the $k−ℓ$ model shows the best agreement with the LES. Interestingly, in the most stable case [SBL $−1$⁠, Fig. 7(c)], the analytical models exhibit the closest agreement with the LES, although this will depend on the selected tuning parameters for the wake spreading rate $kw$ described in Sec. II C.

Figure 10 shows the NMAD for power across the 5  $z0=10$ cm flows. Except for the SBL-0.25 cases, the $k−ℓ$ model outperforms all other models in predicting power. To further test the robustness of the $k−ℓ$ model in its ability to predict the spatial distribution of the wake, we move the hypothetical turbine laterally, such that it is only partially waked. This emphasizes the degree to which the wake models capture the skewing from wind direction shear. The results of this for $εP,x$ are shown in Fig. 11, where the location of $yT$ has been moved from $yh$ to $yh+0.5D$⁠. Overall, the $k−ℓ$ model still performs the best, while the Gaussian and Jensen model results tend to worsen as the Obukhov length decreases. This is because the axisymmetric Jensen and Gaussian wake models do not account for wake skewing. Additionally, in both Figs. 10 and 11, the $k−ℓ$ model consistently outperforms the Scott model, likely due to the spatial information available within the turbulence model.

In Secs. III A–III D, we analyzed four distinct metrics to assess model performance. Here, we aggregate these results to more compactly demonstrate the improved performance of the $k−ℓ$ model over the other three models across the ABL conditions and error metrics considered. In Table II, we show the average NMAD values (⁠ $εs,x$⁠, $εc,x$⁠, and $εP,x$⁠) across the 5  $z0=10$ cm ABL test cases. Across all metrics, the $k−ℓ$ model has the lowest average error. For the hypothetical turbine power production, the $k−ℓ$ model provides a 70% decrease in error relative to the Gaussian model and a 60% decrease relative to the Jensen model. This is important because predicting power for downwind turbines is the typical application of engineering wake models.

We additionally look across all four error metrics, and we quantify the fraction of cases and metrics in which each model provided the lowest, second lowest, second highest, and highest error. These results are shown in Table III, where we see that the $k−ℓ$ model provides the lowest error 80% of the time and the second lowest error 5% of the time. Importantly, the $k−ℓ$ model never has the highest error, which is in contrast to all other models presented here. These results demonstrate the robustness of the proposed fast-running wake model across ABL conditions and error metrics, including both flow field characteristics and power predictions.

The proposed $k−ℓ$ model is a new model for wake-added turbulence, in addition to a new fast-running model for mean streamwise velocity deficit. For the proposed $k−ℓ$ model, we assess the predictions of wake-added turbulence kinetic energy $kwake$⁠. Figure 12 shows slices of $kwake$ in the $y−z$ plane for the CNBL, SBL $−0.25$⁠, and SBL $−1$⁠. The $k−ℓ$ model captures the overall shape of $kwake$⁠, as well as the magnitude. The other models are unable to predict three-dimensional fields for $kwake$⁠, which is an inherent advantage of the proposed $k−ℓ$ model.

Figure 13 shows the maximum value of $kwake$ at each x location in the wake for the $k−ℓ$ and Crespo–Hernández models compared with LES for the CNBL, SBL $−0.25$⁠, and SBL $−1$ starting from the start of the far wake region as defined by the near-wake length. It should be noted that due to the asymptotic nature of the Crespo–Hernández model and the non-monotonic behavior of the near-wake length, which is furthest downstream for SBL $−0.25$ and furthest upstream for SBL $−1$⁠, the maximum value of the model appears to be non-monotonic. However, the model is only changing between the three flows via $I0$⁠, so the behavior is actually monotonic as would be expected. The maximum value for the empirical model visualized in this figure is a function of the near-wake length, which is non-monotonic, rather than the Crespo–Hernández model.

In the CNBL, both models compare well with LES. For the stable cases, the Crespo–Hernández model is unable to predict the changing location of the peak in the maximum $kwake$ profile. While the $k−ℓ$ model does not perfectly capture the magnitude or location of the peak (particularly in SBL $−1$⁠), this model does provide results that compare more favorably with the LES.

As a final note on this comparison between the $k−ℓ$ and Crespo–Hernández models, it is important to emphasize that the current approach in most engineering wake models is to predict the wake expansion rate based on the turbulence intensity incident to the turbine. Therefore, the turbine wake is often not affected by the wake-added turbulence generated by the turbine itself. As a specific example of this, the Gaussian and Jensen model predictions shown throughout this study do not even use a wake-added turbulence model. Several recent engineering wake models have recognized this and incorporated the effect of wake-added turbulence on wake recovery,^22–24 but the influence of stratification on the wake-added turbulence has not been a focus. The $k−ℓ$ model provides dynamically consistent mean velocity predictions, which are based on how the spatially evolving wake-added turbulence, explicitly depending on stability, influences wake recovery.

In this work, a novel wake turbulence model is introduced that utilizes a one-equation $k−ℓ$ model for the turbulent eddy viscosity within the curled wake model framework that solves the parabolized RANS equations for the wake velocity deficit in an efficient manner. It takes $≈3.5$ s on a standard computer to run the $k−ℓ$ model for a single turbine case on the given domain, which has $106$ grid points. A model equation for the wake-added turbulence kinetic energy is developed with a mixing length based on the wake width and near-wake length. The model is validated against LES data for five ABL test cases corresponding to a sweep over the Obukhov length. The model was also compared against existing engineering wake models to assess performance. Four error metrics (pointwise, streamtube, and centerline wake deficits and turbine power production for a downwind, waked turbine) are used to evaluate the models. The proposed $k−ℓ$ model exhibited the best performance across the five ABL conditions investigated that are out-of-sample from the calibration of parameters, with a 65% mean reduction in prediction error for the power production of a downwind turbine. The proposed $k−ℓ$ model also yielded the lowest error 78% of the time across all error metrics and ABL conditions. This illustrates the advantage of this model for wakes in diverse stratified, turbulent, atmospheric boundary layer flow. While the $k−ℓ$ model was tested in wind conditions out-of-sample from the calibration of parameters, the model calibration was able to benefit from data that come from the present LES using different roughness lengths, which could give the model some advantage over wake models from the literature. However, we anticipate that the qualitative conclusions are robust to this, namely, that incorporating a model for wake-added TKE that depends on stability into the wake momentum solution improves predictions. By improving the fidelity of the turbulence model, especially in its ability to capture three-dimensional effects and stratification, we were able to better capture the wake behavior, even in the most extreme of stratified cases (lowest Obukhov length). This improvement comes from the nonlinear interaction between the mean wake momentum deficit and the stability-dependent wake-added turbulence in space. This spatial modeling is in contrast to the common approach in-state-of-practice engineering wake models where the wake recovery only depends on the turbulence intensity incident to the turbine generating the wake without the feedback inherent to the dynamics or where wake recovery depends on a wake-added turbulence model that does not depend on stability.

Future work should investigate several important avenues for improvements to the given model. Most importantly, a near-wake length model should be developed that accounts for stability so that the near-wake length is not only dependent upon ambient turbulence intensity. Furthermore, the proposed steady model could be extended to dynamic wake modeling. The model should be extended to incorporate Coriolis effects²⁸ to capture lower Rossby number regimes associated with larger wind turbines or lower wind speeds. Finally, the model must be extended and evaluated in predictions of wind farm flows, beyond a single wake, where wake superposition, blockage, and farm-level entrainment play important roles in the turbulence and wind farm dynamics. Future investigation should study appropriate mixing length models for wind farm flows. Future work could consider coupling this parabolized RANS based wake model with a top-down wind farm model to jointly capture turbine wakes and wind farm array effects.⁵⁸ 

The authors gratefully acknowledge support from the National Science Foundation (Fluid Dynamics program, Grant No. FD-2226053, Program Manager: Dr. Ronald D. Joslin). All simulations were performed on the Stampede3 supercomputer under the NSF ACCESS project ATM170028.

The authors have no conflicts to disclose.

Kerry S. Klemmer: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal). Michael F. Howland: Conceptualization (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

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

The $k−ℓ$ model has three tunable parameters: $Cν$⁠, $Ck1$⁠, and $Ck2$⁠. As stated in Sec. II B, $Ck1$ and $Ck2$ are set to 1 and $Cν$ is set to 0.04 in the present modeling framework. To arrive at these values, we performed sensitivity analysis in which samples were independently drawn from uniform distributions for all three parameters, and we used the streamtube-averaged velocity deficit as our performance metric. The sampling distributions for $Ck1$ and $Ck2$ are given by $U(0,2)$ and the sampling distribution for $Cν$ is given by $U(0,0.1)$⁠. These distributions assume little prior knowledge, though the support of the $Cν$ distribution is dictated by known values of this parameter in the literature.³⁵ Across all ABL cases, the results were found to be relatively insensitive to changes in $Ck1$ and $Ck2$⁠, so these were set to their mean values. An example of this insensitivity is shown in Fig. 14, which illustrates that in the two cases shown, the error is most strongly affected by changes in $Cν$⁠. For $Cν$⁠, we used the range of values across the $z0=1,50$ cm cases and used the average value that minimized the error across these cases.

To further assess the parameter values, we tested $Cν$ values of 0.03 and 0.05. The results of these perturbed parameter values are given in Table IV for the averaged NMAD values for streamtube-averaged wake deficit, centerline wake deficit, and hypothetical turbine power production and are compared to the baseline models for reference. Overall, we find that the $k−ℓ$ model in all cases has the lowest errors for streamtube-averaged wake deficit and power production, and the second lowest error for centerline wake deficit, with the Jensen model yielding the lowest error in this metric when $Cν=0.05$⁠, but with the $k−ℓ$ model yielding lower error than the Jensen model when $Cν=0.03$⁠. These sensitivity experiments show that the $k−ℓ$ model maintains performance over a range of parameter values.

The wake width is a common choice for the mixing length $ℓ$ in wind turbine wake models.^37–40 This mixing length is an approximation of wakes in stratified ABLs, where the mixing length is often influenced by a variety of factors that limit mixing, including the wall normal distance and stability.⁴¹ To evaluate this mixing length model quantitatively, we compute the mixing length $ℓ$ from LES for three different stabilities as $ℓ2=−(∂Δu¯/∂z· ∂Δu¯/∂z)/u′w′¯wake$⁠. The mixing length $ℓ$ is computed using least squares regression in a domain that is $y∈[−2.5D,2.5D]$ and $z∈[zh−0.6D,zh+D]$⁠, where $zh$ is the hub-height. The results are shown in Fig. 15. The simple and computationally efficient choice of the wake width as the mixing length $ℓ$ provides reasonably accurate predictions in neutral and moderate stabilities. The errors are larger in stronger stability, likely because of more severe effects on the mixing length from the stratification. A comparison of prediction error from the $k−ℓ$ model using various mixing length models, including the wake width used in this study or constant mixing lengths based on the turbine radius, is shown in Fig. 16. The results show that using the wake width as the mixing length has lower error than using a constant value for the mixing length. More complex mixing length models, especially in the context of wind farm flows with multiple turbine wakes, should be investigated in future work.