An analytical formulation for turbulence kinetic energy added by wind turbines based on large-eddy simulation

Wind turbine wakes are plume-like regions characterized by reduced wind speed and enhanced turbulence kinetic energy (TKE) that form downstream of wind turbines. Numerical mesoscale models, like the Weather Research and Forecasting (WRF) model, are generally effective at reproducing the wind speed deficit but lack skills at simulating the TKE added by wind turbines. Here we propose an analytical formulation for added TKE by a wind turbine that reproduces, via least-squares error parameter fitting, the main features of the three-dimensional structure of added TKE as simulated in previous large-eddy simulation (LES) studies, including a streamwise peak at x= 4D–6D (where D is the turbine diameter), a vertical peak near the upper-rotor region, and an annular Gaussian-like distribution along the rotor edge. Validation of the proposed formulation against independent LES results and wind tunnel observations from the literature indicates a promising performance in the case of a single wind turbine wake. The ultimate goal is to insert the proposed formulation, after further improvements, in the WRF model for use within existing or new wind farm parameterizations.

1 Introduction

Over the last 2 decades, there has been a significant surge in renewable energy, especially wind energy, due to its vast potential worldwide (Archer and Jacobson, 2005) and the global shift towards low- or no-carbon sources of energy to fight global climate change. For the first time in history, wind has provided over 10 % of the total electricity production of the world, with a global total installed capacity of approximately 824 GW (International Renewable Energy Agency (IRENA), 2022).

When wind turbines are installed in a wind farm, their layout is a compromise between two opposite needs. On one hand, the distance between turbines should be as small as possible to minimize the total lease area and to reduce cable length costs. On the other hand, their distance should be as large as possible to minimize so-called “wake losses” and maximize total power output. Wake losses are the undesirable consequence of wakes, which are plume-like regions characterized by reduced wind speed and enhanced turbulence kinetic energy (TKE) that form downstream of wind turbines (Stevens and Meneveau, 2017). As a result of the reduced wind speed and the enhanced turbulence in the wakes, not only is the performance of downstream turbines in the same wind farm reduced (Archer et al., 2018), but the power production of neighboring wind farms is also affected (Nygaard, 2014). As such, wake losses are one of the most important issues affecting wind farm production today (Göçmen et al., 2016), and the correct estimation and prediction of wind turbine wakes and their evolution are fundamental to ultimately ensure reliable power production from a wind farm (Ye et al., 2023).

Historically, most analytical wake models were developed to address the issue of wake impacts on power production of downstream turbines; thus they focused on the most relevant parameter, the wind speed deficit, and only a few were developed to address the issue of added turbulence (discussed later). Jensen (1983) proposed the first analytical wake model to predict wind speed deficit behind a wind turbine, based on conservation of momentum. The Jensen model predicts a uniform top-hat wind speed deficit across an expanding rotor area that grows with downstream distance x at a specified constant rate. Since then, a series of analytical wake models have been proposed with different assumptions (Katic et al., 1986; Larsen, 1988; Frandsen et al., 2006), including the model by Barthelmie et al. (2004), who fitted field data with an analytical formula for the hub-height wind speed deficit as an inverse power function of x. The first Gaussian model for the wind speed deficit was developed by Bastankhah and Porté-Agel (2014) and later improved by Xie and Archer (2015), who proposed two different expansion rates along z and y, rather than one for both as in Bastankhah and Porté-Agel (2014). These values were found by fitting high-fidelity simulation results. The Gaussian model was also the foundation of more recent formulations (Gao et al., 2016; Wang et al., 2023). A review of the performance of analytical wake loss models for wind speed deficit caused by wind turbines can be found in Archer et al. (2018).

Only a few analytical wake models provide estimates of the added turbulence by wind turbines, and they all use turbulence intensity (TI) as the variable of interest. Quarton and Ainslie (1990) suggested the first empirical formulation for the maximum added turbulence intensity by a wind turbine (ΔTI_max) by considering the effect of freestream turbulence intensity (TI_∞) and the thrust coefficient of the wind turbine (C_T):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathrm{TI}}_{\max }=\mathrm{4.8}{C}_{\mathrm{T}}^{\mathrm{0.7}}{\mathrm{TI}}_{\mathrm{\infty }}^{\mathrm{0.68}}{\left({\displaystyle \frac{x}{x_N}}\right)}^{-\mathrm{0.57}},\end{array}$$

where x_N is a scaling parameter that represents the length of the near wake. Several modifications to this formulation have been proposed in the literature, for example by Ainslie (1990) and Xie and Archer (2015).

One of the most successful variations of Eq. (1) was that introduced by Crespo and Hernández (1996), who divided the wake into two different regions: the near wake (x<3D, where D is the rotor diameter) and the far wake (x≥3D); they then developed a different equation for each:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathrm{\Delta }{\mathrm{TI}}_{\max }\\ & =\left\{\begin{array}{ll}\mathrm{0.362}(\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}})& x<\mathrm{3}D\\ \mathrm{0.73}{\left({\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}}}{\mathrm{2}}}\right)}^{\mathrm{0.8325}}(T{I}_{\mathrm{\infty }}{)}^{-\mathrm{0.0325}}{\left({\displaystyle \frac{x}{D}}\right)}^{-\mathrm{0.32}}& x\ge \mathrm{3}D.\end{array}\right.\end{array}\end{array}$$

The reason for the two different formulations in the two regions is that, in the near wake, the influence of the rotor aerodynamics, such as blade aerodynamics, stalled flow, and the presence of tip vortices, is predominant on the wake (Vermeer et al., 2003). Tip vortices originate from the blade tips and roots, propagating downstream in helical trajectories over a short distance (Sherry et al., 2013). If the inclination angle is minimal, these tip vortices will be interpreted as cylindrical shear layers (Crespo et al., 1999). These layers expand within the wake due to turbulent diffusion, forming a ring-shaped region characterized by high turbulence intensity and substantial velocity gradients. However, the tip vortices will break down because of instability within a short distance downstream. In the far-wake region, the rotor effects are less important. The wake is fully developed; turbulent diffusion of momentum becomes dominant; and an increased level of turbulence is found, fueled by shear production. It is generally assumed that velocity and turbulence intensity should remain self-similar and axi-symmetric in the absence of ambient wind shear (Johansson et al., 2003; Xie and Archer, 2015; Cafiero et al., 2020; van der Laan et al., 2023).

The study by Ishihara and Qian (2018), hereafter referred to as IQ2018, is the first in the literature to provide a three-dimensional analytical equation for the turbulence intensity added by a wind turbine in the downstream wake. Their empirical formula contains a Gaussian function in the radial direction r (from the center of the rotor) multiplied by an amplitude function (an inverse function of x) as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }\mathrm{TI}(x,y,z)& ={\displaystyle \frac{\mathrm{1}}{d+e\frac{x}{D}+f{\left(\mathrm{1}+\frac{x}{D}\right)}^{-\mathrm{2}}}}\\ & \times \left\{k_{\mathrm{1}}\exp \left[-{\displaystyle \frac{{\left(r-\frac{D}{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right]\right.\\ & \left.+k_{\mathrm{2}}\exp \left[-{\displaystyle \frac{{\left(r+\frac{D}{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right]\right\}-\mathit{\delta }\left(z),\right.\end{array}\text{(4)}& {\displaystyle }\mathit{\delta }\left(z\right)=\left\{\begin{array}{ll}\mathrm{0}& z\ge H\\ {\displaystyle }{\mathrm{TI}}_{\mathrm{\infty }}\sin {\left(\mathit{\pi }{\displaystyle \frac{H-z}{H}}\right)}^{\mathrm{2}}& z<H,\end{array}\right.\end{array}$$

where H is the turbine hub height, and $d,e,f,k_{\mathrm{1}},k_{\mathrm{2}},$ and σ are functions of C_T and TI_∞ of the form $a{C}_{\mathrm{T}}^b{\mathrm{TI}}_{\mathrm{\infty }}^c$. We note that Eq. (1) can also be reduced to this same form and Eq. (2) to a close form (with $(\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}}{)}^b$ instead of $C_{\mathrm{T}}^b$).

Li et al. (2022) later hypothesized that the added turbulence intensity ΔTI in the streamwise direction has a similar self-similarity property as the velocity deficit and proposed a three-dimensional analytical formula for added turbulence intensity similar to that by IQ2018. Tian et al. (2022), hereafter referred to as TIAN2022, developed a three-dimensional cosine-shaped model to estimate the wake turbulence intensity; they assumed that the wake has a similar growth rate in the spanwise and vertical directions and that the maximum added turbulence intensity is redistributed along the radial direction with a dual-cosine-shaped function.

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }\mathrm{TI}(x,y,z)=I_{+,m}\left(x\right)\cdot \left\{\begin{array}{ll}{\displaystyle }\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}.{\displaystyle \frac{r-\frac{D}{\mathrm{2}}}{r_w}}\right)& \mathrm{0}<r\le \frac{D}{\mathrm{2}}\\ {\displaystyle }\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}.{\displaystyle \frac{r-\frac{D}{\mathrm{2}}}{r_w-\frac{D}{\mathrm{2}}}}\right)& \frac{D}{\mathrm{2}}<r\le r_w\end{array}\right.\end{array}$$

Here r_w and $I_{+,m}\left(x\right)$ are functions of C_T, TI_∞, and $\left(\frac{x}{D}\right)$ of the form $a{C}_{\mathrm{T}}^b{\mathrm{TI}}_{\mathrm{\infty }}^c(\frac{x}{D}{)}^d$.

Wind turbine wakes have also been studied with numerical wake models. The first numerical wake model was proposed by Lissaman (1979), who pre-dated Jensen (1983) and developed a numerical program to solve the complex problem of overlapping wakes in an array of multiple wind turbines. Lissaman (1979) was the first to recognize the importance of ambient turbulence in overlapping wakes, which, as stated, may have a greater impact than that due to the momentum deficit generated by the individual turbines. Over 10 years later, Ainslie (1988) proposed a numerical solution for wake development by simplifying the Navier–Stokes equations for the turbulent boundary layer and introduced ambient turbulence and its effect on the wake decay. More recently, Yang et al. (2015) proposed a modeling framework of added TKE by a wind turbine (ΔTKE) in vertical planes downwind as a function of inlet velocity at hub height and thrust coefficient, which can be used to estimate the TKE of turbine wakes in complex terrain.

Whether purely analytical or numerical and whether predicting wind speed deficit or added turbulence intensity or both, wake models are useful to understand the behavior of a wind turbine wake, but they cannot provide any information on the effects of the wake on the surrounding environment, such as changes in vertical mixing or surface temperature or heat and momentum fluxes at the surface. Large-eddy simulation (LES) has been a successful numerical approach to study wind turbine wakes (Breton et al., 2017) and their effects on the surrounding environment (Wu et al., 2023) because of their high spatial and temporal resolutions (order of a few meters and tens of seconds, respectively), as well as the accuracy of the actuator disk (Sørensen and Myken, 1992; Madsen, 1996; Mikkelsen, 2003) and actuator line (Sorensen and Shen, 2002) models used to incorporate the effects of the rotating blades. Many LES studies have been conducted to capture wind speed and TKE properties in wind turbine wakes (Eriksson et al., 2015; Vanderwende et al., 2016; Lee and Lundquist, 2017; Deskos et al., 2019; Siedersleben et al., 2020; Feng et al., 2022). Notably, Wu et al. (2023) conducted LES that included the effect of atmospheric stability to show that the wind speed deficit behaves differently from ΔTKE (e.g., the wind speed deficit reaches the ground within 8D, while added TKE remains aloft) and that the two are not co-located in the wake region (e.g., the wind speed deficit peaks at hub height, while added TKE peaks near the rotor tip). However, LESs are computationally demanding and are therefore not used for medium- or long-term wind farm power predictions but rather for temporal horizons of the order of a few hours to a day.

Numerical weather prediction (or mesoscale) models, like the Weather Research and Forecasting (WRF) model (Skamarock et al., 2021), are the preferred tool to predict weather over longer temporal horizons, from several days to several years. However, due to the coarser spatial resolution than that of LES, ranging between 1 and 100 km, numerical weather prediction models cannot resolve the details of the wind turbine wakes, which therefore need to be “parameterized”. A parameterization is a way to include the effects of a process of interest that cannot be resolved directly by the numerical model, typically because the spatial resolution of the numerical model is not fine enough to explicitly treat that process. A parameterization is basically a model within a model that uses the resolved variables at each grid cell to calculate the effects of the process of interest on the resolved variables in that cell (but not the process itself). In the WRF model, several processes are parameterized, including convection, boundary layer turbulence, and radiation, to name a few. The wind farm parameterization (WFP) available by default in WRF is that by Fitch et al. (2012), which treats the wind turbines in a grid cell as sinks of momentum and sources of TKE. As shown in the literature (Pan and Archer, 2018; Archer et al., 2019, 2020; Fischereit et al., 2022), the Fitch parameterization ignores wake effects within a grid cell and treats ΔTKE in an overly simplistic way.

In summary, most studies in the literature have focused on predicting the velocity deficit caused by wind turbines; far fewer have focused on added turbulence. In the far-wake region, which is the most relevant portion of the wake for long-term impacts on the environment, TKE is formed due to the increased shear caused by the wind speed deficit in the upper part of the rotor area. If ΔTKE is not accounted for properly, inaccurate predictions of the turbulent fluxes of heat and momentum near the surface may occur, which ultimately may cause inaccurate predictions of near-surface temperature and moisture. Here, we aim at developing an analytical formulation for ΔTKE by wind turbines that is designed to be ultimately incorporated in a future wind farm parameterization for the WRF and other numerical weather prediction models. With the understanding that any parameterization is, by definition, an approximation, our goal in this paper is to propose a reasonable analytical formulation for added TKE that avoids overprediction. Avoiding overprediction of added TKE is crucial for a (future) parameterization because a numerical mesoscale model, like the WRF, will add some TKE on its own via the production term in the TKE equation due to the weak, additional, and resolved vertical shear caused by the reduced wind speed in the grid cell of the turbines. If the parameterization overestimates TKE in addition to the TKE added from the resolved shear, then an excess of TKE will occur. This issue of potential “double counting” of TKE was first addressed by Ma et al. (2022b, a). Our formulation for ΔTKE is inspired by that proposed by Ishihara and Qian (2018), as their analytical formula for added turbulence intensity by a wind turbine depends on atmospheric stability (via ambient turbulence intensity) and turbine technical specifications (via the thrust coefficient).

2 Methods

TKE is the kinetic energy per unit mass associated with eddies in a turbulent flow, defined as half of the sum of the variances (squares of standard deviations) of the three velocity components $u,v,$ and w (along $x,y,$ and z, respectively) as follows:

$$\begin{array}{}\text{(6)}& \mathrm{TKE}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\sigma }}_u^{\mathrm{2}}+{\mathit{\sigma }}_v^{\mathrm{2}}+{\mathit{\sigma }}_w^{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{v^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}\right),\end{array}$$

where a bar indicates a mean, and a prime refers to a fluctuating component, i.e., the difference between the instantaneous and the mean wind component, e.g., ${\displaystyle }{u}^{\prime }=u-\stackrel{\mathrm{‾}}{u}$.

Turbulence intensities are defined along each direction as follows (Arya, 2001; Burton et al., 2011):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\mathrm{TI}}_x={\displaystyle \frac{{\mathit{\sigma }}_u}{\stackrel{\mathrm{‾}}{U}}}={\displaystyle \frac{\sqrt{\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}}}{\stackrel{\mathrm{‾}}{U}}},{\mathrm{TI}}_y={\displaystyle \frac{{\mathit{\sigma }}_v}{\stackrel{\mathrm{‾}}{U}}}={\displaystyle \frac{\sqrt{\stackrel{\mathrm{‾}}{{v^{\prime }}^{\mathrm{2}}}}}{\stackrel{\mathrm{‾}}{U}}},\\ & {\mathrm{TI}}_z={\displaystyle \frac{{\mathit{\sigma }}_w}{\stackrel{\mathrm{‾}}{U}}}={\displaystyle \frac{\sqrt{\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}}}{\stackrel{\mathrm{‾}}{U}}},\end{array}\end{array}$$

where $\stackrel{\mathrm{‾}}{U}$ is the mean wind speed, generally taken at hub height and often only horizontal. If turbulence was truly isotropic, then the 3 standard deviations should be approximately equal to one another. In the real atmosphere, however, where the x, y, and z axes are aligned with the west–east, south–north, and bottom-up directions, respectively, this is not true. Typically the largest one is σ_u, followed by σ_v (approximately 0.64–0.75σ_u in neutral conditions) and then by σ_w (approximately 0.50–0.52σ_u in neutral conditions) (Arya, 1988; Stull, 2017).

In the wind energy community, since wind turbines are generally yawed to face the mean wind, the x direction is set to coincide with the longitudinal (or streamwise) direction (which is not necessarily aligned west–east), and the y and z directions are denoted as lateral and upward. The International Electrotechnical Commission (IEC) standard for wind turbines (International Electrotechnical Commission, 2019) defines turbulence intensity as “the ratio of the wind speed standard deviation to the mean wind speed”, but it effectively considers only σ_u (referred to as σ₁ or “turbulence standard deviation”, where 1 is the index for the x axis) in its definition of TI. However, the IEC standard recognizes that the 3 standard deviations in Eq. (7), even with the convention of aligning the x axis streamwise, should be different from one another and recommends that any wind velocity field for turbulence models used for standard turbine classes satisfies the following conditions: σ_v≥0.7σ_u and σ_w≥0.5σ_u (International Electrotechnical Commission, 2019, Table C.1, Eq. C.10).

There is not a straightforward relationship between σ_U and the standard deviations of the individual wind components σ_u, σ_v, and σ_w because of the non-linear relationship $U=\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}+w^{\mathrm{2}}}$. It can be shown that, to a first approximation, the standard deviation of the wind speed σ_U is very close to that of the streamwise component σ_u (Larsén, 2022; Archer, 2024), which supports the convention used by the IEC. In light of these considerations, here we define TI as the ratio of the wind speed standard deviation over the mean wind speed at hub height, but we use the following approximations:

$$\begin{array}{}\text{(8)}& \mathrm{TI}={\displaystyle \frac{{\mathit{\sigma }}_U}{\stackrel{\mathrm{‾}}{U}}}\approx {\displaystyle \frac{{\mathit{\sigma }}_u}{\stackrel{\mathrm{‾}}{U}}}\approx {\displaystyle \frac{\sqrt{\frac{\mathrm{2}}{\mathrm{3}}\mathrm{TKE}}}{\stackrel{\mathrm{‾}}{U}}}.\end{array}$$

Traditionally, only the horizontal wind speed is used for U and $\stackrel{\mathrm{‾}}{U}$, possibly because only the horizontal components of wind velocity contribute to the rotation of the wind turbine blades. The approximate relationship between TI and TKE in Eq. (8) (Wilcox, 2006), which is based on the assumption of isotropic turbulence and is therefore likely to overestimate TI, is used to convert between TI and TKE when needed (e.g., in Sect. 3.2). In particular, the relationship used in this study between added TI (ΔTI) and added TKE (ΔTKE) is

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }\mathrm{TI}={\displaystyle \frac{\sqrt{\frac{\mathrm{2}}{\mathrm{3}}\mathrm{\Delta }\mathrm{TKE}}}{\stackrel{\mathrm{‾}}{U}}}={\displaystyle \frac{\sqrt{\frac{\mathrm{2}}{\mathrm{3}}(\mathrm{TKE}-{\mathrm{TKE}}_{\mathrm{\infty }})}}{\stackrel{\mathrm{‾}}{U}}},\end{array}$$

where TKE_∞ is, broadly speaking, the freestream turbulence kinetic energy. The exact definition of TKE_∞ depends on the type and distribution of the available data. If three-dimensional simulation data are available from a run without turbines (i.e., a precursor run) and a run with turbines, then the point-by-point difference in the time-averaged TKE of the two runs is used to calculate ΔTKE, e.g., for the validation LES datasets described in Sect. 2.2. If only a simulation with turbines is available, as is the case for the validation LES datasets described in Sect. 3.2, then the vertical profile of TKE at an upstream distance of $x=x_{\mathrm{0}}-\mathrm{2}D$ is obtained by calculating at each level the average of TKE over $-\mathrm{3}D\le y-y_{\mathrm{0}}\le +\mathrm{3}D$, where x₀ and y₀ are the coordinates of the turbine. The value of TKE_∞ to use at each point downstream is, then, the value of TKE in the upstream vertical profile at the same vertical level.

A notable difference between Eqs. (6) and (8) is that TKE includes the vertical fluctuations of the wind field, which makes TKE more suitable than TI for applications where vertical mixing is important, like wind turbine wake effects. In addition, in mesoscale models like the WRF, TKE is often a prognostic variable that is directly simulated, whereas TI and the standard deviations of the wind components are not (i.e., there is no equation to directly obtain the individual components of the Reynolds stress tensor). In the field, however, it has been common to compute TI from measurements because TI can be calculated from simple and relatively inexpensive two-dimensional cup anemometers (thus only the horizontal components are considered), whereas sophisticated and expensive three-dimensional sonic anemometers are necessary to measure all three components of the wind in order to derive TKE.

In the rest of this study, x is the downstream distance from the wind turbine, z is the vertical distance from the ground, and y is the lateral distance from the wind turbine (positive to the left of the turbine, facing the turbine).

2.1 Proposed formulation

The equation for added TI by Ishihara and Qian (2018) is the starting point of the proposed formulation because it captures well the annular distribution of turbine-induced turbulence and its evolution into a single-peak Gaussian with distance. A few features, however, are not well resolved: the peak of added turbulence at hub height occurs at about 1D from the turbine position, rather than at 4D–6D where shear production is the highest; the peak in the vertical at the rotor top is too strong; and a spurious peak forms below the rotor. Inspired by the formulation of IQ2018, but with the intent of improving upon the issues above, here we propose that ΔTKE, normalized by the square of the upstream undisturbed hub-height wind speed U_∞, can be modeled as the product of three functions: a streamwise function A(x), a radial function G(r), and a vertical function W(z):

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{\Delta }\mathrm{TKE}}{U_{\mathrm{\infty }}^{\mathrm{2}}}}=\mathit{\alpha }\times A\left(x\right)\times G\left(r\right)\times W\left(z\right).\end{array}$$

We note that the formulation by IQ2018 was similar, except it did not include a W(z) function. The scalar α is a tuning parameter that ensures that the amplitude of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ in the wake is of the right magnitude (i.e., matches the LES data, as described later).

The streamwise function A(x) should not decrease exponentially, as is often assumed for turbulence intensity (Quarton and Ainslie, 1990; Crespo and Hernández, 1996; Xie and Archer, 2015; Ishihara and Qian, 2018), because ΔTKE is known to peak at a distance x_max between 4D and 8D from the turbine's streamwise location x₀ (Xie and Archer, 2015; Wu et al., 2023), not at x₀ (Fig. 1a). Here we propose a Weibull-like distribution for A(x) as follows:

$$\begin{array}{}\text{(11)}& A\left(x\right)={\left({\displaystyle \frac{x-x_{\mathrm{0}}}{{\mathit{\lambda }}_A}}\right)}^{k_A-\mathrm{1}}\exp \left[-{\left({\displaystyle \frac{x-x_{\mathrm{0}}}{{\mathit{\lambda }}_A}}\right)}^{k_A}\right],\end{array}$$

where λ_A and k_A are the scale and shape parameters of the Weibull distribution. The Weibull function is chosen because it is non-symmetric and because it is one tailed, as is the observed distribution of added TKE along x. The Weibull distribution was also recently proposed for the x dependency of added TI by Delvaux et al. (2024, their Eq. 3). We set k_A=2 to reduce the overall number of parameters to fit and to obtain a function with the desired properties, i.e., equal to zero at x₀ (thus k_A>1) and rapidly increasing past x₀ but not too rapidly (which would be the case for k_A<2). An example of the evolution of A(x) is shown in Fig. 1a in blue.

Figure 1 Comparison of the proposed fitting functions against the WRFLES-N results by Wu et al. (2023) (labeled “LES”): (a) streamwise function A(x) at z=H and $y=y_{\mathrm{0}}+D/\mathrm{2}$ versus the fit by Ishihara and Qian (2018) (“IQ2018”) and Crespo and Hernández (1996) (“CH1996”) and (b) radial and vertical functions G(r) and W(z) at $x=x_{\mathrm{0}}+\mathrm{5}D$ and y=y₀. Note that G(r) and W(z) are re-scaled for display purposes. The thin gray lines in (b) indicate the rotor top and bottom.

The radial function G(r) is assumed to be a Gaussian that peaks at the tip annulus of the rotor, inspired by the ϕ(r) function of IQ2018, as follows:

$$\begin{array}{}\text{(12)}& G\left(r\right)=\exp \left[-{\displaystyle \frac{{\left(r-D/\mathrm{2}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_r^{\mathrm{2}}}}\right],\end{array}$$

where r is

$$\begin{array}{}\text{(13)}& r=\sqrt{{\left(y-y_{\mathrm{0}}\right)}^{\mathrm{2}}+{\left(z-H\right)}^{\mathrm{2}}},\end{array}$$

y₀ is the spanwise location of the turbine, and σ_r is a linear function of x:

$$\begin{array}{}\text{(14)}& {\mathit{\sigma }}_r\left(x\right)=k_r\left(x-x_{\mathrm{0}}\right)+{\mathit{\epsilon }}_{r}D,\end{array}$$

where k_r is the radial expansion rate (i.e., $\frac{\partial {\mathit{\sigma }}_r}{\partial x}$), and ε_r is a multiplying factor to the rotor diameter that sets the initial width of the Gaussian distribution of the added TKE along the annulus of the rotor disk.

Lastly, the vertical function W(z) is also assumed to be Weibull-like:

$$\begin{array}{}\text{(15)}& W\left(z\right)={\left({\displaystyle \frac{z}{{\mathit{\lambda }}_W}}\right)}^{k_W-\mathrm{1}}\exp \left[-{\left({\displaystyle \frac{z}{{\mathit{\lambda }}_W}}\right)}^{k_W}\right],\end{array}$$

where the shape parameter k_W is set equal to 4 after a trial-and-error process to ensure a steeper decrease in ΔTKE above the top tip than below it, as shown in the LES results. Examples of vertical profiles of the functions G and W are shown in Fig. 1b with orange and purple lines.

Due to the properties of the Weibull distribution, there is an analytical relationship between the point at which the function reaches the maximum, which by definition is the mode x_m, and the value of λ:

$$\begin{array}{}\text{(16)}& \mathit{\lambda }={\displaystyle \frac{x_m}{{\left(\frac{k-\mathrm{1}}{k}\right)}^{\frac{\mathrm{1}}{k}}}}.\end{array}$$

Thus, the values of λ_A and λ_W are directly related to the position of the maximum ΔTKE along x and z, respectively. As such, it is reasonable to expect that λ_A and λ_W depend on D and H.

In summary, the equation for $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ (Eq. 10) contains five unknown parameters: $\mathit{\alpha },{\mathit{\lambda }}_A,{\mathit{\lambda }}_W,k_r$, and ε_r. Before we explain how we obtain their values (by fitting) and why we refer to them hereafter as the five “direct-fit” parameters (in Sect. 2.3), we need to describe the datasets that we used for the fitting (in Sect. 2.2).

2.2 LES datasets

We use two independent LES datasets to calibrate the fitting parameters of our analytical model (Eqs. 10–15): the published results of Wu et al. (2023) and Vahidi and Porté-Agel (2025), described next.

The modeling system used by Wu et al. (2023) included both an LES model (WRF-LES) (Moeng et al., 2007) and a generalized actuator disk parameterization (WRF-LES-GAD) (Mirocha et al., 2014). A single wind turbine and a row of wind turbines were simulated under neutral, stable, and unstable conditions with different undisturbed hub-height wind speeds, turbulence intensities, and thrust coefficients (Table 1). The outer domain with a horizontal grid spacing of 15 m was one-way nested with an inner domain with a finer horizontal grid spacing of 5 m (see their Fig. 1). The vertical resolution started at 2.5 m near the ground and was then stretched by 10 % per level until 35 m, kept constant at 5 m from 35 to 200 m, and finally stretched by 5 % until the model top. The lateral boundary conditions were periodic for the outer domain and time-varying for the inner domain. The time step was $\mathrm{1}/\mathrm{7}$ s in the outer domain and $\mathrm{1}/\mathrm{21}$ s in the inner domain. The wind turbine was the PSU 1.5 MW turbine with H=80 m and D=77 m, placed in the domain 8.2D from the inlet of the inner domain and centered laterally. The temperature profile was uniform at 297.3 K up to the initial boundary layer height (at 200, 500, and 800 m for stable, neutral, and unstable conditions, respectively), with an inversion of 0.01 K m⁻¹. The desired atmospheric stability was imposed by assigning a surface heat flux of −0.07, 0, and +0.07 K m s⁻¹ for stable, neutral, and unstable conditions, respectively, with a surface roughness length z₀= 0.01 m for all simulations.

Table 1 Simulation details and values of the five direct-fit parameters for the 17 LES cases used for the analytical model calibration. The label “WRFLES” refers to the WRF-LES dataset by Wu et al. (2023), where “S”, “N”, and “U” refer to stable, neutral, and unstable conditions; “VPA-TI064” and “VPA-TI107” refer to the neutral-stability simulations by Vahidi and Porté-Agel (2025), with TI set to 0.064 and 0.107.

The LES results of Vahidi and Porté-Agel (2025) were obtained with the WiRE-LES solver (Abkar and Porté-Agel, 2015), with a standard actuator disk model (Meyers and Meneveau, 2010) for the axial force induced by the wind turbine. The computational domain was 3840 m × 1920 m × 640 m with 256 × 256 × 128 grid points in the x, y, and z directions, respectively, with z₀= 0.001 m and no Coriolis force. The wind turbine, with $D=H=$ 80 m, was located at a distance of 15D from the inlet and centered laterally. The boundary conditions in the horizontal direction were periodic for the precursor runs (i.e., with no wind turbine). When the turbine was inserted in the domain, an inflow boundary condition was applied to override the imposed periodic boundary condition in the streamwise direction using the prior periodic results. To smoothly adjust the flow to an undisturbed inflow condition, a buffer zone was introduced upstream of the inflow section. A constant streamwise pressure gradient was used to maintain a velocity of around 8 m s⁻¹ at the center of the actuator disk.

Further details about the LES suites that were used for calibration can be found in the original studies (Wu et al., 2023; Vahidi and Porté-Agel, 2025) and in Table 1.

2.3 Least-squares error fitting procedures

We use the Python least-squares error non-linear fitting function in two steps. First, we run the fitting separately for each of the 15 LES cases over the wake region of the various computational domains using Eqs. (10)–(15) to obtain 15 sets of the five direct-fit parameters, $\mathit{\alpha },{\mathit{\lambda }}_A,{\mathit{\lambda }}_W,k_r$, and ε_r (Table 1). We refer to them as direct-fit parameters to distinguish them from those obtained in the next step. We also attempted to find fitted values for k_A and k_W, but with seven parameters we could never reach convergence of the least-squares error fitting procedure. Then, we run the Python fitting again using simple functions of C_T and TI_∞ with three fitting coefficients each, described below, to obtain five functional relationships between the five direct-fit parameters and the two independent variables C_T and TI_∞.

Because it is unpractical to use a different set of fitting parameters for each case and because we want a smooth (i.e., not stepwise) transition between different values of C_T and TI_∞, we want to identify the functional relationships of the five direct-fit parameters on a few relevant variables that are turbine and stability dependent. There are a few empirical formulations that have been proposed in the literature for ΔTI and that may be applicable for our purposes. These empirical formulations depend on both C_T and TI_∞ with the following general form (Quarton and Ainslie, 1990; Crespo and Hernández, 1996; Ishihara and Qian, 2018):

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }\mathrm{TI}\propto a{C}_{\mathrm{T}}^b{\mathrm{TI}}_{\mathrm{\infty }}^c.\end{array}$$

Inspired by these well-established empirical formulations, here we propose that, of the five fitting parameters in Eqs. (11)–(15) ($\mathit{\alpha },{\mathit{\lambda }}_A,{\mathit{\lambda }}_W,k_r$, and ε_r), three have the form shown in Eq. (17), namely α, k_r, and ε_r, while λ_A and λ_W include an additional dependency on the relevant lengths D and H as follows:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathit{\lambda }}_A=aD{C}_{\mathrm{T}}^b{\mathrm{TI}}_{\mathrm{\infty }}^c,\text{(19)}& {\displaystyle }{\mathit{\lambda }}_W=H+aD{C}_{\mathrm{T}}^b{\mathrm{TI}}_{\mathrm{\infty }}^c.\end{array}$$

Intuitively, for a wind turbine with a small diameter, the peak of added TKE occurs at a downstream distance from the tower that is shorter than that for a wind turbine with a large diameter; thus λ_A should depend on D, as in Eq. (18). Similarly, the TKE peak in the vertical occurs near the rotor top, thus at higher elevations for taller turbines than for shorter ones, which suggests a dependency of λ_W on H and D, as in Eq. (19).

In summary, the five functional relationships for the five direct-fit parameters are shown in the first row of Table 2. We note that there are now 3 × 5 = 15 functional coefficients a_i, b_i, and c_i, the values of which we obtain via another least-squares error fitting. The results are shown in Fig. 2 and Table 2. Although the functional coefficients are empirical in nature, their values provide useful physical indications of the properties of added TKE, as discussed next.

Table 2 Functional relationships for the five direct-fit parameters. The values in parentheses are the original coefficients before manual overwriting to zero.

Figure 2 Performance of the functional relationships in predicting the five direct-fit parameters: (a) α, (b) k_r, (c) λ_A, (d) ε_r, and (e) λ_W.

Starting with α (Fig. 2a), its functional relationship is extremely consistent among all the runs and independent of TI_∞; thus the value of c, originally equal to 0.037, is set to zero in Table 2 to simplify the functional relationship and reduce the number of coefficients required overall. The implication is that the maximum added TKE in the wake of a wind turbine is essentially independent of atmospheric properties (such as turbulence intensity or stability) but depends only on the turbine operation through its thrust coefficient. This finding is consistent with several published studies. Crespo and Hernández (1996) proposed that the maximum added TI depends only on C_T in the near wake and has a very weak dependency on TI_∞ in the far wake (i.e., the power exponent is −0.0325; see Eq. (2)); the models for maximum added TI by both Larsen et al. (1996) and Frandsen (2007) depend only on C_T; and the formulations for the maximum added TI by IQ2018, shown in Eq. (3), and by TIAN2022, shown in Eq. (5), have a weaker dependency on TI_∞ than on C_T. However, other studies have proposed that TI_∞ has a non-minor role (Ainslie, 1988; Xie and Archer, 2015).

By contrast, the fit for k_r (Fig. 2b) is basically independent of C_T, despite a weak and conflicting dependency in the direct-fitting values; thus b is overwritten as zero from the original value of −0.061 in Table 2. This indicates an interesting finding: the radial expansion rate of the wake TKE is independent of the turbine operation but is only a function of the amount of background turbulence. This finding is consistent with the literature, as k_r has a similar meaning to the well-known expansion rate k_w (also known as just k) in the Jensen model (Jensen, 1983), which is typically set to 0.075 onshore, where background turbulence is generally high, and 0.04 offshore, where turbulence is low (Archer et al., 2018). We note that k_r is the expansion of the wake TKE, while k_w is the expansion of the wake wind speed deficit. As a result of the analysis of the functional relationships for α and k_r, the number of fitting coefficients is reduced from 15 to 13.

Next, ε_r is proportional to both C_T and TI_∞, but the dependency on C_T is stronger (b is more than twice as large as c, Table 2). Since ε_r controls the spread of the added TKE distribution along the annulus of the rotor disk, it is not surprising that its value for the unstable case (WRFLES-U) is underpredicted, while for the stable case (WRFLES-S) it is overpredicted in Fig. 2d.

Focusing on λ_A next, since the values of both b and c are negative (Table 2), λ_A is inversely proportional to both C_T and TI_∞. However, since both C_T and TI_∞ are always lower than 1 and b and c are both negative and lower than 1 but c is larger in magnitude, the dependency on TI_∞ ends up being dominant. This finding, too, is physically sound: everything else being the same, the downstream peak of added TKE is expected to be closer to the wind turbine when the background turbulence is high and further downstream when the background turbulence is low because high turbulence causes a shorter wake than low turbulence, increases mixing, and smooths down the peak.

By contrast, λ_W is inversely proportional to C_T and directly proportional to TI_∞ (Fig. 2e), consistent with negative b and positive c in Table 2. However, due to the additional dependency on H and D, the comparison between the lines in Fig. 2e needs to be conducted for the same D and H, thus only among lines of the same LES group. For example, looking at the VPA runs only, a near doubling of TI_∞ from VPA-TI064 to VPA-TI108 causes an approximately 10 m rise in the vertical placement of the TKE peak, while a near doubling of C_T from 0.4 to 0.7 causes a less than 5 m drop in the peak position.

The two-step approach described so far may appear cumbersome, with an initial least-squares error fit to the LES data for the five direct-fit parameters and then a second least-squares error fit with the functional relationships to finally obtain the values of the 15 desired functional coefficients. A more straightforward approach would have been to put the functional relationships shown in Table 2 directly into Eqs. (11)–(15) and then to apply the least-squares error fitting to the LES data for the 15 unknown functional coefficients. However, when we tried it, we were unable to reach convergence, possibly because the number of fitting parameters was too high (15 versus 5). We note that IQ2018 used a total of nine fitting parameters.

To assess the performance of the two-step approach, we calculate the root mean square errors (RMSEs) over the entire wake region of all calibration cases (Table 3, calibration cases only). The entire wake region covers an area with y from −1D to +1D and z from 0–2H. Note that we have different lengths along the x direction as some LES cases and the experimental test only cover a small region downstream. The proposed analytical formulation outperforms that by IQ2018 in all cases, as the RMSEs of both the direct and the final fit are half as large as or lower than theirs on average and up to 6 times smaller. A general trend that emerges is that the RMSE increases as both C_T and TI_∞ increase. For example, the RMSE of the VPA runs approximately doubles when TI_∞ increases from 0.064 to 0.108 and triples when C_T varies from 0.4 to 0.9. Not surprisingly, the RMSE of the final fit is always higher than that of the direct fit, as, by definition, the direct-fit parameters are those that minimize the error. We note that the RMSEs are relatively large and close to or slightly higher than the mean value of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$. For example, the RMSE of the final fit for the WRFLES runs is about $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{3}}$, and the mean value of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ for the same runs is about $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$.

Table 3 RMSE values of the 15 LES calibration cases and the four validation cases against the two benchmark models proposed by Ishihara and Qian (2018) and by Tian et al. (2022), labeled “IQ2018” and “TIAN2022”, respectively, and the direct fit (with direct-fit parameters from Table 1) and the final fit (with the functional relationships from Table 2). The RMSEs are calculated over a longer domain for the validation cases (0D–20D) than the calibration cases (up to 10D; see text for details).

Note that n/a represents not applicable.

3 Results

3.1 Comparison with LES data

First, we look at the individual functions A,G, and W along one relevant dimension for the calibration case WRFLES-N. For A, the relevant dimension is x, and the proposed formulation exhibits the correct features (Fig. 1a): the function is zero at x₀; it rapidly increases and then peaks at about 6D, slightly further downstream than indicated by the LES profile, which peaks at about 5D; and then it slowly decreases to nearly zero at 20D. Both the IQ2018 and the CH1996 curves peak near or at the rotor (at 0D) and retain too much TKE in the far wake after 10D.

The vertical profiles of the two functions G and W (Fig. 1b) are also, for the most part, correctly reproduced by the proposed formulation, including the peak of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ near the rotor tip and the rapid decrease above it. The second peak of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ located below hub height is qualitatively reproduced but is too weak; $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ below the rotor decreases without exhibiting the weak negative peak near the surface. The fit by IQ2018 greatly overestimates TKE near the rotor top (by a factor of approximately 2) and exhibits a spurious peak near the surface.

Next, we compare horizontal and vertical cross-sections of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ from the LES results, direct fit, and final fit for each of the three stability cases. Starting with a neutral case (VPA-TI064 with C_T=0.7) in Fig. 3, the horizontal cross-sections at hub height show that the main features and distribution of the direct and final fit for $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ resemble those of the LES results, with a symmetric distribution around the wake axis and two maxima near 6D (slightly further downstream than the LES maxima at ≃5D, as mentioned earlier for Fig. 1a). The two maxima collapse into one in the upper part of the rotor (Fig. 3, middle).

Figure 3 Cross-sections of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ under neutral conditions: (a) LES (VPA-TI064 with C_T=0.7); (b) direct fit; and (c) final fit in the x–y plane at hub height (left), x–y plane at the rotor top (middle), and x–z plane at y=y₀ (right).

In the x–z plane, the $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ maximum is properly located at about 120 m above ground ($z/H\sim \mathrm{1.5}$); the elongated feature of higher $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ extending towards the lower rotor area is more or less captured by the proposed formulation (Fig. 3, right). While the location and magnitude of the maximum, as well as the overall distribution above the rotor, are captured well by both fittings, the wake extent in the x direction is underestimated and so is the magnitude of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ below hub height. We note that the magnitude of the maximum $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ is better reproduced by the final fit than by the direct fit.

The comparison between $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ from the LES and the direct and final fits under unstable conditions is also encouraging (Fig. A1 for WRFLES-U). While the direct fit closely mimics the maximum value of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ under unstable conditions across all three cross-sections, the final fit underestimates this term (Fig. A1); the opposite occurs under neutral conditions in Fig. 3. However, the final fit reproduces the $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ propagation along x better than the direct fit.

For stable conditions (Fig. A2), the performance of the proposed fits is more complicated. At hub height, the direct fit underestimates the magnitude of the maximum $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$, and the final fit improves it (Fig. A2, left). At the rotor top and in the z–x plane (Fig. A2, middle and right), the opposite happens: the direct fit matches the maximum well, while the final fit underestimates it. It appears that the final fit shifts the maximum of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ further down in elevation, which causes an increase in TKE around hub height and a reduction near the rotor top.

It is important to note that we did not yet include a treatment of the hub, which causes high ΔTKE between 1D and 2D in the LES results. The well-known feature that $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ extends further downstream under stable conditions compared to unstable conditions is correctly captured by our proposed formulation, although the magnitude of the peak is underestimated by the fits in both cases.

Lastly, we compare the vertical and horizontal profiles of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ from the proposed final fit against the formulation of IQ2018 and the LES results of Wu et al. (2023) under neutral, stable, and unstable conditions at different downstream distances (Fig. 4).

Figure 4 Profiles of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ from the proposed formulation (“Fit”) and from the WRFLES runs (Table 1) along (a) z and (b) y at different downstream distances.

In general, IQ2018 overestimates the LES results, especially the magnitude of the $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ peak in the near wake in the vertical (by a factor of 2, Fig. 4a) and in the horizontal (by a factor of 3, Fig. 4b). Large overestimates by the IQ2018 model have also been reported recently by Bastankhah et al. (2024) in the near- and far-wake regions. By contrast, the proposed fit tends to underestimate the maxima in the near wake by up to 30 %. In the far wake, the proposed formulation predictions are closely aligned with the LES results; the overestimation by IQ2018 is reduced, but a secondary spurious maximum appears near the surface. Since the LES results include the effect of the hub, while both the final fit and the IQ2018 formulation do not, they both miss the peak in $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ at x=2D caused by the hub (Figs. 3–A1, left and right sub-figures).

3.2 Validation

In this section, we compare the performance of the proposed formulation for $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ (i.e., with the final fitting coefficients from Table 2) with data from four other independent studies: the LES studies by Xie and Archer (2017) and Archer et al. (2020), referred to hereafter as XA2017 and ARC2020, respectively; a modified version of the LES by Archer et al. (2013) and Ghaisas et al. (2017), referred to hereafter as SOWFA; and the wind tunnel measurements by Aju et al. (2020), referred to as AJU2020. The four studies are introduced and discussed in the following, and their relevant parameters are listed in Table 4. The RMSEs for the validation cases are calculated over their respective entire domains starting at 0D to the last point available, which is different in each case (e.g., 10D for SOWFA but just 4D for AJU2020).

Table 4 Details of the datasets used for the analytical model validation. The label “XA2017” refers to the unstable LES results by Xie and Archer (2017), “ARC2020” the neutral LES results of Archer et al. (2020), “SOWFA” a modified version of the LES results by Archer et al. (2013) and Ghaisas et al. (2017), and “AJU2020” the neutral wind tunnel experiments of Aju et al. (2020).

XA2017 used the Wind Turbine and Turbulence Simulator (WiTTS) (Xie and Archer, 2015), an in-house flow solver developed at the University of Delaware, for the flow simulation. WiTTS solves the unsteady, filtered three-dimensional Navier–Stokes equations in the incompressible form using the fractional method (Kim and Moin, 1985) and can handle non-neutral stabilities (Xie and Archer, 2017). The wind turbine, modeled as an actuator line (and the nacelle), was the REpower 5 MW turbine, with a hub height of 87.6 m and a rotor diameter of 126 m. The simulation was divided into two stages: a precursor stage (without turbines) and a formal stage (with the inclusion of five wind turbines). The five turbines were arranged in a staggered layout with along-wind spacing of approximately 1000 m (which is roughly 8D) and across-wind spacing of roughly 4D. The resolution was 6.25 m in all directions, and the surface roughness length was 0.016 m. Only the wake field of the front-row turbine (labeled “WT1” in their Fig. 7) from 2D to 8D is used here to avoid contamination from the wakes of nearby turbines.

Looking at the results from the unstable run of XA2017, $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ from the proposed fit matches the LES satisfactorily in the vertical and in the horizontal (Fig. 5a and b), with the maximum $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ having the right magnitude and being correctly located at the upper part of the rotor tip at 4D and 6D. However, at 8D, the profile along the z direction shows an overestimation in predicting the maxima near the rotor tip, while correctly reproducing the profile below hub height. The RMSE is larger than that of any of the studies used for calibration, $\mathrm{3.77}\times {\mathrm{10}}^{-\mathrm{3}}$ (Table 3), possibly because the LES included the effect of the nacelle, while both analytical formulations do not. The RMSE of the proposed fit is, however, significantly lower than that of IQ2018 ($\mathrm{12.36}\times {\mathrm{10}}^{-\mathrm{3}}$), which overestimates added TKE at all distances, but especially in the near wake, and introduces, again, an unrealistic peak near the ground.

Figure 5 Profiles of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$, taken at different downstream distances, from the proposed formulation, labeled “Fit”, and from the LES dataset of XA2017 along (a) the z direction at the centerline and (b) the y direction at hub height, under unstable conditions. (c) The wind tunnel particle image velocimetry (PIV) measurements of AJU2020, along the z direction at the centerline, under neutral stability. The same profiles from IQ2018 and TIAN2000 are added for comparison.

The second LES study used here for validation is ARC2020 (Archer et al., 2020), in which the LES flow solver Software for Wind Farm Applications (SOWFA) was used with an actuator line model for the wind turbine blades without any treatment for the hub (Churchfield et al., 2012; Archer et al., 2013). The computational domain size was 3000 m × 3000 m × 1020 m with a single wind turbine in the middle, and the atmospheric stability was neutral. An idealized NREL 5 MW wind turbine was used with D=126 m and H=90 m. The resolution was set to 200 × 200 × 68 grid points in x, y, and z, respectively, corresponding to grid cells of approximately 15 m in both horizontal dimensions.

The vertical profiles of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ from the final fitting match those from the LES just after x=2D very well (Fig. 6a). We note that the slight reduction in TKE near the surface shown in the LES results, which has been observed and simulated in the literature (Archer et al., 2019; Wu and Archer, 2021), is not reproduced with the proposed fit because it is not accounted for in its equations yet. A way to account for it in the future could be via a correction similar to the δ function of IQ2018, shown here in Eq. (4). Conversely, IQ2018 shows some spurious $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ below the rotor, which is not correct. The y profiles (Fig. 6b) are again characterized by a general underestimation and overestimation of the $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ maxima for the proposed fit and for IQ2018, respectively. The proposed fit shows a significantly improved $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ profile compared to IQ2018, providing a more satisfactory prediction that closely matches the LES results starting at x=4D. This validation dataset is well reproduced by the proposed formulation, with an RMSE of $\mathrm{3.63}\times {\mathrm{10}}^{-\mathrm{3}}$, about half of the average value of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ ($\mathrm{6.46}\times {\mathrm{10}}^{-\mathrm{3}}$, Table 3), possibly because this LES study did not include a treatment for the hub.

Figure 6 Profiles of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$, taken at different downstream distances, from the proposed formulation, labeled “Fit”, and from the LES dataset of ARC2020 under neutral stability along (a) the z direction at the centerline and (b) the y direction at hub height. The same profiles from IQ2018 and TIAN2022 are added for comparison.

The third study used for validation, named SOWFA, is a modification of the studies by Archer et al. (2013) and Ghaisas et al. (2017), which used the SOWFA solver over a complex mesh of 4000 m × 4000 m × 1000 m with fine refinement (about 3.5 m) in six blocks around up to 48 wind turbines (Siemens 2.3 MW with D=93 m and H=63.4 m) and coarser (7 m) in the rest of the domain. Various cases were simulated, varying the number of turbines, their layout, and the atmospheric stability. Here we use an additional stable case, with the same temperature decrease at the bottom boundary of −0.25 K h⁻¹ as in Ghaisas et al. (2017) and the same layout as “Stg-2SpaX” in Archer et al. (2013) but with westerly flow. To maximize the extent of the wake at fine resolution, we extracted the flow details from the wake of turbine no. 36 (see Fig. 1f in Archer et al., 2013), approximately 10D in length.

Figure 7 The same as Fig. 6 but for the SOWFA LES dataset under stable conditions.

The proposed fit provides a satisfactory match to the LES profiles, particularly in the z profiles, in the entire domain, while a noticeable gap between the LES data and IQ2018 still exists in the near wake. We note that the stable LES results show an asymmetry in the y profiles (Fig. 7b) due to the Coriolis force and the resulting Ekman-spiral effect, which is not captured by the proposed formulation, and therefore an overestimation by the final fitting can be seen in the far wake. Because this LES dataset did not simulate the effect of the hub, the RMSE was relatively low: $\mathrm{2.86}\times {\mathrm{10}}^{-\mathrm{3}}$, for an average value of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ of about $\mathrm{4.41}\times {\mathrm{10}}^{-\mathrm{3}}$ (Table 3).

The fourth study used for validation is AJU2020 (Aju et al., 2020), which describes the Boundary Layer and Subsonic Tunnel (BLAST) experiments that were conducted at the University of Texas at Dallas. The wind tunnel is 30 m long, 2.1 m high, and 2.8 m wide, with cubic blocks of 2.5 cm height located 0.2 m between each other on the bottom surface of the test section to achieve a well-developed turbulent boundary layer. Horizontal-axis wind turbines, with D=200 mm and H=200 mm, were used in the experiments, based on models from Sandia National Laboratories, and locally manufactured at the university. A particle image velocimetry (PIV) system was used to measure the turbulence along the center axis of each turbine. The tests were conducted in neutral stability, and only vertical profiles at selected downstream distances were available. Values near the surface are not reliable due to ground effects.

Comparison of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ behind the wind turbine against the experimental data by AJU2020 under neutral conditions indicates that the proposed fit is qualitatively correct but exhibits a large underestimation of the upper-rotor peak by over 100 % (Fig. 5c). Not surprisingly, the AJU2020 dataset is associated with the largest RMSE among the four cases used for validation, exceeding $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}$ for an average value of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ of $\mathrm{5.87}\times {\mathrm{10}}^{-\mathrm{3}}$ (Table 3). This is possibly explained by the fact that this experimental dataset just covers a small region behind the turbine, between $x=\mathrm{0}-\mathrm{4}D$, where the effect of the hub is more significant. By contrast, this is the only dataset that compares well against the IQ2018 predictions, despite an overestimation of the maximum at 3D and an overall RMSE that is still larger than that of the proposed fit ($\mathrm{6.67}\times {\mathrm{10}}^{-\mathrm{3}}$, Table 3).

4 Conclusions and recommendations

This study is a first step in addressing the lack of a proper treatment of the turbulence added by wind turbines in current numerical weather prediction models, like the WRF model. An analytical formulation for $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ is presented, comprising five fitting parameters, each with a functional relationship with the thrust coefficient of the turbine, the undisturbed upstream turbulence intensity, the diameter, and the hub height of the wind turbine. The fitting parameters are obtained after a two-step fitting process based on the LES dataset from a previous study by Wu et al. (2023), which used the WRF-LES code for three atmospheric stability cases (stable, neutral, and unstable), and from 15 LES cases in neutral stability with various combinations of TI_∞ and C_T by Vahidi and Porté-Agel (2025).

The proposed formulation compares well with the LES datasets that were used for the parameter calibration, which is to be expected, with RMSEs of the same order of magnitude as the mean $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$, but it is less accurate when compared against four additional and independent datasets used for validation: an LES study of a 5 MW wind turbine under unstable conditions using the WiTTS code (Xie and Archer, 2017), another LES study of the same 5 MW wind turbine under neutral stability using the SOWFA code (Archer et al., 2020), another SOWFA run under stable conditions (Archer et al., 2013; Ghaisas et al., 2017), and a wind tunnel experiment with a model wind turbine under neutral conditions (Aju et al., 2020).

We conclude that the proposed formulation is promising in predicting the distribution of $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ under all stabilities in the far wake, which is the more relevant region for mesoscale studies of the impacts of wind farms on the environment. In the near wake, the blade geometry, rotor tip, and hub effects have a dominant effect on $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$, and the proposed formulation performs worse there than in the far wake. However, the proposed formulation outperforms that by Ishihara and Qian (2018) in all cases.

The ultimate goal of this research is to eventually insert the $\mathrm{\Delta }\mathrm{TKE}/U_{\mathrm{\infty }}^{\mathrm{2}}$ formulation, after further improvements to better capture the near-wake behavior, in numerical weather prediction models to better quantify the possible impacts of wind turbine wakes on the environment. However, in order for the proposed formulation to be effectively used for this purpose, the total ΔTKE in each grid cell of the mesoscale model needs to be calculated, but the volume integral of Eq. (10) cannot be obtained analytically. As such, numerical integration is required, which may add a small computational cost to the simulation. In addition, in the presence of multiple wind turbines with multiple overlapping wakes, the issue of superposition of wake-added TKE needs to be resolved.

We note that we did not include a treatment of either the nacelle or the tower in our formulation, which, combined, were found to have an effect on the wake dynamics that was larger than previously thought, in the case of a model wind turbine in a wind tunnel (Santoni et al., 2017). Since only XA2017 included the nacelle in their simulation, while none of the validation or calibration studies included the tower, we neglected both for now.

Another limitation of our formulation is that its calibration relies on LES results, which introduce several uncertainties, from the sub-grid turbulence closure model to the sampling method for the actuator line model (Martínez-Tossas et al., 2017). Using only real measurements would not remove all uncertainties either, as measurements have their own intrinsic uncertainties, and each experiment tends to be specific to the chosen setup and therefore difficult to generalize. Even if we used direct numerical simulation (DNS), which we cannot do yet due to the high Reynolds number of the wind flow, resolving the blades correctly would still require an actuator line model or similar parameterization, which would add some uncertainty.

Appendix A: Additional figures

Appendix A contains Figs. A1 and A2 that are briefly discussed in the main text.

Figure A1 As in Fig. 3 but for unstable conditions (WRFLES-U).

Figure A2 As in Fig. 3 but for stable conditions (WRFLES-S).