Articles | Volume 9, issue 1
Research article
18 Jan 2024
Research article |  | 18 Jan 2024

Breakdown of the velocity and turbulence in the wake of a wind turbine – Part 2: Analytical modelling

Erwan Jézéquel, Frédéric Blondel, and Valéry Masson

This work aims to develop an analytical model for the streamwise velocity and turbulence in the wake of a wind turbine where the expansion and the meandering are taken into account independently. The velocity and turbulence breakdown equations presented in the companion paper are simplified and resolved analytically, using shape functions chosen in the moving frame of reference. This methodology allows us to propose a physically based model for the added turbulence and thus to have a better interpretation of the physical phenomena at stake, in particular when it comes to wakes in a non-neutral atmosphere. Five input parameters are used: the widths (in vertical and horizontal directions) of the non-meandering wake, the standard deviation of wake meandering (in both directions) and a modified mixing length. Two calibrations for these parameters are proposed: one if the users have access to velocity time series and the other if they do not. The results are tested on a neutral and an unstable large-eddy simulation (LES) that were both computed with Meso-NH. The model shows good results for the streamwise velocity in both directions and can accurately predict modifications due to atmospheric instability. For the axial turbulence, the model misses the maximum turbulence at the top tip in the neutral case, and the proposed calibrations lead to an overestimation in the unstable case. However, the model shows encouraging behaviour as it can predict a modification of the shape function (from bimodal to unimodal) as instability and thus meandering increases.

1 Introduction

The CPU cost of classical computational fluid dynamic models is too high to deal with all the different cases needed to estimate and optimize the performances of a wind farm. Thus, so-called engineering models have been developed to estimate the power loss due to wakes at a low computational cost (e.g. Jensen1983; Larsen et al.2008; Bastankhah and Porté-Agel2014). These design tools are based on physical considerations and are often calibrated and validated against numerical results or measurements. Among these tools, analytical models are the quickest: they consist of a single formula that can be directly applied to the wind farm setup and atmospheric conditions, leading to fast results even for a whole farm. A very commonly used model is the one developed by Bastankhah and Porté-Agel (2014), who assumed an axisymmetric and self-similar Gaussian velocity deficit in the wake and solved the mass and momentum conservation equations to find a relation between the amplitude and width of the Gaussian. It can be adapted for a non-axisymmetric wake (Xie and Archer2014):


where U is the mean velocity field; U is the mean velocity upstream of the turbine; C(x) is the maximum velocity deficit; CT is the thrust coefficient; D is the turbine diameter; (x,y,z) denotes the streamwise, lateral and vertical coordinates, centred at the turbine's hub; and σy,z denotes the wake widths in the lateral and vertical directions. In the present work, the vertical and horizontal axes are centred at the hub position. Here and in the following, the Reynolds decomposition is used to write any unsteady field X(t) as a sum of a mean and a varying part: X(t)=X+X(t).

The stability of the atmospheric boundary layer (ABL) influences the wake recovery (Abkar and Porté-Agel2015), and the large-scale eddies carried in this region of the atmosphere are often associated with wake meandering, i.e. oscillations in the instantaneous wake around its mean position (Larsen et al.2008). To model the meandering, the concepts of fixed and moving frames of reference (denoted FFOR and MFOR respectively) defined in the dynamic wake meandering (DWM) model are used herein (Larsen et al.2007). The FFOR is bound to the ground: it is the frame of reference in which we want to compute the turbulence and velocity fields. In the FFOR the effects of meandering are not differentiated from the wake expansion caused by turbulent mixing. The MFOR moves with the wake centre at each time step: in this frame of reference, only the wake expansion due to turbulent mixing is represented, making the fields in this frame of reference easier to interpret. The instantaneous streamwise velocity can be changed from one frame to another according to the following relation:

(3) U MF ( x , y , z , t ) = U FF x , y + y c ( x , t ) , z + z c ( x , t ) , t ,

where subscripts MF and FF denote the velocity fields in the MFOR and FFOR respectively, and yc(x,t) and zc(x,t) are the time series of the wake centre's coordinates at the downstream position x. The concept of MFOR and FFOR can be used to write an analytical wake model for the velocity deficit as in the work of Braunbehrens and Segalini (2019):

(4) Δ U FF ( y , z ) = C 1 + σ f y σ y - 1 / 2 1 + σ f z σ z - 1 / 2 exp - y 2 2 σ y 2 + 2 σ f y 2 - z 2 2 σ z 2 + 2 σ f z 2 ,

where σfy,fz(x) denotes the standard deviations of the wake centre's coordinates in the lateral and vertical directions respectively, σy,z(x) denotes the wake widths in the MFOR, and C(x) is the maximum velocity deficit in the MFOR. Such a model allows for the independent calibration of the effects of meandering (through the variables σfy,fz) and of wake expansion due to turbulent mixing (through the variables σy,z). The former parameters are a function of atmospheric stability through lateral and vertical turbulence (Braunbehrens and Segalini2019; Du et al.2021; Brugger et al.2022), whereas the latter parameters can be a function of axial turbulence as in Eq. (1) (Fuertes et al.2018; Niayifar and Porté-Agel2016) or turbine operating conditions such as CT and atmospheric shear (Braunbehrens and Segalini2019).

For the turbulent kinetic energy (TKE), it is common to model only the maximum value of added turbulence, which can be computed with the Crespo model (Crespo and Hernandez1996) or the Frandsen model (Frandsen2007) as in the IEC 61400-1 standard. Their approach is mainly empirical and can be extended to describe the whole profile of turbulence instead of the maximum value alone (Ishihara and Qian2018). This widely used model (hereafter denoted I&Q2018) is simple, since it only requires the knowledge of the thrust coefficient and the upstream turbulence intensity, but it is totally empirical and does not account for atmospheric stability. Moreover, it has been shown that the wake in an unstable ABL dissipates faster than in a neutral ABL, even at the same level of turbulence intensity (Du et al.2021). This behaviour cannot be taken into account in the I&Q2018 model due to the limited number of inputs.

The present work aims to propose a physically based model that predicts both the mean and the variance (i.e. turbulence) of the axial velocity in the wake of a wind turbine. The advantage of basing our model on physical interpretations is that it gives more room for further improvements, as we know which assumptions were made and how it degrades the results. Moreover, the proposed model is dependent on atmospheric stability, since it influences both the velocity and the turbulence fields in the wake (see companion paper). Many models, such as the I&Q2018 model, do not take atmospheric stability into account, assuming that stable and unstable cases compensate for each other and that a calibration on neutral cases is thus sufficient. This approach is valid for monthly or yearly estimations of wind farms' performances. But some applications of the future wind industry, such as digital twins, need estimations over a day, an hour or even smaller periods. In such cases, the stability must be taken into account. Since we showed in the companion paper that stability mainly affects the wake meandering, this phenomenon must be decoupled from the wake expansion to take the ABL stability into account. To do so, the breakdowns described in the companion paper are reused and briefly explained in the following.

A field in the MFOR can be written as an unsteady translation of the same field in the FFOR through Eq. (3). To shorten this equation, the notation a(y,z)^=a(y-yc(t),z-zc(t)) for any field a is used. For the present work, it is also important to note that for any field a,

(5) a ^ = a * * f c ,

where ** denotes a 2D convolution and fc is the probability density function (PDF) of the wake centre position. In the companion paper, it was shown that the velocity (Eq. 6) and turbulence (Eq. 7) in the FFOR can be expressed as a function of their counterparts in the MFOR. This is achieved by dividing these quantities into several terms, (I) and (II) in Eq. (6) and (III) to (VII) in Eq. (7).

(6) U FF = U MF ^ ( I ) + U MF ^ ( II )

(7) k FF = U MF ^ 2 - U MF ^ 2 k m = ( III ) + k MF ^ k a = ( IV ) + 2 cov U MF ^ , U MF ^ ( V ) + U MF 2 ^ ( VI ) - U MF ^ 2 ( VII )

These terms are thoroughly described and quantified in the companion paper, where they are separated into pure terms (I, III and IV) and cross terms (II, V, VI and VII).

Term (I) is the convolution of UMF with fc. It is a pure mean velocity term: it is 0 only if the mean velocity is 0. Conversely, term (II) is a cross term because it can be equal to 0 either if there is no meandering (operator .^ has no effect) or if there is no turbulence in the MFOR (UMF=0). Term (III), also written as km in the following to be consistent with notation from Keck et al. (2013) and Conti et al. (2021), is the turbulence purely induced by meandering: in the case of a meandering steady wake, i.e. UMF=0, Eq. (7) reduces to this term only. Term (IV) is the rotor-added turbulence, which is also written as ka for consistency with other works (Conti et al.2021). It is the turbulence purely induced by the rotor: in the absence of meandering, the equation reduces to this term only, also written as ka in the following for consistency with the literature. Term (V) is the covariance of UMF^ and UMF^, term (VI) can be viewed as the varying part of the MFOR turbulence, and term (VII) is the square of term (II). It is a pure dissipation term as it is always negative. Like term (II), they are cross terms since they are equal to zero if either the turbulence in the MFOR or the meandering is null. The companion paper showed that terms (II) and (VII) are negligible in their respective equations. In the breakdown of the turbulence equation, term (V) is of lesser importance than (III) and (IV) but drives the vertical asymmetry of the turbulence profiles.

The proposed analytical model is based on the velocity and turbulence breakdowns (Eqs. 6 and 7). Similar to Eq. (4) (Braunbehrens and Segalini2019), the reasoning starts by writing the wake properties in the MFOR and the wake meandering with different parameters to take into account meandering due to atmospheric stability independently of the expansion due to turbulence mixing. It is common in wake modelling to assume that meandering can be entirely accounted for by increasing the wake expansion. However, it is a phenomenon of a different nature, and it leads to velocity and turbulence profiles of different shapes. In the present model, these phenomena are modelled separately, and it is assumed that they do not interact. This is equivalent to neglecting cross terms in Eqs. (6) and (7), which have been shown to take smaller values than pure terms in the companion paper. However, in the future, modelling these cross terms might be necessary to improve the results. The main added value of this work is to propose a new framework that can be used with different shape functions in the MFOR to propose other turbulence models. Nevertheless, two calibrations (one requiring the inflow time series and one that does not) are proposed for the model to demonstrate how it can be tuned and to test the model.

In the second section of this work, the datasets are presented: for the calibration of the model, a dataset from the MOMENTA project is used, and for the validation the neutral and unstable cases obtained from the large-eddy simulations (LESs) from the companion paper are reused. The third section presents the derivation of the model. The fourth section shows the chosen calibration methods, and the fifth section presents the corresponding results. All these results are discussed in the sixth section, followed by the conclusion.

2 The LES datasets

2.1 Description of the LES code

The analytical model developed in this work is based on LES datasets generated with the Meso-NH solver (Lac et al.2018). It is a finite-volume and finite-difference research code for ABL simulations where the Navier–Stokes equations and the energy conservation equation are resolved on an Arakawa C grid. This solver models the stability of the ABL with a buoyancy term in the momentum equation. The Coriolis force and large-scale forcing are also taken into account. The effect of the wind turbine on the surrounding flow is modelled with an actuator line method, i.e. rotating source terms in the momentum equation.

To close the set of equations, the subgrid TKE equation is resolved, and all the subgrid quantities are written as a function of this subgrid TKE, the resolved variables and a Deardorff mixing length. A grid nesting method allows for simultaneously having a vertical and horizontal mesh size of 1.5 and 0.5 m in the wake region for the two datasets and a domain large enough to compute the largest eddies of the atmosphere. The model and numerical parameters are described in more detail in the companion paper.

2.2 Simulation setup

Two different LES datasets are used in this work: the first one for creating and calibrating the model and the second one for testing the model's results. Inflow conditions of these datasets can be found in Table 1. For both datasets, only the wake mean streamwise velocity (Ux, written as Ux in the following) and the streamwise turbulence (kx=uu) are computed. The proposed model thus only deals with the streamwise velocity and turbulence.

Table 1List of LES cases.

Download Print Version | Download XLSX

The calibration dataset contains six simulations, with four different ABL stabilities and three different thrust values. The simulated turbine is 92 m in diameter and has a hub height of 80 m. The turbine's data were obtained in the context of the MOMENTA project (Jézéquel2023).

To perform such simulations, a precursor without heat flux is first simulated in a domain of 19 km × 15 km (with a horizontal resolution of 37.5 m) for 25 h to let the turbulence establish itself and to allow the system to reach a quasi-steady state. Then, a ground heat flux is applied for 4 h: 0, 30, 60 and 120 W m−2 for “neutral”, “weakly unstable”, “unstable” and “strongly unstable” cases respectively. This allows us to simulate three different levels of atmospheric instability, starting from the same neutral state. No stable case was simulated because of the induced veer (gradient of inflow wind direction) that leads to a deformed wake. The veer could have been modelled as in Abkar et al. (2018), but it would significantly complicate the present derivations. Moreover, meandering and meandering turbulence are negligible in a stably stratified ABL (see companion paper), and there is thus little interest in using the approach presented herein. Developing the model for veered cases is a challenge that is out of the scope of this work.

After these two steps, the coarsest computational domain (horizontal resolution of 37.5 m) is ready: two grid nestings are then applied to reach a resolution of 1.5 m in the most refined domain. Then, 10 min of dynamics is used to let the flow establish itself in the wake of the wind turbine, and the post-processing is performed on the following 50 min of dynamics. The data are sampled at 1.2 Hz, which is the approximate limit between resolved and subgrid TKEs for these simulations (equivalent to 4 times the mesh size).

The wind turbine rotational speed and pitch are set according to the controller's database and the calculated upstream velocity. Since all the cases are computed at a similar inflow velocity, similar values of the thrust coefficient are obtained in the simulations. To include the influence of the thrust coefficient in the model, two additional cases with a degraded thrust coefficient are also computed, with the same inflow as the neutral case. To reduce the thrust, the pitch value is increased from 0 to 3 and 4.5 respectively.

The second set of simulations, hereafter called the validation dataset, is based on the neutral and unstable cases that are described in the companion paper. The simulated turbine is a modified version of the Vestas V27: it is a three-bladed rotor with a diameter D=27 m and a hub height of 32.1 m. The simulation methodology is quite similar, as described in the paragraph above, except that one additional nesting is required to reach the targeted mesh size. In the validation dataset, the velocity is sampled at 1 Hz, and the simulations last for 80 and 40 min for the neutral and unstable cases respectively. This was due to benchmark requirements and computational limitations. A statistical convergence of our datasets is proposed in the Appendix of the companion paper. Overall, it concluded that increasing the duration of the simulation for the unstable case would improve the reliability of the simulations.

2.3 Wake tracking

For the validation simulations, the wake centre's yc(x,t) and zc(x,t) coordinates are computed at each time step and each downstream position with the constant-flux wake-tracking algorithm, which is described in the companion paper. To facilitate the wake tracking, a reference simulation is also run. It is a simulation with the same inflow and boundary conditions but without the wind turbine. The corresponding velocity field noted as Uref thus represents a developing ABL without the perturbations of a wind turbine.

Figure 1Time series of the wake centre's lateral (a) and vertical (b) coordinates with the constant-flux method and the pollutant method. Weakly unstable case at x/D=6 between 1000 and 2500 s.


Another method is proposed here to compute the unsteady wake centres in the calibration dataset. A passive scalar (similar to a pollutant) is emitted at the rotor disc with a concentration value of 1 at each time step. This new variable is only driven by the advection scheme, in accordance with the passive tracer of the DWM theory, and only marginally impairs the code's performance. By supposing that this variable follows the wake, the unsteady wake centre is deduced from the centre of mass of this pollutant at each downstream position. The results lead to a low-frequency behaviour similar to the constant-flux method used in the companion paper but with fewer outliers (see Fig. 1). Since the post-process is more straightforward and the results seem better, this method has been used for the calibration dataset.

2.4 Inflow conditions

Table 1 shows the hub height velocity, thrust coefficients and turbulence intensities at hub height for each of the cases. The directional turbulence intensities are defined as

(8) I x , y , z = k x , y , z U , hub ,

and the global turbulence intensity is defined as

(9) I = 1 3 I x 2 + I y 2 + I z 2 .

Figures 2 and 3 show the profiles of some inflow variables for the calibration and validation cases respectively. The profiles are taken 2.5 diameters upstream of the wind turbine and are averaged along the y direction (the direction transverse to the wind turbine).

Figure 2Inflow conditions for the calibration cases. (a) Mean velocity profile, (b) mean TKE profile and (c) mean kx to shear ratio profile. Solid lines: LES results; dotted lines: fit with the Monin–Obukhov law.


Figure 3Inflow conditions for the validation cases. (a) Mean velocity profile, (b) mean TKE profile and (c) mean kx to shear ratio profile. Solid lines: LES results; dotted lines: fit with the Monin–Obukhov law.


In the left panel the mean velocity is plotted. The calibration dataset (Fig. 2) was built in order to have similar hub height velocities between the cases (around 7 m s−1), whereas the validation dataset comes from simulations that reproduced the SWiFT benchmark where the hub height velocities differed. The Monin–Obukhov profiles are plotted as follows using dotted lines:

(10) U ( z ) = u * κ ln z / z 0 + ψ ( z , L MO ) ,

where κ=0.41 is the von Kármán constant (Cheng et al.2019),

(11) ψ z , L MO = - 2 ln 1 + x u / 2 - ln 1 + x u 2 / 2 + 2 arctan x u - π / 2 ,

and xu=(1-15z/LMO)0.25. Since z0 is known from the simulations (0.17 in the calibration dataset and 0.014 in the validation dataset), the profiles are found by fitting Eq. (10) on the corresponding velocity profile, with parameters u* and LMO. The results, in dotted lines, match the inflow profiles well, showing that it respects the Monin–Obukhov similarity theory around the turbine's height.

The middle panels of Figs. 2 and 3 show the inflow TKE, defined as

(12) k = 1 2 k x + k y + k z = 3 2 I U , hub 2 .

In the calibration dataset, the amount of TKE increases as the imposed heat flux increases. In the validation dataset, this is not the case since the neutral case is at a higher velocity at hub height, but the TI of the unstable case is indeed higher than that of the neutral case.

The right panels of Figs. 2 and 3 show the modified mixing length lm,* upstream of the wind turbine. This quantity is discussed and used in Sect. 3 to compute the mixing length in the MFOR. Here, the value is computed as the ratio of turbulence and shear:

(13) l m , * = k x , U z .

However, in the unstable cases, the velocity profile becomes nearly constant above a given height, leading to low values of U/z and thus very chaotic behaviour of lm,*. To have a more reliable curve, the derivative of U is resolved analytically using Eq. (10):

(14) U z = u * κ z 1 - 15 z / L MO - 0.25 ,

with LMO and u* fitted from the velocity profile. The resulting curve, in dotted lines, gives a more useful quantity in the turbulence to shear ratio, while still being of the order of magnitude of the directly computed ratio (in solid lines).

3 Model derivation

In this section, we derive an analytical model for the dominating terms of Eqs. (6) and (7). First, an analytical form is proposed for the velocity deficit in the MFOR ΔUx,MF, the turbulence in the MFOR kx,MF and the meandering distribution fc. Then, some terms are neglected, and the convolutions of Eqs. (6) and (7) are resolved analytically to get a model for the velocity deficit and added turbulence in the FFOR. To help the reader, the main variable notations and subscripts used in this section and afterwards are summarized in Table 2.

Table 2Description of the most-used notations in this part and the following.

Download Print Version | Download XLSX

3.1 Independent modelling of the wake in the MFOR and meandering

3.1.1 Wake velocity deficit in the MFOR

Based on the literature (Bastankhah and Porté-Agel2014; Xie and Archer2014), the mean velocity deficit is modelled with the long-established Gaussian velocity deficit (see Eq. 1):

(15) Δ U x , MF , am ( x , y , z ) = C ( x ) exp - y 2 2 σ y 2 ( x ) - z 2 2 σ z 2 ( x ) ,

where subscript .am stands for the analytical model; C(x) is defined in Eq. (2); and σy, σz denotes the wake widths in the MFOR. The overline is dropped because this analytical model is static. In the literature, it has been shown that double-Gaussian (Keane et al.2016) or super-Gaussian (Blondel and Cathelain2020) shapes provide more accurate results, but here the Gaussian shape allows for a straightforward computation of the convolutions in our model and is still pertinent in the far wake. It is shown later that this approximation leads to discrepancies in the near wake.

3.1.2 Wake-added turbulence in the MFOR

To model term (IV) or kx,a, one needs an analytical form for the turbulence in the MFOR kx,MF. It was first thought that it would be better to model the added turbulence in the MFOR, i.e. Δkx,MF=kx,MF-kx,, in order to separate the rotor-added turbulence Δkx,MF from the ambient turbulence. This procedure was done in the companion paper; however it leads to negative values of Δkx,MF (in particular near the ground), i.e. smaller turbulence in the wake compared to the turbulence upstream of the wind turbine. This is not compatible with a model that predicts only increased turbulence in the wake of a turbine (as is the case here or in I&Q2018), and this approach has thus been abandoned.

The derivation of a model for kx,MF is not as straightforward as for ΔUMF because turbulence comes from the unsteadiness of the flow, whereas an analytical model is by definition steady. In the DWM, the Madsen formulation (Madsen et al.2010) is used to scale the velocity profile with an empirical function of the wake-generated shear. One could also assume self-similarity of the Δkx,MF function and try to derive a model as was done for the velocity in Bastankhah and Porté-Agel (2014). The main issue here is that an analytical form of the model is needed in the FFOR; i.e. the convolution of fc,am with the chosen shape function for Δkx,MF,am must have an analytical solution, which is not trivial for the aforementioned models.

It is proposed here to assume that the turbulence in the MFOR is solely driven by wake-generated shear, as in Madsen et al. (2010). To relate the turbulence in the MFOR to mean gradients, two models for the velocity scale u0 are combined. In the first, it is assumed to be proportional to the square root of the TKE (Pope2000). However in the present work, the 3D TKE is not computed, so it is replaced with the axial turbulence kx:

(16) u 0 = C μ 1 / 4 k x 1 / 2 ,

where Cμ is a constant and lm is the mixing length. In the second method, the velocity scale is defined from the norm of the strain-rate tensor |S|

(17) u 0 = l m | S | = l m U x x 2 + U y y 2 + U z z 2 + 1 2 U x y + U y x 2 + 1 2 U x z + U z x 2 + 1 2 U y z + U z y 2 1 2 .

From the literature (Iungo et al.2017), it appears that in the wake of a wind turbine, the dominating term (in cylindrical coordinates) is Ur. It is supposed herein that these results can be transposed in Cartesian coordinates and are applicable in the MFOR. The velocity scale can thus be written as a function of the derivatives of the axial velocity.

(18) u 0 = l m 1 2 U x y 2 + 1 2 U x z 2

Combining Eqs. (16) and (18) leads to

(19) k x , MF , am = u 0 C μ 1 / 4 2 = l m 2 2 C μ 1 / 2 [ U x , MF y 2 + U x , MF z 2 ] = l m 2 2 C μ 1 / 2 [ - U ( z ) Δ U MF y 2 + - U ( z ) Δ U MF z + 1 - Δ U MF U ( z ) z 2 ] .

In Eq. (19), the last term (1-ΔU)U(z)z represents the produced turbulence due to the interaction between wake-generated shear and atmospheric shear. It is this term that induces a maximum of turbulence at the top tip in cases of high atmospheric shear such as neutral or stable ABLs. Even though an analytical form of this term can be found by assuming U(z) to be a log law or a power law, the convolution product with fc in Eq. (7) did not lead to any analytical solution.

Figure 4Profiles of the modified mixing length (turbulence to shear ratio) for the different simulations.


It was thus decided to neglect shear in the formulation and to add the contribution of the inflow turbulence with a maximum function. This is a strong assumption that impacts the results (see Sect. 5) but allows us to compute the total added turbulence,

(20) ( IV ) am = k x , a , am = max k x , , f c * * k x , MF , am ,


(21) k x , MF , am = U ( z ) l m * 2 Δ U MF , am y 2 + Δ U MF , am z 2 = K MF ( x , z ) y σ y 2 ( x ) 2 + z σ z 2 ( x ) 2 exp - y 2 σ y 2 ( x ) - z 2 σ z 2 ( x ) ,

where KMF=(UClm*)2 and lm* is the modified mixing length lm*=lm/2Cμ1/4. In other words, the modified mixing length is the ratio of the axial turbulence to the quadratic sum of the vertical and horizontal gradients of the axial velocity deficit.

Figure 4 shows the profiles of the modified mixing length in the wake normalized by the modified mixing length upstream of the turbine at hub height: lm*/lm,*(zhub), where lm,* is defined in Eq. (13) and lm* is computed as

(22) l m * = k x , MF / U Δ U x y 2 + Δ U x z 2 .

One can see that there are two distinct values: one inside the wake and one outside the wake. Inside the wake, the value is fairly constant (except in the bottom of the wake where it increases chaotically, probably due to the effect of the ground). It only seems to vary with the streamwise distance, and it was thus assumed that lm*/lm,*(zhub) is only dependent on x. Theoretically, it could be possible to develop a model with two mixing lengths (one for the wake and another for the ambient turbulence), but with such an assumption, no analytical solution of Eq. (7) could be achieved.

Note that in Eq. (21), the error in the near wake due to the Gaussian shape assumption for velocity deficit in the MFOR propagates onto Δkx,MF,am. Using a Gaussian instead of a super-Gaussian function leads to an underestimation of the wake-generated shear and thus to a much weaker but more spread axial turbulence around the blade's tips. Moreover, the model does not account for the atmospheric-generated shear. This phenomenon, which leads to a smaller value of wake-generated turbulence at the bottom tip compared to the top tip, cannot be represented in our model. Finally, the model imposes that kx,MF,am=0 at the centre of the wake, a condition that is not fulfilled in the calibration dataset. Another possible improvement would be to add the streamwise gradient Ux/x in Eq. (18). Despite these flaws, this expression has been chosen since it has an analytical solution of its convolution with the wake centre position distribution fc,am.

3.1.3 Wake meandering

For the PDF of wake meandering, the central limit theorem leads to a Gaussian distribution (Braunbehrens and Segalini2019). The distribution of the wake centre fc is non-axisymmetric, and its variance σf is thus defined in both dimensions:

(23) f c , am ( x , y , z ) = 1 2 π σ f y ( x ) σ f z ( x ) exp - y 2 2 σ f y 2 ( x ) - z 2 2 σ f z 2 ( x ) .

3.2 Velocity in the FFOR

In the following, the dependency of the variables on coordinate x is omitted to shorten the equations.

In Eq. (6), the velocity in the wake is written under its dimensional form, whereas the model chosen in Eq. (15) is written under the velocity deficit form. To relate the velocity to the velocity deficit, it is needed to assume that despite its dependency on z due to the atmospheric shear, the upstream velocity U can be considered to be a constant when applying the 2D convolution product with the wake centre distribution. For any function g(y,z), this simplification can be written as

(24) f c , am ( y , z ) * * U ( z ) g ( y , z ) = U ( z ) f c , am ( y , z ) * * g ( y , z ) .

An analytical form of term (I) can then be deduced from Eqs. (15) and (23):

(25) ( I ) am ( y , z ) = f c , am ( y , z ) * * U ( z ) 1 - Δ U MF , am ( y , z ) = U ( z ) ( 1 - Δ U MF , am y - y c , z - z c f c , am y c , z c d y c d z c ) = U ( z ) 1 - Δ U FF , am .

Since it has been shown in the companion paper that term (II) of Eq. (6) is negligible, we approximate that Ux,FF,am=(I)am. The velocity deficit in the FFOR ΔUFF,am is thus the convolution product of two Gaussian functions. It is known that the convolution product of two normalized Gaussian functions of variance σa2 and σb2 is a normalized Gaussian function of variance σa2+σb2 (Teitelbaum2022). Equation (25) can be written as the product of two convolution products, leading to

(26) Δ U x , FF , am = 2 C π σ y σ z [ 1 2 π σ y exp - y - y c 2 2 σ y 2 1 2 π σ f y exp - y c 2 2 σ f y 2 d y c 1 2 π σ z exp - z - z c 2 2 σ z 2 1 2 π σ f z exp - z c 2 2 σ f z 2 d z c ] = C σ y 2 σ y 2 + σ f y 2 σ z 2 σ z 2 + σ f z 2 exp ( - y 2 2 σ y 2 + 2 σ f y 2 - z 2 2 σ z 2 + 2 σ f z 2 ) .

Even though the reasoning of Braunbehrens and Segalini (2019) is different, it is shown here that their model (Eq. 4) can be found by neglecting term (II) and assuming Eq. (24) as well as Gaussian shapes for the velocity deficit in the MFOR and the wake centre's distribution. This is still a Gaussian form, i.e. Eq. (1), with a FFOR wake width defined as σty,tz=σy,z2+σfy,fz2 and a maximum velocity deficit of

(27) C FF = C σ y 2 σ y 2 + σ f y 2 σ z 2 σ z 2 + σ f z 2 .

To fulfil the conservation of momentum as in Eq. (2), one would need CFF=1-1-CT/(8σtyσtz/D2), which is not the case here. In fact, with this methodology, the conservation of momentum can only be enforced in the MFOR or the FFOR. This is the consequence of neglecting term (II) in the velocity breakdown; however, the induced error is relatively low since term (II) is negligible. Combining Eqs. (25) and (26) leads to our model for the velocity in the wake of a wind turbine:

(28) U x , FF , am ( y , z ) = U ( z ) 1 - C σ y 2 σ y 2 + σ f y 2 σ z 2 σ z 2 + σ f z 2 exp - y 2 2 σ y 2 + 2 σ f y 2 - z 2 2 σ z 2 + 2 σ f z 2 .

3.3 Model for the turbulence in the FFOR

For the turbulence, a model has been found for terms (III) (Eq. 30) and (IV) (Eq. 33). Even though the contribution of the three cross terms of Eq. (7) is not always negligible (see companion paper), the two modelled terms are predominant, and the result of the model limited to these two terms can be compared to the turbulence in the FFOR. The total modelled turbulence is computed here as

(29) k x , am = k x , m , am + k x , a , am .

3.3.1 Meandering term

With the same assumptions as for term (I), it is possible to derive an analytical formulation for term (III) of Eq. (7), i.e. the turbulence induced by wake meandering. The assumption of Eq. (24) must be used again to get U2 out of the convolution product, and Eq. (26) is reused to compute the right-hand side of term (III): UMF^2. On the left-hand side, there is a convolution of the Gaussian function fc,am with ΔUx,MF,am2, which is also a Gaussian function of widths 0.5σy and 0.5σz. It is thus possible to use the fact that the convolution of two Gaussian functions is a Gaussian function (Teitelbaum2022).

(30) ( III ) am = k x , m , am ( y , z ) = f c , am * * U x , MF , am 2 - U x , FF , am 2 = U 2 ( z ) 1 - Δ U x , MF , am y - y c , z - z c 2 f c , am ( y c , z c ) d y c d z c - U 2 ( z ) 1 - Δ U x , FF , am 2 = U 2 ( z ) Δ U x , MF , am 2 y - y c , z - z c f c , am y c , z c d y c d z c - U 2 ( z ) Δ U x , FF , am 2 = C U ( z ) 2 [ σ y 2 σ y 2 + 2 σ f y 2 σ z 2 σ z 2 + 2 σ f z 2 exp - y 2 σ y 2 + 2 σ f y 2 - z 2 σ z 2 + 2 σ f z 2 - σ y 2 σ y 2 + σ f y 2 σ z 2 σ z 2 + σ f z 2 exp - y 2 σ y 2 + σ f y 2 - z 2 σ z 2 + σ f z 2 ]

The shape of term (III) is thus not a double-Gaussian, as it may be interpreted in the literature (Stein and Kaltenbach2019; Ishihara and Qian2018), but rather a Gaussian of width 0.5σ2+σf2 minus a thinner and less pronounced Gaussian of width 0.5σ2+0.5σf2. It can be verified that this expression is always larger than 0; i.e. the meandering only produces turbulence and does not dissipate it.

3.3.2 Rotor-added turbulence term

Combining the chosen models for the wake meandering distribution and the added turbulence in the MFOR (Eqs. 21 and 23) in Eq. (20) leads to an analytical form of the axial rotor-added turbulence:

(31) ( IV ) am ( y , z ) = k x , a , am ( y , z ) = max k x , ; k x , MF , am * * f c , am = max ( k x , ; K MF 2 π σ f y σ f z [ y c σ y 2 2 + z c σ z 2 2 ] exp - y c 2 σ y 2 - z c 2 σ z 2 exp - y - y c 2 2 σ f y 2 - z - z c 2 2 σ f z 2 d y c d z c ) = max ( k x , ; K MF 2 π σ f y σ f z [ y c σ y 2 2 exp - y c 2 σ y 2 - y - y c 2 2 σ f y 2 d y c exp - z c 2 σ z 2 - z - z c 2 2 σ f z 2 d z c + z c σ z 2 2 exp - z c 2 σ z 2 - z - z c 2 2 σ f z 2 d z c exp - y c 2 σ y 2 - y - y c 2 2 σ f y 2 d y c ] ) .

At this point, the added turbulence in the FFOR is the sum of two terms that are identical if the coordinates y and z are swapped. It is the product of two convolutions: the first of f:yy2exp(-y2/σy2) with a Gaussian function and the second of two Gaussian functions. The first convolution product has been solved with a computer algebra tool (Scherfgen2022), and the other has already been solved in Eq. (30). It gives

(32) y c σ y 2 2 exp - y c 2 σ y 2 - y - y c 2 2 σ f y 2 d y c exp - z c 2 σ z 2 - z - z c 2 2 σ f z 2 d z c = 2 π σ f y σ y 2 y 2 + σ f y 4 σ y 2 + 2 σ f y 4 σ y σ y 2 + 2 σ f y 2 5 / 2 exp - y 2 σ y 2 + σ f y 2 2 π σ f z σ z σ z 2 + 2 σ f z 2 exp - z 2 σ z 2 + 2 σ f z 2 = 2 π σ f y σ f z σ y σ z σ y 2 + 2 σ f y 2 σ z 2 + 2 σ f z 2 σ y 2 y 2 + σ f y 2 σ y 2 + 2 σ f y 4 σ y 2 σ y 2 + 2 σ f y 2 2 exp - y 2 σ y 2 + σ f y 2 - z 2 σ z 2 + 2 σ f z 2 .

From Eq. (32), we need to add the same quantity with yz and zy, then factorize and simplify to deduce the model for kx,a,am as follows:

(33) k x , a , am = max [ k x , ; K FF ( y 2 σ y 2 + σ y 2 σ f y 2 + 2 σ f y 4 σ y 2 σ y 2 + 2 σ f y 2 2 + z 2 σ z 2 + σ z 2 σ f z 2 + 2 σ f z 4 σ z 2 σ z 2 + 2 σ f z 2 2 ) exp ( - y 2 σ y 2 + 2 σ f y 2 - z 2 σ z 2 + 2 σ f z 2 ) ] ,


(34) K FF = K MF 1 + 2 σ f y / σ y 2 1 + 2 σ f z / σ z 2 .

It can be noted that in the absence of meandering, i.e. for σfy=σfz=0, the model retrieves its MFOR form (Eq. 21). As for terms (I) and (III), the expression of kx,a,am is based on a Gaussian velocity deficit hypothesis, even in the near wake where the LES wake exhibits a shape closer to a top-hat function. The velocity gradient that is the source of the rotor-added turbulence is thus lower and more spread in the model compared to the actual values. Another issue of the model is that it inadequately takes shear into account due to the assumptions of Eqs. (20) and (24). Indeed, the only source of vertical asymmetry in Eq. (33) is U2, i.e. the velocity shear upstream of the turbine.

4 Model's calibration

The model's equations are based on five variables: the wake widths in the MFORs σy and σz, the modified mixing length lm*, and the standard deviations of the meandering distributions σfy and σfz. Each of these variables needs to be calibrated from the inflow conditions to have a usable model. To do so, the results from the calibration dataset are used. Two versions of the wake meandering calibration are used: the “base” calibration, used when the time series of the upstream velocities are known, and the “engineering” calibration, used when they are not.

4.1 Wake width in the MFOR

As described in Sect. 3, we assume that the wake in the MFOR follows a Gaussian shape function. Moreover, we assume here that the wake is axisymmetric (σy=σz), thus reducing the number of parameters in the model from five to four. The width of the wake in the MFOR is deduced from fitting the function of Eq. (35) on the velocity deficit ΔUMF through a non-linear least-squares method.

(35) f y , z , C 0 , y 0 , z 0 , σ = C 0 + C exp ( - y - y 0 2 2 σ 2 - z - z 0 2 2 σ 2 )

C is fixed as a function of σ (Eq. 2 with σy=σz=σ), and the optimization is run on parameters {C0,y0,z0,σ}, where y0, z0 are the mean wake centres; σ is the wake width (the parameter of interest); and C0 is an offset to help the algorithm.

The resulting wake widths in the MFOR as a function of the downstream distance are plotted with solid lines in Fig. 5 for the six cases of the calibration dataset. Except in the near wake, the wake width evolves linearly with the distance to the turbine. Moreover, the greater the instability (and thus the level of turbulence; see Fig. 2), the greater the slope of this linear relation. Finally, the simulations with degraded thrust seem to have the same slope as the neutral case but with a different origin.

Figure 5Wake width in the MFOR for the different cases of the calibration dataset. Solid lines: results from the LES (Eq. 35); dotted lines: proposed calibration.


For all these reasons, the chosen function for the calibration is the following:

(36) σ / D = ( a I + b ) x D + c β ,

where a, b and c are parameters to fit; I is the total turbulence intensity (Eq. 9); and β=0.51+1-Ct/1-Ct (Bastankhah and Porté-Agel2014). A least-square fit method on the six different curves allowed us to compute the best values of a, b and c (see Table 3). Note that this fit is in the end very similar to what can be found in the literature (e.g. Fuertes et al.2018), except that the slope (parameter a) is smaller because the models of the literature implicitly assume that the meandering is included in the wake expansion.

Table 3Parameters for the wake width in the MFOR.

Download Print Version | Download XLSX

4.2 Modified mixing length

The modified mixing length lm* in Eq. (33) directly drives the amount of turbulence added by the turbine. In Sect. 3, it is shown that this variable in the upper part of the wake is independent of the simulation case when normalized with the upstream modified mixing length. Therefore, the evolution of lm*/lm,* is plotted in Fig. 6, and it shows an approximately linear behaviour with the downstream distance.

Figure 6Normalized modified mixing length for the different cases of the calibration dataset.


In the first approach, it is thus decided to fit the mixing length with a linear function of x/D:

(37) l m * = l m , * d x D + e ,

where lm,* is deduced from Eqs. (13) and (14), in which u* and LMO can be found from a fit of the inflow velocity profile. A least-square fit method on the different curves from Fig. 6 is used to fit Eq. (37). The resulting parameters d and e can be found in Table 4, and the corresponding fitted function is plotted with a dashed black line in Fig. 6. The results are quite satisfying even though all the curves are not perfectly superimposed.

Table 4Parameters for the mixing length.

Download Print Version | Download XLSX

4.3 Wake meandering

The widths of the wake centre's distributions, σfy and σfz, are computed as the standard deviations of the wake centre's coordinates, yc(x,t) and zc(x,t):

(38) σ f y ( x ) = y c ( x , t ) 2 ; σ f z ( x ) = z c ( x , t ) 2 .

The resulting amount of meandering in the horizontal (top figure) and vertical directions (bottom figure) for the six cases of the calibration dataset can be found in Fig. 7. The LES results are plotted with solid lines. Overall, the more unstable the case, the more meandering is found. However, the meandering does not solely depend on the lateral turbulence intensity. In particular, the weakly unstable case has greater vertical meandering than the unstable case, despite having a lower Iz value (see Table 1). It is also worth noting that the reductions in the thrust coefficient have little to no effect on the meandering (all the neutral cases are equivalent).

Figure 7Normalized standard deviation of the wake centre from the LES (solid lines) and results from the base calibration (dashed line) and the engineering calibration (dotted line).


To model the amount of meandering, Braunbehrens and Segalini (2019) propose the following formula:

(39) σ f y , f z ( x ) 2 = 2 k y , z 0 x / U c x U c - ζ A v , w ( ζ ) d ζ ,

where Uc is an advection velocity and 𝒜 is the autocorrelation function of the velocity (the lateral and vertical one respectively). For each case, the results of Eq. (39) are plotted with dashed lines in Fig. 7 for Uc=0.8U. This model works fairly well for the amount of meandering, with the right order of magnitude in each case, and it predicts the different behaviours of the vertical and lateral directions for the unstable and weakly unstable cases. However, such calibration for σfy and σfz is not appropriate for analytical wake modelling because it requires time series of wind velocities at hub height, whereas usually only the mean values are available.

Therefore, we hereby propose (dotted lines in Fig. 7) an engineering-oriented solution to approximate the amount of meandering without access to the unsteady time series of velocities upstream of the turbine. In the first attempts to model the meandering (Ainslie1988), it was proposed that the wake meandering should be a linear function of the inflow wind direction's variance. However, more recent work (Doubrawa et al.2018) has shown that the amount of meandering decreases with the rotor size. Indeed, following the theory of the DWM model, only the eddies larger than the size of the rotor are energetic enough to induce wake meandering. Thus the idea is to calibrate the amount of wake meandering only with eddies larger than this size as follows:

(40) σ f y = k y D U x D

and similarly for σfz. In Eq. (40), kyD is the lateral turbulence with a size that is larger than the diameter of the turbine, i.e. the variance of the wind velocity averaged over a circle of 2 rotor diameters and centred at the hub. Note that the time variance is performed after the spatial averaging.

The issue is that kyD and kzD are not known a priori, and since the stability of the ABL modifies the low-frequency range of the turbulence spectrum, it is expected that the share of the turbulence with a larger size than the rotor to the total turbulence is dependent on the atmospheric stability. This can be observed in Fig. 8 where the ratio between the turbulence larger than a disc of diameter ddisc, ky,zddisc to the total turbulence is computed for ky and kz for each case.

Figure 8Ratio of turbulence averaged over a disc to the total turbulence, for different disc sizes.


Figure 8 highlights two distinct behaviours, depending on the stability conditions: the unstable cases (orange and purple curves) decrease much slower than the near-neutral cases (red, grey, brown and blue), and this phenomenon is particularly marked for the lateral turbulence. It shows that the unstable cases have (in proportion) more low-frequency or large-size eddies than the near-neutral cases.

Even though a fully physical approach would require a measure of the stability and an in-depth study of the turbulence spectrum as a function of the ABL conditions, the objective here is to propose an analytical model that is easy to implement and use. It is thus proposed to model the ratios kyddisc/ky and kzddisc/kz with an analytical function:

(41) k y d / k y = exp - d disc / Γ y , k z d / k z = exp - d disc / Γ z .

A least-square fit is used to determine the value of the parameter Γ. Two different fits are used in order to have one result for unstable cases and one for near-neutral cases. The results are given in Table 5, and Eq. (41) is plotted in Fig. 8 with dotted black and dash-dotted black lines for the neutral and unstable values respectively.

Table 5Parameters for the large-scale turbulence function.

Download Print Version | Download XLSX

Γ can be interpreted as a measure of the large-scale eddies of the atmosphere, even though it is not defined as the integral length scale. The combination of Eqs. (40) and (41) with values from Table 5 is plotted with dotted lines in Fig. 7. Even though the model cannot predict the non-linear behaviour in the far wake, the results remain quite good. Only the weakly unstable case gives poor results, possibly because it is at the edge between the near-neutral and unstable case, and necessitates a Γ value of its own.

Figure 9Results of the analytical velocity model for the different calibrations (dashed blue, dotted red and dash-dotted orange lines) in the neutral case compared to Meso-NH (solid black line) and the I&Q2018 model (dotted red line). Lateral (a–c) and vertical (d–f) profiles are plotted for different positions downstream.


5 Results

In this section, we analyse the results of the new model described in the precedent sections. For the streamwise velocity, the model is described with Eq. (28) and for streamwise turbulence with the sum of Eqs. (30) and (33). This validation is done with the two validation cases (see Table 1), i.e. with the unstable and neutral SWiFT simulations. Three versions of the calibration of σ, σfy, σfz and lm* are shown:

  • The base calibration is defined with Eqs. (36), (38) and (39) as well as values for a, b, c, d and e from Tables 3 and 4. This calibration makes more sense physically but requires the time series of the inflow velocities to determine the autocorrelation necessary to compute σf from Eq. (39). It is plotted with dashed blue lines in Figs. 912.

  • The engineering calibration uses the same equations except for the wake meandering, where Eqs. (40) and (41) are used instead of Eq. (38) and parameter Γ is taken from Table 5. It is plotted with dotted red lines in Figs. 912.

  • Finally, we also propose the “best” version of the model. Knowing that the calibration produces errors, it was interesting to see what would be the results of the “best calibration possible”, i.e. with parameters σ, σfy, σfz and lm* directly taken from the LES of the SWiFT simulation (and not from the calibration deduced from Sect. 4). Obviously, this version of the model cannot be used, but it is helpful to determine if the discrepancies between our model and the LES come from the calibration or the construction of the model itself. It is plotted with dash-dotted orange lines in Figs. 912.

Additionally, in Figs. 912 the reader will find the results directly from Meso-NH (solid black line) and the result from a widely used model of the I&Q2018 model (dash-dotted-dotted purple line), one of the few in the literature that predicts the profiles of both mean streamwise velocity and streamwise turbulence.

5.1 Velocity field

The results for the streamwise velocity field in the FFOR can be found in Figs. 9 and 10 for the neutral and unstable cases respectively. The horizontal (top) and vertical (bottom) profiles of velocity are plotted for the reference LES, results from the literature and the three versions of the aforementioned model's calibration. The three columns are three different positions downstream of the wind turbine: x/D=2, x/D=5 and x/D=8.

Figure 10Results of the analytical velocity model for the different calibrations (dashed blue, dotted red and dash-dotted orange lines) in the unstable case compared to Meso-NH (solid black line) and the I&Q2018 model (dotted red line). Lateral (a–c) and vertical (d–f) profiles are plotted for different positions downstream.


Our model (with any calibration) behaves very similarly to the I&Q2018 model in the neutral case (Fig. 9), and both are accurate compared to the LES data in black. The only discrepancy is in the near wake, where both models assume a Gaussian shape, whereas a super-Gaussian shape (Blondel and Cathelain2020) would be more appropriate. These overall good results confirm that the hypotheses made in Sect. 3 for the velocity in the MFOR and the wake centre distribution are good and that meandering has been correctly computed.

In the unstable case (Fig. 10), the literature model underestimates the wake dissipation, whereas the proposed model is more accurate. This is because the I&Q2018 model only uses CT and Ix as parameters. As shown in Table 1, these values are very similar in the neutral and unstable cases of the validation case, and the I&Q2018 results are thus very similar between the neutral and unstable cases. It cannot predict the increase in meandering under unstable ABL conditions due to higher values of large-scale turbulence in the lateral and vertical directions.

The proposed model is better in that regard, showing a larger wake expansion due to the higher predicted meandering compared to the neutral case. It shows that the determination of the velocity deficit in non-neutral cases necessitates more than only the total streamwise turbulence. In this case, one can note a discrepancy between the best version and the two calibrations of our model. It is due to an overestimation of the wake width in the MFOR for the unstable case (not shown here). Indeed, the neutral and unstable cases have similar wake widths in the MFOR, while they have different total turbulence intensities I (Table 1), and therefore Eq. (36) gives accurate results for the neutral case but overestimates the MFOR wake width in the unstable case. As a result, there is a compensation of error, where the calibration underestimates the velocity deficit, whereas the best version is supposed to slightly overestimate it, resulting in a very good match. Nevertheless, even without this error compensation, the best version still outperforms the literature model.

5.2 Turbulence field

With the same plotting convention as in Figs. 9 and 10, the profiles of turbulence in the horizontal and vertical directions are plotted in Figs. 11 and 12 for the neutral and unstable cases respectively.

Figure 11Results of the analytical streamwise turbulence model for the different calibrations (dashed blue, dotted red and dash-dotted orange lines) in the neutral case compared to Meso-NH (solid black line) and the I&Q2018 model (dotted red line). Lateral (a–c) and vertical (d–f) profiles are plotted for different positions downstream.


Figure 12Results of the analytical streamwise turbulence model for the different calibrations (dashed blue, dotted red and dash-dotted orange lines) in the unstable case compared to Meso-NH (solid black line) and the I&Q2018 model (dotted red line). Lateral (a–c) and vertical (d–f) profiles are plotted for different positions downstream.


In the neutral case (Fig. 11), the I&Q2018 model performs remarkably well. It correctly predicts the location of the double peak in the horizontal direction and the top tip peak in the vertical direction. The proposed model shows less good results: despite the order of magnitude being accurate, the shape of the function is not, and the top tip maximum is not correctly positioned. Since the calibrations do not significantly differ from the best version of the model, this is attributed to modelling errors and not to the calibration. The authors suggest that this deviation originates from the omission of shear in the modelling of rotor-added turbulence, as described in Eq. (33).

The unstable case shows the main shortcomings of the I&Q2018 model and the added value of our model. As shown previously, the I&Q2018 model gives similar results between the unstable and neutral SWiFT cases because they have similar inflow Ix and CT values. However, the Meso-NH simulations show significant differences, in particular the fact that around x=5D, the turbulence profile is unimodal in the unstable case and bimodal in the neutral case. This difference cannot be predicted by the I&Q2018 model as it always assumes a bimodal shape with a maximum at the top tip. However, this change in shape can be predicted by our model since both Eqs. (30) and (33) are bimodal when σσf and unimodal when σσf.

Except for the upstream turbulence profiles, the inflow conditions used in the I&Q2018 are very similar between the neutral and unstable cases. Consequently, the purple profiles are alike in Figs. 11 and 12, whereas the stronger meandering in the unstable case leads to a Gaussian-like turbulence profile, even in the vertical direction. The maximum turbulence is thus no longer located at the top tip but rather at hub height. This property is well predicted by our model but not by the I&Q2018 model, which does not take meandering into account and which predicts quasi-identical behaviours between the neutral and unstable cases. As shown in Figs. 5 and 7, the amount of meandering starts lower but grows faster than the wake width in the MFOR, in particular in unstable conditions. Hence, one can expect a bimodal shape in the near wake and a unimodal shape in the far wake, as seen in Figs. 11 and 12.

However, the calibration of our model leads to an overestimation of the streamwise turbulence, in particular in the near wake. Since there are not many differences between the basic and engineering calibrations, it is not attributed to the meandering calibration (these two calibrations only differ by the meandering modelling) but rather to the overestimated σ in the MFOR, as well as an overestimated lm*. When computed directly from the simulation, the values of lm* are very similar between the neutral and unstable cases, whereas the values of lm,* are much greater in the unstable case (see Fig. 3). Therefore, the value of lm* is overestimated by the model, leading to an overestimation of the rotor-added turbulence and thus to the total turbulence.

The best version of the model gives interesting results, showing that if a better calibration was achieved, in particular for the modified mixing length, the results of the model would be better. This question is further detailed in the next section.

6 Discussion

The previous section shows the results of the model developed in this paper. It is quite good for the streamwise velocity field but can be improved for the turbulence, where the fully empirical model of Ishihara and Qian (2018) shows overall better results in neutral cases but has shortcomings in unstable cases. However, since it is physically based, we know the assumptions of the present model and thus have clear opportunities for improvements. The main ones known by the authors are listed below. Moreover, this work shows that the modification of the velocity and turbulence fields when the ABL stability is modified (and not the Ix or CT) can be predicted. This is a crucial point, as future applications of analytical models, such as digital twins, will require an estimation of the wake velocity and turbulence over small time lapses and not a yearly average like annual energy production (AEP) calculations.

The authors want to emphasize that the presented work is a first step toward a fully physically based model for turbulence profiles that depend on atmospheric stability. In the companion paper, it was shown that the turbulence in the wake of a wind turbine is the sum of several terms, and here we present a methodology to analytically model the most important of these terms. Even though a fully usable calibration is proposed for anyone who would like to test the model, the main purpose of this work is to demonstrate how the rotor-added turbulence and meandering turbulence can be modelled from simple functions.

Figure 13Results of the axial turbulence analytical velocity model in the neutral case, for different values of parameter d in the calibration of lm*.


6.1 Calibration improvement

In Figs. 912, there are discrepancies between the best version of the model and our proposed calibrations. This is particularly true for turbulence, and it is attributed to the calibration of lm*. Contrarily to σ and σf, which can be computed on a wake no matter what, our computation of the modified mixing length lm* makes sense only if it is assumed that the rotor-added turbulence only comes from the wake shear. Additionally, the vertical velocity gradient of the ABL U/z is voluntarily omitted in Eq. (20).

On one hand, all of these assumptions make the measure of lm* a hardly reliable variable. On the other hand, our model is strongly dependent on this parameter. Indeed, the rotor-added turbulence is proportional to the square of lm*. Therefore, a small over- or underestimation of lm* is likely to happen, and it leads to large differences. Figure 13 shows the effect of multiplying parameter d of the mixing length (Table 4) by a factor of 0.8 (dotted red line), 1.2 (dash-dotted orange line) and 1.5 (dash-dotted-dotted purple line) for the basic calibration in the neutral case. It results in large differences from one result to another, showing that even small differences in lm* can drastically change the conclusions.

6.2 Modellization improvements

Besides a better calibration, the model could benefit from conceptual improvements. Indeed, the best version of the model (orange curve in Figs. 11 and 12) does not match the LES results. In other words, even with a “perfect” calibration, the model still misses some features of the turbulence in the wake.

At several points of the reasoning, the atmospheric shear, i.e. the dependence of U on z, is neglected (Eqs. 24 and 20). The first improvement that comes to mind is to model the interaction between atmospheric and wake shear. By doing so, it would be possible to have the reduction in shear near the ground and an increase in shear at the top tip, leading to a smaller value of turbulence at the bottom tip compared to the top tip, as observed in the LES datasets and modelled in the I&Q2018 model. In the current form of the model, the shear is only accounted for through U2 as a factor in km,am and ka,am. This small contribution is compensated for by the upstream turbulence k that is larger at the bottom than at the top, leading to almost symmetric vertical profiles for the model, whereas LES profiles clearly show higher values at the top tip of the wake (at least for the neutral case and in the near wake of the unstable case).

Table 6Calibration parameters of the model.

Download Print Version | Download XLSX

A second improvement that could be done concerns the near wake. As mentioned in Sect. 5, instead of using a simple Gaussian function, a super-Gaussian function would be more accurate. This generic function takes a top-hat form in the near wake and progressively transitions to a Gaussian function as it travels downstream. It was shown in Blondel and Cathelain (2020) that it gives more accurate results in the near wake. Such a function would improve not only the velocity model but also the meandering and rotor-added turbulence terms, which are built upon the velocity model. The latter in particular is a function of the spatial derivative of ΔU using the Gaussian function instead of the super-Gaussian function as done in this work, thus leading to an underestimation of the shear at the edge of the turbine.

For both of these improvements, some solutions were tried: not neglecting the U/z in the derivation of the rotor-added turbulence and using a super-Gaussian function instead of a Gaussian function for the velocity in the MFOR. In both cases, no analytical solution for the models was reached. If such a fully analytical resolution is indeed impossible, an approximated form (for instance, based on LES results) could be proposed in the future.

Finally, modelling the additional terms of Eq. (7), in particular the covariance term (V), could further improve the model. It was shown in the companion paper that this term can represent about 10 % of the total turbulence in the wake and redistributes the turbulence vertically. Given the order of magnitude, this is of lesser importance than the aforementioned points but would also improve the results, or at least the physical accuracy of the model.

7 Conclusions

This work is the second part of a two-step study that aims to model the turbulence in the wake of a wind turbine based on the meandering phenomenon. In the companion paper, the velocity and turbulence in the FFOR were broken down into different terms, some of which were shown to be negligible. In the present work, an analytical model is proposed for the dominating terms of the velocity and turbulence breakdowns, i.e. the meandering turbulence and the rotor-added turbulence. The originality of this work is that it allows for the independent modelling of the effects of meandering (and thus of the ABL stability) and the wake expansion and that it gives the whole turbulence profile rather than only the maximum value. For the velocity, it writes

(42) U x , am ( y , z ) = U ( z ) 1 - C σ y 2 σ y 2 + σ f y 2 σ z 2 σ z 2 + σ f z 2 exp - y 2 2 σ y 2 + 2 σ f y 2 - z 2 2 σ z 2 + 2 σ f z 2

and for the turbulence

(43) k x , am = max [ k x , , C U ( z ) l m * ( x ) 2 1 + 2 σ f y / σ y 2 1 + 2 σ f z / σ z 2 ( y 2 σ y 2 + σ y 2 σ f y 2 + 2 σ f y 4 σ y 2 σ y 2 + 2 σ f y 2 2 + z 2 σ z 2 + σ z 2 σ f z 2 + 2 σ f z 4 σ z 2 σ z 2 + 2 σ f z 2 2 ) exp ( - y 2 σ y 2 + 2 σ f y 2 - z 2 σ z 2 + 2 σ f z 2 ) ] + C U ( z ) 2 [ σ y 2 σ y 2 + 2 σ f y 2 σ z 2 σ z 2 + 2 σ f z 2 exp - y 2 σ y 2 + 2 σ f y 2 - z 2 σ z 2 + 2 σ f z 2 - σ y 2 σ y 2 + σ f y 2 σ z 2 σ z 2 + σ f z 2 exp ( - y 2 σ y 2 + σ f y 2 - z 2 σ z 2 + σ f z 2 ) ] ,

where C=1-1-CT/(8σyσz/D2), CT is the thrust coefficient, D is the turbine diameter, and kx,∞ and U are the variance and mean values of the upstream axial velocity. The model's parameters are the wake widths σy, σz; the amount of meandering σfy,σfz; and the modified mixing length lm*. Two calibrations of these parameters are proposed in Table 6: the first one (base calibration) can be used if time series of the wind velocity are available and the second one (engineering calibration) if they are not. In this table, 𝒜ϕ is the autocorrelation of ϕ, Uc=0.8U and lm,* is found by fitting the inflow velocity profile (Eq. 14). The expressions of velocity and added turbulence in the MFOR used to build Eqs. (42) and (43) can also be used as inputs to the DWM: combined with a synthetic turbulence generation, the unsteady effects of meandering can be modelled.

The model has been tested on two LESs of a single wind turbine wake in a neutral and unstable atmosphere. For the velocity, the results are satisfactory, either in the vertical or in the lateral direction. The model performs better than the model from Ishihara and Qian (2018) in the unstable case as it correctly predicts the increased dissipation due to the increase in meandering. For the turbulence profiles, however, the results are not as good. Since the atmospheric shear was neglected in several steps of the model, the maximum turbulence at the top tip in the neutral case could not be predicted. In the unstable case, the modified mixing length lm* was overestimated, and since the model is very sensitive to this parameter, it resulted in too large values of added turbulence. However, the model of Ishihara and Qian (2018) does not correctly predict the turbulence in the unstable case either. In particular, it still predicts a bimodal shape with a maximum at the top tip in the whole wake, whereas the proposed model successfully transitions from a bimodal to a unimodal shape, according to the LES results.

This is the first step toward a fully analytical, physically based model for turbulence and velocity profiles in the wake of a wind turbine that takes into account atmospheric stability. For future work, the treatment of shear must be improved to model vertical turbulence profiles more realistically. The MFOR velocity deficit function could be replaced by a more accurate function in the near wake to improve the model's results in this region. It would also be interesting to derive an analytical model for the other terms of the turbulence breakdown.

Finally, this model can currently only be used for one turbine, as it predicts only the streamwise velocity and turbulence but necessitates the upstream lateral and vertical turbulence. For the model to be usable for multi-turbines, an expression for every term of the Reynolds stress tensor (or at least the diagonal terms to get the total TKE) would be needed, which implies a model for the lateral and vertical velocities Uy and Uz. This also implies more advanced studies of wake meandering from a turbine working in waked conditions, as most of the wake meandering studies are performed in freestream conditions.

Code and data availability

The Meso-NH code is open source and can be downloaded on the dedicated website. The authors can provide the source code of the modified version 5-4-3 that was used in this work. The data used for the plot presented here and in Part 1 (Jézéquel et al.2024) are available through the following online deposit: (Jézéquel2022). The data for the calibration can be found under the following DOI: (Jézéquel2023). The model equations were written in Python and can be found through the following online deposit: (Jézéquel and Blondel2023).

Author contributions

EJ wrote the analytical model with FB. All the authors worked on the interpretation of the results. The paper was written by EJ with feedback from FB and VM.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


The authors would like to thank the different stakeholders of the MOMENTA project (Jézéquel2023) for allowing them to use the LES data for the calibration of the present model.

Review statement

This paper was edited by Raúl Bayoán Cal and reviewed by three anonymous referees.


Abkar, M. and Porté-Agel, F.: Influence of atmospheric stability on wind-turbine wakes: A large-eddy simulation study, Phys. Fluids, 27, 035104,, 2015. a

Abkar, M., Sørensen, J., and Porté-Agel, F.: An Analytical Model for the Effect of Vertical Wind Veer on Wind Turbine Wakes, Energies, 11, 1838,, 2018. a

Ainslie, J. F.: Calculating the flowfield in the wake of wind turbines, J. Wind Eng. Indust. Aerodynam., 27, 213–224,, 1988. a

Bastankhah, M. and Porté-Agel, F.: A new analytical model for wind-turbine wakes, Renew. Energy, 70, 116–123,, 2014. a, b, c, d, e

Blondel, F. and Cathelain, M.: An alternative form of the super-Gaussian wind turbine wake model, Wind Energ. Sci., 5, 1225–1236,, 2020. a, b, c

Braunbehrens, R. and Segalini, A.: A statistical model for wake meandering behind wind turbines, J. Wind Eng. Indust. Aerodynam., 193, 103954,, 2019. a, b, c, d, e, f, g

Brugger, P., Markfort, C., and Porté-Agel, F.: Field measurements of wake meandering at a utility-scale wind turbine with nacelle-mounted Doppler lidars, Wind Energ. Sci., 7, 185–199,, 2022. a

Cheng, Y., Zhang, M., Zhang, Z., and Xu, J.: A new analytical model for wind turbine wakes based on Monin-Obukhov similarity theory, Appl. Energy, 239, 96–106,, 2019. a

Conti, D., Dimitrov, N., Peña, A., and Herges, T.: Probabilistic estimation of the Dynamic Wake Meandering model parameters using SpinnerLidar-derived wake characteristics, Wind Energ. Sci., 6, 1117–1142,, 2021. a, b

Crespo, A. and Hernandez, J.: Turbulence characteristics in wind-turbine wakes, J. Wind Eng. Indust. Aerodynam., 61, 71–85,, 1996.  a

Doubrawa, P., Martínez-Tossas, L. A., Quon, E., Moriarty, P., and Churchfield, M. J.: Comparison of Mean and Dynamic Wake Characteristics between Research-Scale and Full-Scale Wind Turbines, J. Phys.: Conf. Ser., 1037, 072053,, 2018. a

Du, B., Ge, M., Zeng, C., Cui, G., and Liu, Y.: Influence of atmospheric stability on wind turbine wakes with a certain hub-height turbulence intensity, Phys. Fluids, 33, 055111,, 2021. a, b

Frandsen, S.: Turbulence and turbulence-generated structural loading in wind turbine clusters, PhD thesis, risø-R-1188(EN), DTU, ISBN 87-550-3458-6, 2007. a

Fuertes, F. C., Markfort, C., and Porté-Agel, F.: Wind Turbine Wake Characterization with Nacelle-Mounted Wind Lidars for Analytical Wake Model Validation, Remote Sens., 10, 668,, 2018. a, b

Ishihara, T. and Qian, G.-W.: A new Gaussian-based analytical wake model for wind turbines considering ambient turbulence intensities and thrust coefficient effects, J. Wind Eng. Indust. Aerodynam., 177, 275–292,, 2018. a, b, c, d, e

Iungo, G. V., Santhanagopalan, V., Ciri, U., Viola, F., Zhan, L., Rotea, M. A., and Leonardi, S.: Parabolic RANS solver for low-computational-cost simulations of wind turbine wakes, Wind Energy, 21, 184–197,, 2017. a

Jensen, N.: A note on wind turbine interaction, techreport,Tech. Rep. Risø-M-2411, Risø National Laboratory, Denmark, 16 pp., (last access: 5 December 2022), 1983. a

Jézéquel, E.: Figures data from papers “Breakdown of the velocity and turbulence in the wake of a wind turbine”, parts 1 and 2 [Data set], Zenodo [data set],, 2022. a

Jézéquel, E.: SIMU ATMO, Aeris [data set],, 2023. a, b, c

Jézéquel, E. and Blondel, F.: Implementation of the analytical model deduced from the velocity and turbulence breakdown, version of the reviewed paper, Zenodo [code],, 2023. a

Jézéquel, E., Blondel, F., and Masson, V.: Breakdown of the velocity and turbulence in the wake of a wind turbine – Part 1: Large-eddy-simulation study, Wind Energ. Sci., 9, 97–117,, 2024. a

Keane, A., Aguirre, P. E. O., Ferchland, H., Clive, P., and Gallacher, D.: An analytical model for a full wind turbine wake, J. Phys.: Conf. Ser., 753, 032039,, 2016. a

Keck, R.-E., Maré, M. D., Churchfield, M. J., Lee, S., Larsen, G., and Madsen, H. A.: Two improvements to the dynamic wake meandering model: including the effects of atmospheric shear on wake turbulence and incorporating turbulence build-up in a row of wind turbines, Wind Energy, 18, 111–132,, 2013.  a

Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Aouizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.-M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.-B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.-P., Brossier, C. L., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P., Pergaud, J., Perraud, E., Pianezze, J., Redelsperger, J.-L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V., and Wautelet, P.: Overview of the Meso-NH model version 5.4 and its applications, Geosci. Model Dev., 11, 1929–1969,, 2018. a

Larsen, G. C., Madsen Aagaard, H., Bingöl, F., Mann, J., Ott, S., Sørensen, J. N., Okulov, V., Troldborg, N., Nielsen, N. M., Thomsen, K., Larsen, T. J., and Mikkelsen, R.: Dynamic wake meandering modeling, Risø National Laboratory, ISBN 978-87-550-3602-4, 2007. a

Larsen, G. C., Madsen, H. A., Thomsen, K., and Larsen, T. J.: Wake meandering: a pragmatic approach, Wind Energy, 11, 377–395,, 2008. a, b

Madsen, H. A., Larsen, G. C., Larsen, T. J., Troldborg, N., and Mikkelsen, R.: Calibration and Validation of the Dynamic Wake Meandering Model for Implementation in an Aeroelastic Code, J. Sol. Energ. Eng., 132, 4,, 2010. a, b

Niayifar, A. and Porté-Agel, F.: Analytical Modeling of Wind Farms: A New Approach for Power Prediction, Energies, 625, 012039,, 2016. a

Pope, S. B.: Turbulent Flows, Cambridge University Press,, 2000. a

Scherfgen, D.: Integral Calculator, (last access: 29 April 2022), 2022. a

Stein, V. P. and Kaltenbach, H.-J.: Non-Equilibrium Scaling Applied to the Wake Evolution of a Model Scale Wind Turbine, Energies, 12, 2763,, 2019. a

Teitelbaum, J.: Convolution of Gaussians is Gaussian, (last access: 29 April 2022), 2022. a, b

Xie, S. and Archer, C.: Self-similarity and turbulence characteristics of wind turbine wakes via large-eddy simulation, Wind Energy, 18, 1815–1838,, 2014. a, b

Short summary
Analytical models allow us to quickly compute the decreased power output and lifetime induced by wakes in a wind farm. This is achieved by evaluating the modified velocity and turbulence in the wake. In this work, we present a new model based on the velocity and turbulence breakdowns presented in Part 1. This new model is physically based, allows us to compute the whole turbulence profile (rather than the maximum value) and is built to take atmospheric stability into account.
Final-revised paper