Analytical modelling of wind-turbine wake turbulence in neutral atmospheric boundary layers

TL;DR

This paper introduces a new physically-based model for wake-added turbulence in wind farms, improving accuracy over empirical models by incorporating turbulence budgets and matching LES and wind tunnel data.

Contribution

A novel wake-added turbulence model based on turbulence budgets, offering a more physically consistent and accurate approach for wind farm flow simulations.

Findings

Model agrees well with LES and wind tunnel data.

Provides a practical and physically consistent turbulence prediction.

Enhances wind farm flow modelling accuracy.

Abstract

So-called engineering or analytical wind farm flow solvers typically build upon two submodels: one for the velocity deficit and one for the wake-added turbulence intensity. While velocity deficit modelling has received considerable attention, wake-added turbulence models are less prevalent in comparison. Yet, accurate estimates of local turbulence intensity are essential for predicting flow interactions and energy yield, as turbine wakes are both sensitive to, and sources of turbulence. Existing wake-added turbulence models are typically empirical or assume axial symmetry despite the inherently three-dimensional nature of turbulent wake fields. In this work, we present a new model for wake-added turbulence intensity. Our approach is based on the analysis of the TKE and the streamwise Reynolds stress budget, incorporating classical RANS modelling assumptions and far-wake approximations. The resulting model maintains a simple and practical form, demonstrating strong agreement with LES and wind tunnel measurements. Our model provides a more physically consistent and predictive tool for wind farm flow modelling and performance estimation.

Figures

Tables

Table 1. Table 1: Description of the analytical models used in section 6 to predict the wake velocity and the wake-added turbulence.

Model Component	Equation / Definition	Constants & Sub-parameters
1. Wake Velocity Model
Streamwise velocity $\overline{u}$	$\overline{u}{}_{\infty }{}^{}\left(z\right)(1-C\left(x\right)\exp (-\frac{1}{2}\frac{r^n}{{\sigma }^2}))$	$r=\sqrt{{\left(y-y_h\right)}^2+{\left(z-z_h\right)}^2}$
Inflow Profile $\overline{u}{}_{\infty }$	$\frac{\kappa }{u_{\ast }}\log \left(\frac{z}{z_0}\right)$
Super-Gaussian Order $n$	$a_f{e}^{b_{f}x}+c_f$	$a_f=-8.2635C_T^3+8.5939C_T^2$ $-8.9691C_T+10.7286$ $b_f=1.68\exp \left(-25.98I_x\right)-1.06$ $c_f=2$
Wake Width $\sigma$	$\left(a_s{I}_x+b_s\right)x+c_s\sqrt{\frac{1}{2}\frac{1+\sqrt{1-C_T}}{\sqrt{1-C_T}}}$	$a_s=0.28$ $b_s=0.01$ $c_s=0.1C_T+0.1$
Maximum velocity deficit $C\left(x\right)$	$2^{2/n-1}-\sqrt{2^{4/n-2}-\frac{n{C}_T}{16\Gamma \left(2/n\right){\sigma }^{4/n}}}$
2. Wake-Added Turbulence Model
Wake added $\overline{u^{\prime }{u}^{\prime }}{}_w$	${\zeta }_{\overline{u^{\prime }{u}^{\prime }},w}\left({\nu }_{t,\infty }\left[{\left(\frac{\partial \overline{u}}{\partial y}\right)}^2+{\left(\frac{\partial \overline{u}}{\partial z}\right)}^2-{\left(\frac{\partial {\overline{u}}_{\infty }}{\partial z}\right)}^2\right]-{\gamma }_{\overline{u^{\prime }{u}^{\prime }}}\Delta \overline{u}\right)$
Parameter ${\gamma }_{\overline{u^{\prime }{u}^{\prime }}}$	$5\times {10}^{-3}{C}_T^2\left({\overline{u}}_h^2/D\right)$
Parameter ${\zeta }_{\overline{u^{\prime }{u}^{\prime }},w}$		$a_{\zeta }=0.65$ $b_{\zeta }=-0.8$
Inflow Eddy Viscosity ${\nu }_{t,\infty }$	${\left(\frac{\kappa z}{1+\frac{\kappa z}{l_{m,\infty }^{\text{max}}}}\right)}^2\frac{\partial {\overline{u}}_{\infty }}{\partial z}$	$l_{m,\infty }^{\text{max}}=0.35\delta$ $\delta =3D$ (default)

Equations

$$f_q(x+c_q\Delta t,t+\Delta t)-f_q(x,t)={\Omega }_q(x,t)+S_q(x,t).$$

$$\underset{{\epsilon }_k}{\underbrace{\nu \overline{\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_i^{\prime }}{\partial x_j}}}}=\underset{{\mathcal{A}}_k}{\underbrace{-\overline{u}{}_j\frac{\partial k}{\partial x_j}}}\underset{{\mathcal{D}}_k}{\underbrace{-\frac{1}{\rho }\frac{\partial \overline{u_i^{\prime }{p}^{\prime }}}{\partial x_i}}}\underset{{\mathcal{T}}_k}{\underbrace{-\frac{1}{2}\frac{\partial \overline{u_j^{\prime }{u}_j^{\prime }{u}_i^{\prime }}}{\partial x_i}}}\underset{{\mathcal{V}}_k}{\underbrace{+\nu \frac{{\partial }^{2}k}{\partial x_j^2}}}\underset{{\mathcal{P}}_k}{\underbrace{-\overline{u_i^{\prime }{u}_j^{\prime }}\frac{\partial \overline{u}{}_i}{\partial x_j}}},$$

$$\underset{{\epsilon }_{k,w}}{\underbrace{{\epsilon }_k-{\epsilon }_{k,\infty }}}=\underset{{\mathcal{A}}_{k,w}}{\underbrace{{\mathcal{A}}_k-{\mathcal{A}}_{k,\infty }}}+\underset{{\mathcal{D}}_{k,w}}{\underbrace{{\mathcal{D}}_k-{\mathcal{D}}_{k,\infty }}}+\underset{{\mathcal{T}}_{k,w}}{\underbrace{{\mathcal{T}}_k-{\mathcal{T}}_{k,\infty }}}+\underset{{\mathcal{V}}_{k,w}}{\underbrace{{\mathcal{V}}_k-{\mathcal{V}}_{k,\infty }}}+\underset{{\mathcal{P}}_{k,w}}{\underbrace{{\mathcal{P}}_k-{\mathcal{P}}_{k,\infty }}}.$$

$$\underset{{\epsilon }_{k,w}}{\underbrace{{\epsilon }_k-{\epsilon }_{k,\infty }}}=\underset{{\mathcal{A}}_{k,w}}{\underbrace{{\mathcal{A}}_k-{\mathcal{A}}_{k,\infty }}}+\underset{{\mathcal{T}}_{k,w}}{\underbrace{{\mathcal{T}}_k-{\mathcal{T}}_{k,\infty }}}+\underset{{\mathcal{P}}_{k,w}}{\underbrace{{\mathcal{P}}_k-{\mathcal{P}}_{k,\infty }}}.$$

$${\mathcal{A}}_{k,w}\approx -\overline{u}\frac{\partial k}{\partial x}\approx -\overline{u}\frac{\partial k_w}{\partial x}.$$

$$-\frac{1}{2}\overline{u_i^{\prime }{u}_i^{\prime }{u}_j^{\prime }}\approx \left(\nu +\frac{{\nu }_t}{{\text{Pr}}_t}\right)\frac{\partial k}{\partial x_j}\approx \frac{{\nu }_t}{{\text{Pr}}_t}\frac{\partial k}{\partial x_j},$$

$${\mathcal{T}}_{k,w}\approx \frac{\partial }{\partial x_j}\left(\frac{{\nu }_t}{{\text{Pr}}_t}\frac{\partial k}{\partial x_j}\right)\approx \frac{\partial }{\partial y}\left(\frac{{\nu }_t}{{\text{Pr}}_t}\frac{\partial k_w}{\partial y}\right)+\frac{\partial }{\partial z}\left(\frac{{\nu }_t}{{\text{Pr}}_t}\frac{\partial k_w}{\partial z}\right).$$

$${\mathcal{P}}_{k,w}={\nu }_t\left[{\left(\frac{\partial \overline{u}}{\partial y}\right)}^2+{\left(\frac{\partial \overline{u}}{\partial z}\right)}^2\right]-{\nu }_{t,\infty }{\left(\frac{\partial {\overline{u}}_{\infty }}{\partial z}\right)}^2.$$

$${\epsilon }_{k,w}=\frac{k_w}{{\tau }_{k,w}}.$$

$${\color[rgb]{0,0,0}\overline{u}{}\dfrac{\partial k_{w}}{\partial x}=\nu_{t}\left[\left(\frac{\partial\overline{u}{}}{\partial y}\right)^{2}+\left(\frac{\partial\overline{u}{}}{\partial z}\right)^{2}\right]-\nu_{t,\infty}\left(\frac{\partial\overline{u}_{\infty}}{\partial z}\right)^{2}+\dfrac{\partial}{\partial y}\left(\frac{\nu_{t}}{\text{Pr}_{t}}\frac{\partial k_{w}}{\partial y}\right)+\dfrac{\partial}{\partial z}\left(\frac{\nu_{t}}{\text{Pr}_{t}}\frac{\partial k_{w}}{\partial z}\right)-\dfrac{k_{w}}{\tau_{k,w}}.}$$

$${\nu }_t=C_{\mu }\frac{k^2}{{\epsilon }_k},$$

$${\nu }_t=C_{\mu }k{\tau }_k.$$

$$f_p=\frac{2f_0}{1+\sqrt{1+4f_0\left(f_0-1\right){\left(\frac{\sigma }{\tilde{\sigma }}\right)}^2}}.$$

$$\sigma =\frac{k}{{\epsilon }_k}\sqrt{{\left(\frac{\partial u_i}{\partial x_j}\right)}^2}={\tau }_k\sqrt{{\left(\frac{\partial u_i}{\partial x_j}\right)}^2}\approx {\tau }_k\sqrt{{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial z}\right)}^2}.$$

$${\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial z}\right)}^2-{\left(\frac{\partial u_{\infty }}{\partial y}\right)}^2>{ϵ}_{\text{wake}}.$$

$${\epsilon }_{k,w}=\frac{k_w}{{\tau }_{k,w}}={\epsilon }_k-{\epsilon }_{k,\infty }=\frac{k}{{\tau }_k}-\frac{k_{\infty }}{{\tau }_{k,\infty }}.$$

$$\frac{k_w}{{\tau }_{k,w}}\approx \frac{k}{{\tau }_{k,\infty }}-\frac{k_{\infty }}{{\tau }_{k,\infty }}=\frac{k_w}{{\tau }_{k,\infty }},$$

$$k_w/{\tau }_{k,w}\approx k/{\tau }_k,\text{and}{\tau }_{k,w}\approx {\tau }_k(1-k_{\infty }/k)\approx {\tau }_k.$$

$${\tau }_{k,w}=\{\begin{array}{ll}{\tau }_{k,\infty }& \text{if }z\le z_{BT}=z_h-D/2,\\ {\tau }_k& \text{if }z>z_{BT}.\end{array}$$

$${\tilde{\tau }}_k=a_{\tau }x+b_{\tau },$$

$${\tau }_{k,\infty }=\frac{k_{\infty }}{l_{m,\infty }^2}{\left(\frac{\kappa }{u_{\ast }}z\right)}^3.$$

$${\tau }_{k,w}=\{\begin{array}{ll}\frac{z}{z_{BT}}\left(a_{\tau }x+b_{\tau }\right),& \text{if }0<z\le z_{BT}\\ a_{\tau }x+b_{\tau },& \text{if }z>z_{BT}.\end{array}$$

$${\mathcal{A}}_{k,w}=-\overline{u}\frac{\partial k_w}{\partial x}=-\upsilon (x){\overline{u}}_h\frac{\partial k_w}{\partial x}.$$

$${\mathcal{A}}_{k,w}\approx \underset{{\mathcal{A}}_{k,w}^a}{\underbrace{-{\overline{u}}_h\frac{\mathrm{d}K(x)}{\mathrm{d}x}g(\xi )}}+\underset{{\mathcal{A}}_{k,w}^b}{\underbrace{{\overline{u}}_h\frac{K(x)}{\sigma }\xi \frac{\mathrm{d}\sigma }{\mathrm{d}x}\frac{\mathrm{d}g}{\mathrm{d}\xi }}}.$$

$${\mathcal{A}}_{k,w}^a=-{\overline{u}}_h\frac{\mathrm{d}K(x)}{\mathrm{d}x}g(\xi )=-{\overline{u}}_h\frac{\mathrm{d}K(x)}{\mathrm{d}x}\frac{1}{K(x)}{k}_w=-{\overline{u}}_h\frac{\alpha }{(x-x_0)}{k}_w.$$

$${\mathcal{A}}_{k,w}\approx -{\overline{u}}_h\frac{\alpha }{x-x_0}{k}_w+{\mathcal{A}}_{k,w}^b.$$

$${\mathcal{T}}_{k,w}\approx -{\mathcal{A}}_{k,w}^b-a_k{\mathcal{P}}_{k,w}-b_k\Delta \overline{u}.$$

$$k_w\left({\overline{u}}_h\frac{\alpha }{x-x_0}+\frac{1}{{\tau }_{k,w}}\right)={\nu }_t\left(1-a_k\right)\left[{\left(\frac{\partial \overline{u}}{\partial y}\right)}^2+{\left(\frac{\partial \overline{u}}{\partial z}\right)}^2-\frac{{\nu }_{t,\infty }}{{\nu }_t}{\left(\frac{\partial {\overline{u}}_{\infty }}{\partial z}\right)}^2\right]-b_k\Delta \overline{u}.$$

$$\frac{1}{{\zeta }_{k,w}}:=\frac{1}{1-a_k}\left({\overline{u}}_h\frac{\alpha }{x-x_0}+\frac{1}{{\tau }_{k,w}}\right),$$

$${\gamma }_k=\frac{b_k}{1-a_k}.$$

Full text

1 Introduction

Wind farm layout design or real-time control typically rely on simplified flow models, often based on analytical wake formulations. For yield assessment and power production modelling of finite-size wind farms, the common practice is to model cumulative wake effects using analytical wake-deficit models. A key parameter in these models is the wake recovery rate, which depends heavily on the local turbulence intensity at the wake-emitting rotor position. The velocity-deficit models are thus often coupled with a wake-added turbulence model (see Niayifar & Porté-Agel, 2016; Bastankhah et al., 2021; Schmidt et al., 2023; Bay et al., 2023; Pedersen et al., 2023; Blondel, 2023). While velocity deficit models (Katic et al., 1987; Bastankhah & Porté-Agel, 2014; Keane et al., 2016; Blondel & Cathelain, 2020; Ishihara & Qian, 2018) have received significant attention, wake-added turbulence models have remained considerably less examined, despite their importance to overall flow solver accuracy.

The topic of axisymmetric wakes behind bluff bodies subjected to uniform inflows is extensively discussed in turbulence textbooks as a classic case of free-shear flows (e.g., Tennekes & Lumley, 1972; Pope, 2000). Yet, wind turbines typically operate within sheared and turbulent atmospheric boundary layers. The change in flow conditions introduces additional complexity and leads to asymmetric wake shapes, particularly in the vertical profiles of wake turbulence (Chamorro & Porté-Agel, 2008), where even the interaction of turbines with the inflow boundary layer may reduce the turbulence level close to the ground (Xie & Archer, 2015; Bastankhah & Porté-Agel, 2017).

Wake-added turbulence models for wind turbines are mostly either empirical or provide limited information on the turbulence structure in turbine wake flows. Crespo & Hernández (1996), for instance, approximated the streamwise evolution of turbulence in turbine wakes, neglecting any cross-flow spatial variations. The other model, which provides only a single representative value of the wake-added turbulence at each streamwise position, is the one developed by Frandsen (1992). It is derived based on the top-down approach of modelling wind farms (Frandsen et al., 2006). More recent empirical three-dimensional models by Ishihara & Qian (2018) or Khanjari et al. (2025) include shear effects via correction terms or fitted shape functions. Similarly, Stein & Kaltenbach (2019) proposed axisymmetric analytical models of the three velocity variances based on wind tunnel measurements behind three-bladed wind turbines. Jézéquel et al. (2024) showed that the wake-added turbulence is sensitive to atmospheric stability and can be modelled by splitting turbine-induced and meandering-induced added turbulence.

More recently, Bastankhah et al. (2024) proposed a model for axisymmetric wake flows. This model mathematically solves a simplified version of the turbulent kinetic energy (TKE) transport equation. In addition to assuming axisymmetry, which does not hold for turbine wakes in boundary layer flows, the model relies on numerical integration, making it less suitable for computationally intensive tasks such as wind farm layout optimisation.

In this work, we propose a new three-dimensional analytical model for the TKE and the streamwise normal Reynolds stress $\overline{u^{\prime}u^{\prime}}$ in wind turbine wakes. Section 2 describes the datasets generated using a large eddy simulation (LES) solver based on the lattice Boltzmann method (LBM) for a range of ground roughness and wind turbine thrust coefficients. Afterwards, we analyse the TKE transport equation budget in section 3 and derive a simplified partial differential equation (PDE) for TKE transport that is solved numerically. Using far-wake assumptions, we then propose an analytical model for the TKE written in an explicit and closed-form expression in section 4. We extend the discussion in section 5 to predictions of $\overline{u^{\prime}u^{\prime}}$. Lastly, we validate model predictions against wind tunnel measurements in section 6.

2 LES wake-added TKE and $\overline{u^{\prime}u^{\prime}}$ evolution behind a porous disk in boundary layers

Throughout this article, we employ the following notations. Three axes, $x_{i}$ for $i=1,2,3$, span the Cartesian coordinate system of the domain. Here, $x_{i}$ corresponds to the streamwise axis $x$, the lateral axis $y$, or the vertical axis $z$, respectively. Similarly, $u_{i}$ denotes the corresponding velocity components $u$, $v$ and $w$. Temporally averaged quantities read $\overline{\{\cdot\}}$ and the deviation from the mean $\{\cdot\}^{\prime}$. The Reynolds decomposition, which splits the instantaneous velocity $u_{i}$ into the mean velocity $\overline{u}_{i}$ and the velocity fluctuation $u_{i}^{\prime}$, then reads $u_{i}=\overline{u}_{i}+u_{i}^{\prime}$.

2.1 Large eddy simulation: the waLBerla-wind flow solver

To establish and validate the developed analytical model, we use data obtained from the LES of porous disks immersed in boundary layers. For this, we use waLBerla-wind, a solver based on the LBM. Unlike traditional approaches, the LBM does not aim to solve the Navier-Stokes equations for macroscopic quantities, such as pressure and velocity, but solves the lattice Boltzmann equation:

[TABLE]

Here, $f_{q}$ are the particle distribution functions (PDFs), which represent the density of particles at a given time $t$ and position $x$ with a specific discrete velocity $c_{q}$, $\Omega_{q}$ is the collision operator that relaxes the PDFs towards their equilibrium, and $S_{q}$ denotes a volumetric source term. The mass and momentum density follow from the PDFs as $\rho=\sum_{q}f_{q}$ and $\rho u=\sum_{q}f_{q}c_{q}$.

The discrete velocity set $\{c_{q}\}$ emerges during the discretisation of the continuous velocity space that adds to the temporal and spatial discretisation on a Cartesian grid. In this study, we use a three-dimensional velocity set with $27$ discrete velocities, represented by a $D3Q27$ stencil, which ensures high precision, stability, and rotational invariance (Bauer & Rüde, 2018). Furthermore, we employ the cumulant lattice Boltzmann method (Geier et al., 2015) along with a fourth-order correction (Geier et al., 2017). The QR minimum-dissipation subgrid-scale model (Verstappen, 2011; Verstappen et al., 2014) accounts for unresolved eddies. For further details about waLBerla and its extension to wind turbines and atmospheric flows, see Bauer et al. (2021) and Schottenhamml et al. (2022, 2024).

In the context of analytical wake-added turbulence modelling, mainly purely empirical models have been introduced so far. Our primary goal is to propose an approach that more realistically captures the physics of wind turbine wakes and their interaction with the ABL. To achieve this, we deliberately adopted a simplified setup to first address the truly neutral situation, similar to, e.g., see Wu & Porté-Agel (2011); Stevens et al. (2018); Lin & Porté-Agel (2019). Due to neglected Coriolis forcing, no wind veer is present in the flow. Atmospheric stratification is likewise neglected. Stratification is well known to affect wake structure and recovery, see, e.g., Hancock et al. (2014); Jézéquel et al. (2024). Our approach accounts for how wake-added turbulence is impacted by ambient turbulence, which in turn is influenced by thermal stratification. However, the more complex impact of thermal stratification on the size of atmospheric turbulent structures, wake meandering, and related phenomena is neglected in this work for simplicity.

The simulation setup is as follows: the domain has dimensions of $36D$ in length, $12D$ in width, and $5D$ in height, where $D=$0.15\text{,}\mathrm{m} is the diameter of the wind turbine. For simplicity, we model the turbine as a non-rotating, uniformly loaded actuator disk, as is often done in experimental (Aubrun et al., [2013](#bib.bib3)) and numerical studies (Wu & Porté-Agel, [2011](#bib.bib75)). We follow Stevens et al. ([2014](#bib.bib67)) and impose a modified thrust coefficient $C_{T}^{{}^{\prime}}=C_{T}/(1-a)^{2}$ with the thrust coefficients $C_{T}=\qtylist{0.4;0.6;0.8}{}$ and the wind turbine induction factor $a=1/2(1-\sqrt{1-C_{T}})$. The mean velocity at hub-height $z_{h}$ is $\overline{u}_{h}=$2.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}. Appendix A shows the Reynolds number independence of our results, by reproducing the simulation with a utility scale rotor of diameter $D=$150\text{,}\mathrm{m} and a reference wind velocity at hub height of $u_{h}=$8\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}.

Following Munters et al. (2016), we apply shifted periodic boundary conditions in the streamwise direction to prevent the development of locked large-scale turbulent structures. A free-slip boundary condition at the top of the domain represents the finite height of the atmospheric boundary layer (ABL), and a wall boundary condition (Han et al., 2021) adapted to the $D3Q27$ stencil models the ground of the domain. The wall boundary condition calculates the shear stress based on the Monin-Obukhov similarity theory. A concurrent methodology is applied, as described in Dhamankar et al. (2018). We allow the ABL in the “empty” simulation run to develop over $150$ flow-through times before starting the main simulation, which includes the wind turbine. Afterwards, both simulations run for an additional $300$ flow-through times. Data is collected during the last $150$ flow-through times to exclude wake build-up effects. The concurrent simulation provides the inflow plane for the main simulation, and a non-reflective outflow condition models the domain outlet.

For a domain height of $\delta=5D$, the ground roughness is varied over $z_{0}/\delta=$1\text{\times}{10}^{-7}$,$1\text{\times}{10}^{-6}$,$1\text{\times}{10}^{-5}$,$1\text{\times}{10}^{-4}$,$1\text{\times}{10}^{-3} in order to obtain time-averaged hub-height turbulence intensities $I_{x}=\sqrt{\overline{u^{\prime}u^{\prime}}{}}/\overline{u}_{h}$ between $5\%$ and $13\%$. This range represents typical offshore conditions, as well as moderately rough onshore conditions (Ib & Lundtang Petersen, [1989](#bib.bib33)). In addition, further simulations were performed using $z_{0}/\delta=$1\text{\times}{10}^{-7} and varying domain heights, $\delta=$3$D,$4$D,$6$D$. These cases are analysed in Section 4.4.1.

Figure 1 presents the resulting inflow velocity profiles (i.e., unaffected by the wake) and streamwise normal Reynolds stress profiles extracted four rotor diameters upstream of the actuator disk and laterally aligned with the rotor centre. In the remainder of this manuscript, the subscript $\infty$ denotes such time-averaged vertical inflow profiles. The corresponding time-averaged hub-height turbulence intensities are reported in the figure legend.

Among these simulations, we will refer to the cases with $C_{T}=0.8$ and $z_{0}/\delta=$1\text{\times}{10}^{-7} as the smooth boundary layer (SBL) case and $C_{T}=0.8$ and $z_{0}/\delta=$1\text{\times}{10}^{-4} as the rough boundary layer (RBL) case. For these cases, the turbulence intensities at hub height are approximately \qtylist5.5;9.5%. The other cases are less extensively presented in the analysis, but are included in the model calibration routines. This LES setup will be validated against laboratory experiments later in section 6.

In the remainder of this paper, the wake is divided into two regions: the near wake and the far wake. Following Bastankhah & Porté-Agel (2016), the near-wake length $x_{0}$ is estimated to range between $3$ and $5$ rotor diameters downstream of the disk, depending on the configuration. In this work, we adopt $x_{0}=4D$ as the transition between the near-wake and far-wake regions.

2.2 Velocity deficit, TKE and $\overline{u^{\prime}u^{\prime}}$ evolution

Figure 2 shows lateral and vertical profiles of velocity deficit $\Delta\overline{u}{}=\overline{u}-\overline{u}_{\infty}$, wake-added TKE $k_{w}=k-k_{\infty}$ and wake-added streamwise turbulence $\overline{u^{\prime}u^{\prime}}{}_{w}=\overline{u^{\prime}u^{\prime}}{}-\overline{u^{\prime}u^{\prime}}{}_{\infty}$ at different streamwise positions in the wake. The subscript $w$ represents the wake quantities, whereas the absence of a subscript stands for the total quantities. At each streamwise position, the quantities are normalised by their minimum or maximum values; spatial dimensions by the wake width in the vertical and lateral directions, denoted by $\sigma_{z}$ and $\sigma_{y}$, respectively.

Our simulations predict the same self-similar behaviour of the velocity deficit in the vertical and lateral directions as reported in experimental observations (Bastankhah & Porté-Agel, 2014). The wake-added TKE and $\overline{u^{\prime}u^{\prime}}{}_{w}$, however, show a transition from a bi-modal distribution in the near-wake towards a more uniform distribution in the far wake, as confirmed by the budget analysis in section 3.2. True self-similarity occurs only in the very far wake (for $x>9D$), as also observed in the experiments of Stein & Kaltenbach (2019). Furthermore, strong three-dimensionality is observed: lateral and vertical profiles differ significantly. Consequently, the assumption of self-similarity is not deemed justified for the development of an analytical model for wake-added turbulence.

Generally, the wake-added TKE and $\overline{u^{\prime}u^{\prime}}{}_{w}$ show very similar trends. Under the hypotheses we will introduce in section 5.2, both of them differ only in the presence of the pressure‐strain correlation term in the $\overline{u^{\prime}u^{\prime}}{}_{w}$ transport equation (see Lumley, 1975, among others). Hence, we will treat $\overline{u^{\prime}u^{\prime}}{}_{w}$ as an extension of the wake-added TKE models later in this study.

3 A simplified TKE transport equation

We begin by simplifying the wake-added TKE transport equation to derive a relatively simple PDE, which can be solved numerically using a one-dimensional marching-forward scheme in the streamwise direction. In the next section, this is further simplified in the far wake to develop an explicit analytical model for wake-added turbulence, eliminating the need for numerical schemes and reducing the computational time.

3.1 TKE transport equation

For an incompressible fluid with no external forces, the time-averaged TKE transport equation is (Pope, 2000)

[TABLE]

where $k=1/2\,\overline{u_{i}^{\prime}u_{i}^{\prime}}$. $\mathcal{A}_{k}$ is the advection by the mean flow, $\mathcal{D}_{k}$ is the transport by turbulent pressure fluctuations, $\mathcal{T}_{k}$ is the transport by turbulent velocity fluctuations, $\mathcal{V}_{k}$ is the viscous diffusion, where $\nu$ is the dynamic viscosity, $\mathcal{P}_{k}$ is the shear production, and $\mathcal{\varepsilon}_{k}$ is the TKE dissipation. The regularised cumulant LES model that we use is a mixture of explicit/implicit LES, as discussed in Gehrke & Rung (2022); Spinelli et al. (2023). In such a context, the TKE dissipation is a combination of numerical dissipation, the implicit part, and the contribution of the subgrid-scale model (i.e., the energy transfer from the filtered velocity field to the sub-grid scale motions), which is the explicit part. We thus cannot compute the dissipation directly from the subgrid-scale model, as in Klemmer & Howland (2024). Instead, the dissipation is taken as the residual of equation (2), as in Heinze et al. (2015) or Colombié et al. (2021). We also neglect the turbulent diffusion due to sub-grid stresses as their impact on the budget analysis is expected to be small (Klemmer & Howland, 2024). Subtracting equation (2), the total TKE transport equation, from its counterpart based on the inflow conditions yields the wake-added transport equation:

[TABLE]

3.2 Budget of the wake-added TKE transport equation

Figure 3 shows the budgets of the wake-added TKE transport equation. The four dominant terms are the production $\mathcal{P}_{k,w}$, the advection $\mathcal{A}_{k,w}$, the transport $\mathcal{T}_{k,w}$ and the dissipation $\mathcal{\varepsilon}_{k,w}$. The two remaining terms, the viscous diffusion $\mathcal{V}_{k,w}$ and the transport by turbulent pressure fluctuations $\mathcal{D}_{k,w}$ are deemed negligible: their maximum value is at most one fourth of the predominant terms and mostly below 10% in the wake region.

The production term in equation (2) directly depends on the velocity gradients.

The velocity gradients are high in the tip regions, leading to large values of production. However, below the hub height, the wake-generated gradients are typically negative, which counteracts the positive inflow gradient, leading to low or even negative added-turbulence production. Note that negative values of $\mathcal{P}_{k,w}$ do not indicate a sink of TKE but rather a lower production compared to the upstream region. The turbulent transport term acts as a diffusion mechanism, redistributing TKE from high- to lower-magnitude regions. We can also observe that the TKE dissipation exhibits a bimodal shape, which evolves with streamwise distance into a more uniform shape. Lastly, it is worth noting that in the wake centre, the advection term initially acts as a sink in the near wake and later becomes a source in the far wake, as it advects the high levels of TKE from the near wake to the far wake.

3.3 Term-by-term analysis

According to the analysis of the wake-added TKE budget in the previous section, the simplified wake-added TKE equation consists of four terms, each of which we will now discuss and simplify for wind turbine wakes:

[TABLE]

In the wake, the streamwise component of the advection term, $-\overline{u}{}\,\partial k/\partial x$, dominates over the other components, as shown in appendix B (figure 24). Considering that $\partial k_{\infty}/\partial x=0$, the wake-added advection simplifies to

[TABLE]

Similar to the advection term, one can argue that the turbulent transport outside of the wake, $\mathcal{T}_{k,\infty}$, is negligible, resulting in ${\mathcal{T}_{k,w}}{}\approx{\mathcal{T}_{k}}{}=-1/2\,\partial(\overline{u_{i}^{\prime}u_{i}^{\prime}u_{j}^{\prime}})/\partial x_{j}$. We apply the gradient diffusion hypothesis (see Harlow & Hirt, 1969; Jones & Launder, 1972; Pope, 2000) on the turbulent transport $\mathcal{T}_{k}$, and thus, approximate the triple velocity correlations in terms of the turbulent Prandtl number $\text{Pr}_{t}$ (assumed to be one) and the eddy viscosity $\nu_{t}$, defined later in section 3.4, as follows:

[TABLE]

where $\nu$ is the kinematic viscosity of the fluid. Using equation (6), we can model the wake-added turbulent transport term as

[TABLE]

The validity of these simplifications is further assessed in figure 25 in appendix B.

For the wake-added production term ${\mathcal{P}_{k,w}}{}={\mathcal{P}_{k}}{}-{\mathcal{P}_{k,\infty}}{}$, we evaluate all Reynolds stresses using Boussinesq’s eddy-viscosity hypothesis, $\overline{u_{i}^{\prime}u_{j}^{\prime}}=2/3k\delta_{ij}-2\nu_{t}S_{ij}$, defined in terms of the Kronecker delta $\delta_{ij}$ and the rate-of-strain tensor $S_{ij}=\frac{1}{2}\left(\partial\overline{u}{}_{i}/\partial x_{j}+\partial\overline{u}{}_{j}/\partial x_{i}\right)$. Generally, only gradients of the streamwise mean velocities contribute significantly to the TKE production. Moreover, the terms $\overline{u^{\prime}u^{\prime}}{}\,\partial\overline{u}{}/\partial x$, $\partial\overline{u}_{\infty}{}/\partial x$ and $\partial\overline{u}_{\infty}{}/\partial y$ are negligible. The validity of these assumptions is again discussed in appendix B (figure 26). The wake-added production term is therefore simplified to

[TABLE]

Lastly, we derive the wake-induced dissipation term based on dimensional analysis in terms of wake-added TKE and a turbulent dissipation time scale $\tau_{k,w}$. The time scale is defined as $\tau_{k,w}\propto 1/\omega$, i.e., proportional to the inverse of the specific dissipation rate used in the classical $k$–$\omega$ turbulence model developed in Wilcox (2008). This choice is consistent with the $k$–$\tau$ turbulence modelling framework introduced in Speziale et al. (1992). The dissipation takes the form

[TABLE]

An alternative approach would be to model the dissipation as: ${\mathcal{\varepsilon}_{k,w}}=ck_{w}^{3/2}/{l_{m}}$ (Pope, 2000), where $l_{m}$ is the mixing length scale and $c$ is a model constant . However, we avoid this formulation for two main reasons: first, we prefer working with a linear equation for $k_{w}$, as it later greatly facilitates the derivation of an analytical model. Second, our LES data indicate that $\tau_{k,w}$ takes a simpler three-dimensional form in the wake (shown later in section 3.5) than a mixing-length-based model, where both $l_{m}$ and $c$ can potentially vary in the wake, particularly for non-equilibrium turbulent flows (Bastankhah et al., 2024). Appendix C supports these hypotheses, showing comparisons of the turbulent time scale and turbulent length scale behaviours within wakes at different thrust coefficients and ground roughness.

Combining equations (5)-(9) results in a simplified PDE based on the eddy viscosities $\nu_{t}$ and $\nu_{t,\infty}$, turbulent dissipation time scale $\tau_{k,w}$, and the mean velocity field:

[TABLE]

The mean velocity field can be taken directly from experimental data, numerical simulations, or an analytical wake model. In the latter case, it is worth noting that since solving equation (10) includes a forward-marching scheme in the streamwise direction, any inaccuracy in near-wake modelling may propagate and affect far-wake predictions. Therefore, using an engineering model that provides realistic near-wake predictions, such as those proposed by Blondel & Cathelain (2020) and Schreiber et al. (2020), is recommended. The dependence of far-wake predictions on near-wake properties is mitigated in Section 4, where we develop an analytical TKE model for the far wake.

Assuming a given velocity field and inflow conditions, equation (10) is closed by estimating the eddy viscosity $\nu_{t}$ and the dissipation time scale $\tau_{k,w}$, which are elaborated in the following sections.

3.4 Estimation of turbulent viscosity $\nu_{t}$

Based on dimensional analysis, the eddy viscosity can be expressed as the product of the turbulent kinetic energy, $k$, and a turbulent time scale, $\tau_{k}$. In the classical formulation of Launder & Spalding (1974), the eddy viscosity is written as

[TABLE]

with the model constant $C_{\mu}$. Introducing the turbulent time scale $\tau_{k}=k/\varepsilon_{k}$, this expression can be equivalently reformulated as

[TABLE]

As noted by van der Laan et al. (2015), this approach yields a model that is too dissipative in the near wake, where the velocity gradients are high. To address this, we adopt the method proposed by van der Laan et al. (2015) and introduce a flow-dependent coefficient $C_{\mu}^{*}=f_{p}C_{\mu}$. The scalar function $f_{p}$ models the effect of non-equilibrium conditions, and reads:

[TABLE]

Here, $\sigma$ is the so-called shear parameter:

[TABLE]

which is given by $\widetilde{\sigma}=1/\sqrt{C_{\mu}}$ outside of the wake. $f_{0}$ is a constant factor, expressed in terms of the Rotta constant: $f_{0}=C_{R}/\left(C_{R}-1\right)$. Note that an alternative approach was recently proposed in Klemmer & Howland (2025), consisting in applying the inflow eddy viscosity in the near-wake region, not used here. Figure 4 shows the performance of this approach, using $C_{\mu}=0.06$ and $C_{R}=2.5$. For the LES data, the eddy viscosity is computed from $\nu_{t}={\mathcal{P}_{k}}/\left({\left({\partial\overline{u}{}}/{\partial y}\right)^{2}+\left({\partial\overline{u}{}}/{\partial z}\right)^{2}}\right)$.

The figure shows that the overall trends are captured by the model. At $x/D=4$, the model predicts a total eddy viscosity lower than $\nu_{t,\infty}$ in the wake region. In the lateral direction, at larger streamwise distances, the model tends to underpredict the eddy viscosity. In the vertical direction, the overall agreement is satisfactory, despite large oscillations in the LES data in regions with small velocity gradients, where a near-zero denominator amplifies the uncertainty in estimating the turbulent eddy viscosity.

To prevent the generation of spurious TKE in the undisturbed flow, the eddy viscosity outside the wake region is constrained to its ambient value. To this end, we introduce a wake sensor based on velocity gradients, defined as

[TABLE]

Outside the region identified by this criterion, we enforce $\nu_{t}=\nu_{t,\infty}$. In practice, a threshold value of $\epsilon_{\text{wake}}=10^{-4}$ is used.

3.5 Estimation of dissipation time scale $\tau_{k,w}$

To model the wake-added dissipation time scale, $\tau_{k,w}$, we start with the definition of the wake-added dissipation ${\mathcal{\varepsilon}_{k,w}}{}$,

[TABLE]

Quantifying $\tau_{k,w}$ directly based on the LES data is challenging since near the wake edges, where both $k_{w}$ and $\varepsilon_{w}$ are small, computing $\tau_{k,w}=k_{w}/\varepsilon_{w}$ is subject to numerical errors. Instead, we use the LES data to estimate $\tau_{k}$ and $\tau_{k,\infty}$ and assess which one is more representative of $\tau_{k,w}$ at different locations. Figure 5 shows vertical and lateral variations of $\tau_{k}$ obtained from the LES data at two different streamwise locations, where $\tau_{k,\infty}$ is also shown as a reference. Additional results are presented in Appendix appendix C, where the evolution of $\tau_{k,\infty}$ is presented along lateral and vertical profiles across a range of thrust coefficients and streamwise distances. Figure 5 shows that in the lateral direction, the value of $\tau_{k}$ remains almost constant in the wake and increases to $\tau_{k,\infty}$ outside the wake.

In the vertical profiles, we have very different behaviour in each of the three regions. Below the wake, we have $\tau_{k}\approx\tau_{k,\infty}$, a behaviour particularly evident in the RBL case. Therefore, from equation (16), we have below the wake

[TABLE]

resulting in $\tau_{k,w}\approx\tau_{k,\infty}$. Inside the wake, both $k>k_{\infty}$ and $\tau_{k}<\tau_{k,\infty}$ hold; hence,

[TABLE]

Lastly, above the wake, the wake-added TKE vanishes quickly, and therefore, the choice of $\tau_{k,w}$ becomes less significant. For simplicity, we assume $\tau_{k,w}\approx\tau_{k}$ above the wake.

Using $z_{BT}$ as the height of the wind turbine bottom-tip, the dissipation time scale results as

[TABLE]

To use equation (19), we still need to model $\tau_{k,\infty}$ and $\tau_{k}$. Figure 6 shows the streamwise evolution of $\tilde{\tau}_{k,w}$ for a selected number of LES, where the tilde denotes the best fit, obtained from a least-square minimisation procedure between the model and the LES data. A linear fit seems to provide a good approximation of $\tilde{\tau}_{k,w}$ in the far wake region. Interestingly, the slopes of the linear fits seem to be independent of the thrust coefficient and turbulence intensities, whereas the y-intercepts differ slightly depending on the flow and actuator disk conditions.

Therefore, $\tilde{\tau}_{k}$ is assumed to change linearly with $x$ within the wake region according to

[TABLE]

where $a_{\tau}$ and $b_{\tau}$ are constants. Numerical experiments indicate that treating $\tau_{k,w}$ as independent of the thrust coefficient yields satisfactory results while preserving the simplicity of the formulation. In the remainder of this article, we use $a_{\tau}=0.33$ and $b_{\tau}=8.1-33.4I_{x}$, where $I_{x}$ is the streamwise turbulence intensity at hub height outside the wake. The validity of this approximation has been assessed over the following ranges of thrust coefficients and turbulence intensities: $0.4\leq C_{T}\leq 0.8$ and $0.05\lesssim I_{x}\lesssim 0.12$.

For the boundary-layer inflow, the two dominant terms in the TKE transport equation are expected to be the turbulent production and dissipation. Equating these two terms at the inflow gives $\tau_{k,\infty}=k_{\infty}/(\nu_{t,\infty}(\partial\overline{u}_{\infty}/\partial z)^{2})$. Assuming a logarithmic velocity profile, $\overline{u}_{\infty}=(u_{*}/\kappa)\textrm{log}(z/z_{0})$, where $u_{*}$ is the friction velocity, $\tau_{k,\infty}$ takes the form of

[TABLE]

Equation (21) can be simplified by assuming that, at low heights ($z<z_{BT}$), variations of $k_{\infty}$ with $z$ are small and that $l_{m,\infty}=\kappa z$. This implies $\tau_{k,\infty}\propto z$. For practical applications, where $k_{\infty}(z)$ is unknown, and to ensure a smooth transition at $z_{BT}=z_{h}-D/2$ for the values of $\tau_{k,w}$ given in equation (19), a linear dependence of $\tau_{k,w}$ on $z$ for $0<z<z_{BT}$ can be defined as:

[TABLE]

3.6 Comparison of $k-(\tau)$ model with LES data and $k-(l)$ model

Before making further assumptions and simplifications to develop a fully closed analytical model, we evaluate the performance of the simplified TKE equation (10) against the LES data. This approach will be referred to as the $k-(\tau)$ model, $k$ referring to the resolved TKE and the parentheses indicating the analytically prescribed closure variable. Our aim here is solely to test the validity of the assumptions made in this section for deriving the simplified TKE equation, rather than to construct a complete model. Therefore, the LES data provide the velocity distribution and inflow eddy viscosity, while equation (22) provides a linear fit of the dissipation time scale; all quantities are then used in equation (10) to compute the wake-added TKE. The advection term of the simplified PDE is discretized using a forward finite difference scheme, while the radial derivatives of velocity and wake-added TKE are computed using a central finite difference scheme. The wake-added TKE is initialised to zero at $x=0$, and the second-order derivatives at the $y-$ and $z-$boundaries are estimated using a one-sided first-order scheme.

Figure 7 compares the modelled wake-added TKE based on equation (10) and the LES wake-added TKE. It shows that, for both inflow conditions, the simplified numerical solution reproduces the LES trends in the entire wake, including the near-wake region. The overall wake shape and the three-dimensional structure of the wake-added TKE are well captured.

For a more quantitative analysis, lateral and vertical TKE profiles predicted by the proposed $k-(\tau)$ model are compared with LES data, and the $k-(l)$ model proposed in Klemmer & Howland (2025). Here, $(l)$ indicates the analytically prescribed closure length scale used in  Klemmer & Howland (2025). The parameters suggested in Klemmer & Howland (2025) are adopted for comparison, whereas the $k-(\tau)$ model uses values directly from the LES. In both cases, figure 8, the agreement between the $k-(\tau)$ model and the LES data is very good, with only an under-estimation of the TKE near the wake centre at $x/D=6$ and $x/D=8$. Despite being calibrated on a different dataset, the $k-(l)$ model also performs well. It captures the relative trends, but tends to overestimate the TKE, in particular in the RBL case. As a conclusion, this section shows that equation (10) is a good approximation of the LES and that an analytical model can be built upon it.

4 An analytical model for the wake-added TKE

To develop an analytical solution for equation (10), we use far-wake assumptions to simplify the advection and transport terms further and obtain an explicit algebraic formulation of $k_{w}$ rather than a PDE. Besides the simplified formulation, solving an algebraic relation for the far-wake region has the advantage that inaccuracies in the near wake are not propagated downstream and do not negatively affect far-wake predictions.

4.1 Advection and transport terms simplification

As shown in Appendix B (figure 24), the advection term can be reduced to its streamwise component. Furthermore, we write the advection velocity as a streamwise varying velocity given by $\upsilon(x)\overline{u}_{h}$, so:

[TABLE]

Later in Section 4.3, $\upsilon(x)$ will be absorbed in a single, streamwise dependent constant $\mathcal{Z}_{k,w}(x)$, so hereafter $\upsilon(x)$ is dropped from the advection term for brevity. To simplify the advection term (and only for this term), we assume self-similarity of the TKE. While it does not hold in the near and intermediate wake regions, self-similarity works reasonably well in the far wake (see figure 2), where the advection term is prominent (see figure 3). In the context of self-similarity, we express the wake-added TKE as the product of the maximum value of the wake-added TKE, $K(x)$, and a shape function $g(\xi)$, i.e., $k_{w}=K(x)g(\xi)$. Here, $\xi=r/\sigma$ is the radial position $r$ normalised by the wake width $\sigma$ assuming wake axisymmetry, $\sigma=\sigma_{y}=\sigma_{z}$. With this, we obtain

[TABLE]

Assuming that the maximum wake-added TKE evolves as $K(x)=a_{K}(x-x_{0})^{\alpha}$ for a virtual origin $x_{0}$ and constants $a_{K}$ and $\alpha$, we write the first term of equation (24), $\mathcal{A}_{k,w}^{a}$, as

[TABLE]

Using equations (24) and (25) the simplified advection term, then, reads

[TABLE]

As will be shown hereafter, $\mathcal{A}_{k,w}^{b}$ will not be explicitly modelled and will be absorbed in the transport term.

Figure 9 compares the modelled advection term, derived using the above discussion, with the original term obtained from the LES data.

The self-similarity-based approach provides a reasonable estimation of the advection term in the very far wake, e.g., $x=10D$, but performs less satisfactorily at shorter downwind distances, e.g., $x=7D$, due to the breakdown of the self-similarity assumption closer to the turbine.

This assumption is considered acceptable because the contribution of the advection term to the budget is primarily significant in the far wake (i.e., for $x/D>7$; see figure 3) and is much less important in the near- and intermediate-wake regions. Improving the formulation of the advection term remains an area for future work, as equation (26) relies on strong underlying assumptions.

Next, we simplify the transport term by analysing its behaviour across three regions of the wake cross-section: the wake core (central region), the intermediate region where shear is maximised, and the wake boundary (edge of the wake).

Figure 10 shows the wake-added production $\mathcal{P}_{k,w}$, dissipation $\mathcal{\varepsilon}_{k,w}$, transport $\mathcal{T}_{k,w}$ and advection $\mathcal{A}_{k,w}$ at $x/D=8$. The transport and advection show similar behaviour at the wake edges. In this region, dissipation cancels production, i.e., ${\mathcal{\varepsilon}_{k,w}}{}\approx-{\mathcal{P}_{k,w}}{}$, leaving advection to cancel transport, as suggested by Tennekes & Lumley (1972). According to figure 9, $\mathcal{A}_{k,w}^{a}$ is small in this region and we therefore assume that $\mathcal{A}_{k,w}^{b}$ cancels the transport, i.e., ${\mathcal{T}_{k,w}}{}\approx-\mathcal{A}_{k,w}^{b}$.

Within the intermediate region, where wake shear reaches its maximum, turbulence production is the dominant term. As shown in figure 10, the transport term follows the distribution of production by acting as a sink of turbulence and transferring the generated turbulence to the neighbouring regions, i.e., wake centre and wake edge. We therefore assume that, in the intermediate region, the transport term can be modelled as ${\mathcal{T}_{k,w}}{}\approx-a_{k}{\mathcal{P}_{k,w}}{}$, where $a_{k}$ is a positive constant. Finally, in the wake core region, the transport term is the primary source of TKE, since the production term vanishes as $\partial\overline{u}/\partial r$ approaches zero. We therefore model the transport term as being proportional to the velocity deficit, which is maximum in this region, i.e., ${\mathcal{T}_{k,w}}{}\approx-b_{k}\Delta\overline{u}$, where $b_{k}$ is a positive constant.

To model the transport term across all three regions, we represent it as the sum of the modelled contributions from each region. This linear superposition is justified by the fact that the contribution of each modelled term becomes small outside its respective region. For example, $\Delta\overline{u}$ diminishes outside the wake core, production contributes mainly within the intermediate region, and $\mathcal{A}_{k,w}^{b}$ is close to zero except near the wake edges. The final equation describing transport across all regions, therefore, reads

[TABLE]

Substituting the advection and the diffusion terms in equation (10) with those derived in this section yields

[TABLE]

4.2 Eddy viscosity ratios simplification

At this point, equation (28) is an algebraic expression for the added TKE that can be solved to find $k_{w}$. However, it requires the calibration of six different parameters, which is cumbersome and increases the risk of overfitting the model. Experimental evidence (Bastankhah et al., 2024) and numerical studies (Iungo et al., 2017) have shown that $\nu_{t}$ approaches its asymptotic value $\nu_{t,\infty}$ in the far wake. For convenience, we therefore replace $\nu_{t}$ in equation (28) with $\nu_{t,\infty}$ in the remainder of this article. It is worth noting that this assumption was not justifiable in the previous section, where equation (10) was solved for the entire wake. Since we develop an algebraic equation applicable only to the far wake, the assumption is expected not to introduce substantial errors. We emphasise that this assumption is not mandatory in order to derive an algebraic model from equation (28) but allows to minimise the number of model parameters.

In the following, the $k-(\tau)$ model is used to evaluate the impact of this hypothesis by enforcing $\nu_{t}=\nu_{t,\infty}$ in the far-wake region. In the near-wake region, the model is used as previously described in order to provide a realistic initial condition to the simplified model and focus on the far wake. The near-wake length $x_{0}$ is computed as introduced in Bastankhah & Porté-Agel (2016). Downstream of $x_{0}$, we consider $\nu_{t}=\nu_{t,\infty}$. Results are shown in figure 11. The impact on the results is considered acceptable for an analytical model, despite a noticeable overestimation of the wake-added  TKE in the RBL case.

4.3 Grouping the advection and transport terms

A final simplification consists in grouping the advection term $\mathcal{A}_{k,w}^{a}\approx\overline{u}_{h}\alpha/(x-x_{0})$, the term $(1-a_{k})$, and the inverse of the dissipation time scale $1/\tau_{k,w}$ on the left-hand side of equation (28) as

[TABLE]

and the parameters $b_{k}$ and $(1-a_{k})$ on the right-hand side of equation (28) as

[TABLE]

The streamwise evolutions of $\zeta_{k,w}$ and $\gamma_{k}$ are currently left undetermined, and will be studied in section 4.4. Below the wake, where $z\leq z_{BT}$, we apply a linear decay on $\zeta_{k,w}$ towards zero at the wall, as we did for the dissipation time scale. $\zeta_{k,w}$ thus consist of a streamwise component, $\mathcal{Z}_{k,w}(x)$, that is damped near the ground:

[TABLE]

Using the term groupings introduced in this section and taking the limit $\nu_{t}\to\nu_{t,\infty}$, equation (28) can be further simplified to develop the final form of the analytical model:

[TABLE]

$\zeta_{k,w}$, $\gamma_{k}$ and $\nu_{t,\infty}$ are the model’s parameters and will be defined in the next sections.

4.4 Model calibration

4.4.1 Inflow eddy viscosity

A common choice for the inflow eddy viscosity is $\nu_{t,\infty}=l_{m,\infty}^{2}\times(\partial\overline{u}_{\infty}/\partial z)$, where $l_{m,\infty}$ is the mixing length in the ABL, and it is estimated based on the model proposed by Blackadar (1962). This leads to

[TABLE]

where $l_{m,\infty}^{\text{max}}$ is the maximum mixing length in the atmosphere.

Consistent with previous studies (Apsley & Castro, 1997; Tomas et al., 2011; Rodier et al., 2017), we set $l_{m,\infty}^{\text{max}}=0.35\delta$, where the boundary-layer thickness $\delta$ is assumed to be comparable to the domain height. Note that because the simulation does not include a capping inversion, the boundary layer grows with less resistance and is expected to cover the entire vertical extent of the domain.

Figure 12 (left) shows that the model provides a reasonable estimate of the mixing length for different values of $I_{x}$ and shear at $\delta=5D$.

As shown in figure 12 (right), the approximation $l_{m,\infty}^{\text{max}}=0.35\delta$ provides consistent results for domain heights ranging from $\delta=3D$ to $\delta=6D$ at a fixed roughness length. The figure also indicates that, within the rotor region, the length scale profiles become independent of the domain height once the condition $\delta>4D$ is satisfied.

4.4.2 Modelling $\zeta_{k,w}$ and $\gamma_{k}$

The evolution of $\zeta_{k,w}$ and, subsequently, $\gamma_{k}$, is computed with a best-fit procedure between the LES data and the model at several streamwise locations downstream of the turbine. The best fit minimises the root mean square error between the model predictions and the LES results over each $y$-$z$ plane. As shown in figure 13, the parameter $\mathcal{Z}_{k,w}(x)=\text{max}(\zeta_{k,w}(x,z))$ is deemed independent of both the inflow turbulence intensity and the turbine thrust coefficient, and $\mathcal{Z}_{k,w}/(D/u_{h})=0.475\,x/D-0.4$ provides a good fit for the data in the far wake region.

The evolution of $\gamma_{k}$, shown in figure 13, exhibits a strong dependence on the thrust coefficient and a weaker dependence on the inflow turbulence intensity.

Figure 13 shows that $\gamma_{k}$ increases in the near-wake until it reaches an asymptotic value in the far wake region. To keep the formulation simple and reduce the number of tuning parameters, streamwise variations of $\gamma_{k}$ are neglected. Due to the predominance of the production term in the region where $\gamma$ is not constant, more advanced formulations that were tested (but not shown here) yielded only marginal improvements in the model’s predictive accuracy. Additionally, $\gamma_{k}$ is assumed independent of the inflow turbulence intensity. Based on the model calibration, the following simple empirical relationship is proposed:

[TABLE]

4.5 Comparison to LES data

Figure 14 compares the wake-added TKE predicted using the proposed model, equation (32), and the LES data. Here, the velocity fields from the LES simulations are used as input for the analytical model.

The model captures the wake shape correctly, showing a horseshoe shape in the near wake that smoothly transitions to a more uniform TKE distribution in the far wake. While the agreement in the far wake is good, we observe an overestimation of the TKE in the near-wake at {\color[rgb]{0,0,0}x/D=2} and {\color[rgb]{0,0,0}x/D=4}, particularly near the top tip of the wake. This is to be expected, as the assumptions for this model are only valid in the far wake. A quantitative comparison is shown in figure 15, where the lateral and vertical profiles of the added TKE are plotted. The analytical model tends to overpredict the TKE near the wake edges in the near wake, i.e., at $x/D\leq 4$. Further downstream, the agreement in both lateral and vertical directions is good. Small deviations are observed, but the overall trends are well captured.

5 Extension to Streamwise Reynolds normal stress $\overline{u^{\prime}u^{\prime}}$

Analytical wake velocity deficit models used in the wind energy industry (e.g., Bastankhah & Porté-Agel, 2014; Peña et al., 2016; Blondel & Cathelain, 2020, among others) typically rely on the streamwise turbulence intensity, instead of TKE, to estimate the wake expansion rate. This section extends the work of sections 3 and 4 to the wake-added normal turbulent stress, $\overline{u^{\prime}u^{\prime}}{}_{w}$.

5.1 Budget of the streamwise turbulent stress

The transport equation for the normal turbulent stress, under the same hypotheses as for the TKE, is (Pope, 2000; Stull, 1988):

[TABLE]

The normal turbulent stress transport equation has the same structure as the TKE transport equation (2), but with an additional term, the pressure-strain correlation term $\mathcal{S}_{u^{\prime}u^{\prime}}$. Our LES flow solver may struggle to assess this term accurately as the SGS term may have a non-negligible contribution to it. Thus, two terms are assumed to be unknown in equation (35): the dissipation term $\mathcal{\varepsilon}_{u^{\prime}u^{\prime}}$ (due to the explicit/implicit nature of our LBM LES solver) and the pressure-strain correlation $\mathcal{S}_{u^{\prime}u^{\prime}}$ (due to the contribution of SGS terms not captured by LES). To determine all terms in equation (35) from our LES data, we model the dissipation term and take the pressure-strain correlation as the residual of equation (35). The dissipation tensor in anisotropic turbulence and the relation between different terms such as $\varepsilon_{u^{\prime}u^{\prime}}$ with $\varepsilon_{k}$ is an active area of research (See Perot & Natu, 2004; Gerolymos & Vallet, 2016, among others). A common approach (Hanjalić & Launder, 1976) is to blend an isotropic assumption, $\varepsilon_{u^{\prime}u^{\prime}}=2/3\varepsilon_{k}$, with Rotta’s approach (Rotta, 1951), where the anisotropy of the dissipation tensor is assumed to be aligned with the anisotropy of the Reynolds stresses, i.e., $\varepsilon_{u^{\prime}u^{\prime}}=\varepsilon_{k}\overline{u^{\prime}u^{\prime}}{}/k$. Rotta’s approach is considered valid only at low Reynolds numbers, such as close to walls. Let $f$ represent the blending parameter, which can take values between [math] and $1$. The dissipation writes

[TABLE]

where $\varepsilon_{k}$ is obtained from the TKE budget analysis discussed in section 3.1. We consider the two limiting cases, $f=0$ and $f=1$, and, therefore, the enclosed range of possible solutions for $\varepsilon_{u^{\prime}u^{\prime}}$ and $\mathcal{S}_{u^{\prime}u^{\prime}}$. Figure 16 shows the resulting budget at different streamwise positions.

The pressure-strain correlation $\mathcal{S}_{u^{\prime}u^{\prime}}$ is part of the dominant terms. Since its role is to redistribute turbulence between the Reynolds stress components, it is a sink for the normal stress in most regions of the wake, except near the ground. Because the shear production of the other normal stresses $\overline{v^{\prime}v^{\prime}}$ and $\overline{w^{\prime}w^{\prime}}$ is negligible, they are primarily sustained by the redistribution of energy from the normal turbulent stress $\overline{u^{\prime}u^{\prime}}{}$. Consequently, the pressure-strain correlation is sometimes referred to as a “return to isotropy” term (see Rotta, 1951; Sarkar & Speziale, 1990; Speziale et al., 1991, among others).

5.2 Pressure-strain correlation model

Under the hypothesis used in section 3 to derive the TKE model, the $\overline{u^{\prime}u^{\prime}}$ production, diffusion and dissipation terms take the same simplified form as their TKE counterparts. Only the additional pressure-strain correlation term requires further investigation. Following Pope (2000), the pressure-strain takes the following form:

[TABLE]

Equation (37) is a linear combination of two models proposed by Launder et al. (1975). The first term stems from Rotta’s model (Rotta, 1951) and accounts for the slow pressure effects. The second term is the so-called isotropisation of production (IP) (Naot et al., 1973) and represents the rapid pressure effects. The IP counteracts the production of $\overline{u^{\prime}u^{\prime}}$. Values of $s_{1}=1.8$ and $s_{2}=0.6$ are suggested in the literature (Pope, 2000).

Figure 17 shows the evolution of the two terms, formulated in terms of wake-added quantities, i.e., $\mathcal{S}_{u^{\prime}u^{\prime},w}=\mathcal{S}_{u^{\prime}u^{\prime}}-\mathcal{S}_{u^{\prime}u^{\prime},\infty}$, in the horizontal and vertical directions for the smooth boundary layer with low surface roughness. The overall behaviour is similar for the rough boundary layer (not shown here).

Both terms seem to vary similarly in the lateral and vertical directions, but the rapid term dominates roughly by a factor of two. Based on these observations, only the rapid term is retained in the following, with a modified value of the $s_{2}$ constant, accounting for the absorption of the slow term into the rapid term. Assuming that the production of the vertical and lateral stresses is negligible, i.e., $\mathcal{P}_{k,w}\approx\mathcal{P}_{u^{\prime}u^{\prime},w}$, the pressure-strain correlation term takes the simplified form

[TABLE]

where $s_{3}$ is a constant.

5.3 Simplified PDE and analytical model

Using the model equations for advection, transport, production, dissipation and pressure-strain correlation, the simplified PDE for the normal Reynolds shear stress reads:

[TABLE]

As in the $k-(\tau)$ model, it is possible to numerically solve equation 39. However, for the sake of brevity, we focus on the derivation of a closed-form expression for $\overline{u^{\prime}u^{\prime}}{}$. Following the same procedure used for the TKE model in section 4.3 and integrating the newly introduced constant $s_{3}$ as well as the factor $2$ on the production term into a function $\zeta_{\overline{u^{\prime}u^{\prime}}{},w}$ and a parameter $\gamma_{\overline{u^{\prime}u^{\prime}}{}}$, the analytical model reads

[TABLE]

Therefore, the analytical models for the TKE and the streamwise turbulent stress are identical. They only differ in the chosen values for the parameters $\gamma_{\overline{u^{\prime}u^{\prime}}{}}$ and $\zeta_{\overline{u^{\prime}u^{\prime}}{},w}$. We remind here that we write:

[TABLE]

5.4 Model calibration

The same procedure as for the TKE analytical model is applied to calculate $\overline{u^{\prime}u^{\prime}}{}_{w}$. The eddy viscosity in the inflow is computed using equation 33, with $l_{m,\infty}^{\text{max}}=0.35\delta$. Figure 18 (left) shows the evolution of $\mathcal{Z}_{\overline{u^{\prime}u^{\prime}}{},w}$ with respect to the streamwise distance. Once again, $\mathcal{Z}_{\overline{u^{\prime}u^{\prime}}{},w}$ appears to be largely independent of the thrust coefficient and turbulence intensity, and is represented by a linear regression of the form $\mathcal{Z}_{\overline{u^{\prime}u^{\prime}}{},w}/(D/u_{h})=0.65x/D-0.8$.

The evolution of $\gamma_{\overline{u^{\prime}u^{\prime}}{}}$ is illustrated in figure 18 (right).

As for the turbulent kinetic energy, it highlights a strong dependence on the thrust coefficient and a weaker dependence on the turbulence intensity. Based on a similar reasoning as in section 4.4.2, we model it as {\color[rgb]{0,0,0}\gamma_{\overline{u^{\prime}u^{\prime}}{}}/(\overline{u}_{h}^{2}/D)=5\times 10^{-3}\,C_{T}^{2}}.

5.5 Comparison to LES data

We compare the analytical model predictions with the LES dataset for both the smooth and the rough boundary layers in figure 19.

We note that, similarly to the previous section, the LES data was used to estimate the velocity field, which then served as input for the $\overline{u^{\prime}u^{\prime}}{}_{w}$ field. Generally, we observe a qualitatively good agreement between the LES and the model predictions. Further validation of the model against experimental data will be performed in the next section.

6 Validations and comparison to two wind-tunnel datasets

The analytical model for the normal turbulent stress is validated against data from two wind-tunnel campaigns covering inflow turbulence intensities from $I_{x}\approx 6.7\%$ to $I_{x}\approx 13.8\%$.

Note that the model predictions rely on parameters calibrated against the LES data and assumed to be universal. Consequently, the wind-tunnel datasets were strictly “unseen” during model calibration and serve solely for validation purposes. Here, to have a complete engineering modelling framework and predict velocity distributions, we use the super-Gaussian wake model (Blondel & Cathelain, 2020), although other wake velocity models (e.g., Bastankhah & Porté-Agel, 2014; Schreiber et al., 2020; Bastankhah et al., 2021; Ali et al., 2024, among others) could also be used. The parameters used in the present comparison are taken from Blondel (2023). Table 1 provides a summary of the velocity deficit and wake-added turbulence models, together with all employed model parameters. The tilde superscript indicates normalization with respect to the turbine diameter $D$, and $\Gamma(\cdot)$ is the gamma function.

6.1 Wind tunnel experiments by Bastankhah & Porté-Agel (2017)

The first test case considers the wind-tunnel data reported in Bastankhah & Porté-Agel (2017). A $0.15\text{\,}\mathrm{m}$ diameter miniature wind turbine is immersed in a turbulent boundary layer with a hub-height turbulent intensity of approximately $6.7\%$. The hub centre is located at $z_{h}=$0.125\text{,}\mathrm{m}$$. In this case, the turbine operates at a high thrust coefficient of $C_{T}\approx 0.77$.

Figure 20 shows the lateral and vertical profiles of normalised velocity deficit and wake-added turbulence intensity, respectively. The figures contain laboratory measurements, the LES data based on our LBM solver, and our proposed analytical model. They also include the added turbulence intensity predictions made by Ishihara & Qian (2018), shown as the I&Q model.

The figures demonstrate that the super-Gaussian model reproduces the velocity deficit in both lateral and vertical directions with good accuracy, enabling its application within the wake-added turbulence model.

In terms of lateral wake-added turbulence intensity profiles, defined as $I_{u,add}=\sqrt{(\sigma_{u}/\overline{u})^{2}-(\sigma_{u_{\infty}}/\overline{u}_{\infty})^{2}}$ for $\sigma_{u}/\overline{u}\geq\sigma_{u_{\infty}}/\overline{u}_{\infty}$ and $I_{u,add}=-\sqrt{(\sigma_{u_{\infty}}/\overline{u}_{\infty})^{2}-(\sigma_{u}/\overline{u})^{2}}$ otherwise, the model follows the experimental trends fairly closely over the entire set of streamwise distances. The agreement between the LES simulations and the analytical model is also very good, both following the experimental trends. In the near wake, at $x/D=3$ and $x/D=4$, an overestimation of the wake-added turbulence is observed in the LES. This could be explained by the absence of a nacelle in the simplified actuator disk approach employed in the simulations. Below hub height, the I&Q model overpredicts the wake-added turbulence intensity, for all downstream distances. As shown in the vertical profiles, there is a mismatch between the measurements and the LES simulations close to the ground. Such a deviation was also observed in Dar et al. (2024). The overall behaviour of the proposed model for this first test case is highly satisfactory.

6.2 Wind tunnel experiments by Stein & Kaltenbach (2019)

The second dataset, based on the experimental campaign reported in (Stein & Kaltenbach, 2019)

, allows for a direct assessment of the influence of the inflow turbulence intensity. Two inflow conditions are examined: a SBL with $TI\approx 8.3\%$ and a RBL with $TI\approx 13.8\%$. The rotor diameter and hub centre are equal, i.e., $z_{h}=D=$0.45\text{,}\mathrm{m}$$.

As shown in figure 21, the super-Gaussian model shows good agreement with the reference data and provides reliable input for the wake-added turbulence model. In the lateral direction, figure 21, the proposed model is able to correctly reproduce the wake-added streamwise normal Reynolds stress behaviour for both ABL s, highlighting its robustness. The I&Q model reproduces the SBL case very accurately, but strongly overestimates the streamwise normal Reynolds stress in the RBL case. In the vertical direction, figure 21, similar observations are drawn. In the near-wake and below the top-tip region, a clear overestimation of the wake-added $\overline{u^{\prime}u^{\prime}}{}$ is observed in the RBL case. In the rest of the wake, the proposed approach slightly overestimates the turbulence at the top tip position, but accurately reproduces the rest of the profile. The I&Q model underpredicts the wake-added streamwise normal Reynolds stress below the top-tip region in the far wake, $x/D\geq 5$, in the RBL case, but accurately reproduces the SBL case. The LES simulation is very close to the measured data, even in the near-wake region.

7 Conclusion and future work

In this work, we derived an analytical model to predict the TKE and the normal turbulent stress and validated the results against large eddy simulations and wind tunnel measurements. It is grounded in the analysis of turbulence transport equation budgets, classical approximations commonly used in eddy-viscosity Reynolds-averaged Navier-Stokes (RANS) closures, and additional simplifications leading to a tractable analytical formulation. The derivation achieved several intermediate results. Firstly, we formulated a simplified transport equation for wake-added TKE, which supports initial assumptions and can be used directly due to its low computational cost. From this equation, we then use far-wake assumptions to obtain an analytical model for wake-added TKE by applying a self-similarity hypothesis for the advection term and a budget-based simplification of the diffusion term. We then applied the same methodology to the $\overline{u^{\prime}u^{\prime}}{}_{w}$ transport equation, resulting in a final model structurally similar to the TKE model. This model can replace empirical turbulence models in analytical wind farm solvers. Its reliance on physical principles and its ability to incorporate inflow shear make it a powerful predictive tool. Future work will focus on extending the model to farm-scale wake predictions and the inclusion of atmospheric stability. Finally, it is worth noting that the wake recovery rate is expected to depend on the local level of turbulence rather than the inflow value (Pedersen et al., 2022). Accounting for this can be achieved in future work by coupling a wake velocity model with the turbulence model developed in this study via an iterative approach.

Appendix A Impact of the Reynolds number on the LES simulations

In the setup used throughout this work, wind-tunnel-scale simulations are performed, based on the solver validation against the experimental data from Bastankhah & Porté-Agel (2017). Here, an additional simulation at utility scale is presented, using a rotor diameter of $D=$150\text{,}\mathrm{m}$$ and a velocity at hub-height of $8\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. The other parameters are kept constant. With this setup, we want to assess the impact of the rotor diameter (more generally, of the Reynolds number) on our simulations. At the wind tunnel scale, the Reynolds number based on the rotor diameter is of roughly $2\times 10^{4}$, while at the utility scale, it reaches $8\times 10^{7}$.

Figure 22 shows that the impact of the Reynolds number on the predicted velocity is largely negligible: the profiles are almost perfectly superimposed.

Similarly, figure 23 shows the impact of the Reynolds number on the TKE and on the streamwise component of the normal Reynolds stress. The resulting profiles are essentially independent of the Reynolds number. In an experimental study with a scaled wind turbine, Chamorro et al. (2012) observed a good level of Reynolds-number independence only for slightly higher Reynolds numbers. However, as noted in (Chamorro et al., 2012), the wake in the wind tunnel might be influenced by the Reynolds-number-dependent local turbine airfoil aerodynamics, which does not occur in our actuator-disk simulations.

Appendix B Verification of the advection, diffusion and production terms treatment

This appendix assesses the validity of some of the assumptions used to simplify the TKE transport equations in section 3.3. Figure 24 shows the full advection term, the simplified advection term, which consists only of $\overline{u}\partial k_{w}/\partial x$, and the linearised term, $\overline{u}_{h}\partial k_{w}/\partial x$. Despite some numerical oscillations at the outer edges of the wake, the overall conclusion remains unaffected: using only one term is sufficient to model the advection accurately. Furthermore, linearising the advection term by replacing $\overline{u}$ with $\overline{u}_{h}$ has only a marginal impact in the far-wake.

Figure 25 compares the proposed model for the turbulent transport term and the complete term. Overall, the agreement between the chosen transport model and the complete term is less good than for the advection and production term. This discrepancy could originate, at least partially, from numerical errors introduced when computing the second-order derivative of the wake-added TKE.

The last term discussed in the appendix is the production term. Figure 26 compares the complete production term to its model based on the Boussinesq eddy viscosity hypothesis and a simplified rate-of-strain tensor. The eddy viscosities in the model is obtained directly from the LES data as $\nu_{t}=\mathcal{P}_{k}/(2S_{ij}S_{ij})$.

The model reproduces the LES data very well.

Appendix C Additional results: turbulent time scale evolution

In section 3.5, a simple model for the turbulent time scale $\tau_{k}$ in a wind turbine wake is introduced, based on the assumption that $\tau_{k}$ remains constant within the wake. The validity of this assumption is further assessed here using additional simulation cases.

Figures 27 and 28 present the evolution of $\tau_{k}$ for three different thrust coefficients $C_{T}$, plotted as a function of the lateral and vertical distances normalised by the Gaussian wake width $\sigma$. The evolution of the turbulent length scale $l=k^{3/2}/\varepsilon$ is also shown. Since $l$ is often assumed to be constant in wake modelling (see, e.g., Klemmer & Howland (2025); Bastankhah et al. (2024)), comparing its behaviour with that of $\tau_{k}$ is of particular interest. Figure 27 corresponds to the SBL case, while figure 28 focuses on the RBL case.

For the lateral profiles, all quantities are normalised by their respective maximum values, as is customary when analysing self-similar behaviour. For the vertical profiles, normalisation is performed using the mean value over the profile. This choice avoids excessive compression of the curves and improves the readability of the figures.

Both figures indicate that treating $\tau_{k}$ as constant within the wake region is a reasonable approximation, regardless of ground roughness. Moreover, the validity of this assumption appears to improve with increasing thrust coefficient.

The same figures also suggest that assuming a constant turbulent length scale $l$ within the wake is more questionable. While this approximation seems acceptable in the lateral direction, significant variations of $l$ are observed in the vertical direction, with a noticeable increase with height.