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

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. Jensen, 1983; Larsen et al., 2008; Bastankhah and Porté-Agel, 2014). 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 Archer, 2014):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }U(x,y,z)={\displaystyle \frac{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{\infty }}-\stackrel{\mathrm{‾}}{U}}{{\stackrel{\mathrm{‾}}{U}}_{\mathrm{\infty }}}}=C\left(x\right)\exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y(x{)}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z(x{)}^{\mathrm{2}}}}\right),\text{(2)}& {\displaystyle }& {\displaystyle }C\left(x\right)=\mathrm{1}-\sqrt{\mathrm{1}-{\displaystyle \frac{C_{\mathrm{T}}}{\mathrm{8}{\mathit{\sigma }}_y\left(x\right){\mathit{\sigma }}_z\left(x\right)/D^{\mathrm{2}}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{U}$ is the mean velocity field; ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{\infty }}$ is the mean velocity upstream of the turbine; C(x) is the maximum velocity deficit; C_T 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\left(t\right)=\stackrel{\mathrm{‾}}{X}+X^{\prime }\left(t\right)$.

The stability of the atmospheric boundary layer (ABL) influences the wake recovery (Abkar and Porté-Agel, 2015), 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:

$$\begin{array}{}\text{(3)}& U_{\mathrm{MF}}(x,y,z,t)=U_{\mathrm{FF}}\left(x,y+y_{\mathrm{c}}(x,t),z+z_{\mathrm{c}}(x,t),t\right),\end{array}$$

where subscripts MF and FF denote the velocity fields in the MFOR and FFOR respectively, and y_c(x,t) and z_c(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):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathrm{\Delta }{U}_{\mathrm{FF}}(y,z)=& C{\left[\mathrm{1}+\left({\displaystyle \frac{{\mathit{\sigma }}_{fy}}{{\mathit{\sigma }}_y}}\right)\right]}^{-\mathrm{1}/\mathrm{2}}{\left[\mathrm{1}+\left({\displaystyle \frac{{\mathit{\sigma }}_{fz}}{{\mathit{\sigma }}_z}}\right)\right]}^{-\mathrm{1}/\mathrm{2}}\\ & \cdot \exp \left[-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right],\end{array}\end{array}$$

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 Segalini, 2019; 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é-Agel, 2016) or turbine operating conditions such as C_T and atmospheric shear (Braunbehrens and Segalini, 2019).

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 Hernandez, 1996) or the Frandsen model (Frandsen, 2007) 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 Qian, 2018). 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 $\widehat{a(y,z)}=a(y-y_{\mathrm{c}}(t),z-z_{\mathrm{c}}(t\left)\right)$ for any field a is used. For the present work, it is also important to note that for any field a,

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{‾}}{\widehat{a}}=\stackrel{\mathrm{‾}}{a}**f_{\mathrm{c}},\end{array}$$

where $**$ denotes a 2D convolution and f_c 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).

$$\begin{array}{}\text{(6)}& \stackrel{\mathrm{‾}}{U_{\mathrm{FF}}}=\underset{\left(\mathrm{I}\right)}{\underbrace{\stackrel{\mathrm{‾}}{\widehat{\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}}}}}+\underset{\left(\mathrm{II}\right)}{\underbrace{\stackrel{\mathrm{‾}}{\widehat{U_{\mathrm{MF}}^{\prime }}}}}\end{array}$$

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{k}_{\mathrm{FF}}& =\underset{k_m=\left(\mathrm{III}\right)}{\underbrace{\stackrel{\mathrm{‾}}{{\widehat{\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}}}^{\mathrm{2}}}-{\stackrel{\mathrm{‾}}{\widehat{\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}}}}^{\mathrm{2}}}}+\underset{k_a=\left(\mathrm{IV}\right)}{\underbrace{\stackrel{\mathrm{‾}}{\widehat{k_{\mathrm{MF}}}}}}+\underset{\left(\mathrm{V}\right)}{\underbrace{\mathrm{2}\mathrm{cov}\left(\widehat{\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}},\widehat{U_{\mathrm{MF}}^{\prime }}\right)}}\\ & +\underset{\left(\mathrm{VI}\right)}{\underbrace{\stackrel{\mathrm{‾}}{\widehat{{\left(U_{\mathrm{MF}}^{{}^{\prime }\mathrm{2}}\right)}^{\prime }}}}}\underset{\left(\mathrm{VII}\right)}{\underbrace{-{\stackrel{\mathrm{‾}}{\widehat{U_{\mathrm{MF}}^{\prime }}}}^{\mathrm{2}}}}\end{array}\end{array}$$

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 $\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}$ with f_c. 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 $\widehat{.}$ has no effect) or if there is no turbulence in the MFOR ($U_{\mathrm{MF}}^{\prime }=\mathrm{0}$). Term (III), also written as k_m 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. $U_{\mathrm{MF}}^{\prime }=\mathrm{0}$, Eq. (7) reduces to this term only. Term (IV) is the rotor-added turbulence, which is also written as k_a 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 k_a in the following for consistency with the literature. Term (V) is the covariance of $\widehat{\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}}$ and $\widehat{U_{\mathrm{MF}}^{\prime }}$, 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 Segalini, 2019), 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 (${\stackrel{\mathrm{‾}}{U}}_x$, written as U_x in the following) and the streamwise turbulence ($k_x=\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$) are computed. The proposed model thus only deals with the streamwise velocity and turbulence.

Table 1 List of LES cases.

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équel, 2023).

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⁻² 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 y_c(x,t) and z_c(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 U_ref thus represents a developing ABL without the perturbations of a wind turbine.

Figure 1 Time 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=\mathrm{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

$$\begin{array}{}\text{(8)}& I_{x,y,z}={\displaystyle \frac{\sqrt{k_{x,y,z}}}{U_{\mathrm{\infty },\mathrm{hub}}}},\end{array}$$

and the global turbulence intensity is defined as

$$\begin{array}{}\text{(9)}& I=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left(I_x^{\mathrm{2}}+I_y^{\mathrm{2}}+I_z^{\mathrm{2}}\right)}.\end{array}$$

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 2 Inflow conditions for the calibration cases. (a) Mean velocity profile, (b) mean TKE profile and (c) mean k_x to shear ratio profile. Solid lines: LES results; dotted lines: fit with the Monin–Obukhov law.

Figure 3 Inflow conditions for the validation cases. (a) Mean velocity profile, (b) mean TKE profile and (c) mean k_x 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⁻¹), 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:

$$\begin{array}{}\text{(10)}& U\left(z\right)={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\left(\ln \left(z/z_{\mathrm{0}}\right)+\mathit{\psi }(z,L_{\mathrm{MO}})\right),\end{array}$$

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

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathit{\psi }\left(z,L_{\mathrm{MO}}\right)& =-\mathrm{2}\ln \left(\left(\mathrm{1}+x_u\right)/\mathrm{2}\right)-\ln \left(\left(\mathrm{1}+x_u^{\mathrm{2}}\right)/\mathrm{2}\right)\\ & +\mathrm{2}\arctan \left(x_u\right)-\mathit{\pi }/\mathrm{2},\end{array}\end{array}$$

and $x_u=(\mathrm{1}-\mathrm{15}\cdot z/L_{\mathrm{MO}}{)}^{\mathrm{0.25}}$. Since z₀ 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 L_MO. 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

$$\begin{array}{}\text{(12)}& k={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(k_x+k_y+k_z\right)={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\left(I\cdot U_{\mathrm{\infty },\mathrm{hub}}\right)}^{\mathrm{2}}.\end{array}$$

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 $l_{\mathrm{m},\mathrm{\infty }}^{*}$ 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:

$$\begin{array}{}\text{(13)}& l_{\mathrm{m},\mathrm{\infty }}^{*}={\displaystyle \frac{\sqrt{k_{x,\mathrm{\infty }}}}{\frac{\partial U_{\mathrm{\infty }}}{\partial z}}}.\end{array}$$

However, in the unstable cases, the velocity profile becomes nearly constant above a given height, leading to low values of $\partial U_{\mathrm{\infty }}/\partial z$ and thus very chaotic behaviour of $l_{\mathrm{m},\mathrm{\infty }}^{*}$. To have a more reliable curve, the derivative of U is resolved analytically using Eq. (10):

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial U_{\mathrm{\infty }}}{\partial z}}={\displaystyle \frac{u_{*}}{\mathit{\kappa }z}}{\left(\mathrm{1}-\mathrm{15}z/L_{\mathrm{MO}}\right)}^{-\mathrm{0.25}},\end{array}$$

with L_MO 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 $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{U}}_{x,\mathrm{MF}}$, the turbulence in the MFOR k_x,MF and the meandering distribution f_c. 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 2 Description of the most-used notations in this part and the following.

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é-Agel, 2014; Xie and Archer, 2014), the mean velocity deficit is modelled with the long-established Gaussian velocity deficit (see Eq. 1):

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{U}_{x,\mathrm{MF},\mathrm{am}}(x,y,z)=C\left(x\right)\exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}\left(x\right)}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}\left(x\right)}}\right),\end{array}$$

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 Cathelain, 2020) 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 k_x,a, one needs an analytical form for the turbulence in the MFOR k_x,MF. It was first thought that it would be better to model the added turbulence in the MFOR, i.e. $\mathrm{\Delta }{k}_{x,\mathrm{MF}}=k_{x,\mathrm{MF}}-k_{x,\mathrm{\infty }}$, in order to separate the rotor-added turbulence Δk_x,MF from the ambient turbulence. This procedure was done in the companion paper; however it leads to negative values of Δk_x,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 k_x,MF is not as straightforward as for ΔU_MF 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 Δk_x,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 f_c,am with the chosen shape function for $\mathrm{\Delta }{k}_{x,\mathrm{MF},\mathrm{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 u₀ are combined. In the first, it is assumed to be proportional to the square root of the TKE (Pope, 2000). However in the present work, the 3D TKE is not computed, so it is replaced with the axial turbulence k_x:

$$\begin{array}{}\text{(16)}& u_{\mathrm{0}}=C_{\mathit{\mu }}^{\mathrm{1}/\mathrm{4}}{k}_x^{\mathrm{1}/\mathrm{2}},\end{array}$$

where C_μ is a constant and l_m is the mixing length. In the second method, the velocity scale is defined from the norm of the strain-rate tensor $\mathrm{|}\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{S}}\mathrm{|}$

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{u}_{\mathrm{0}}& =l_{\mathrm{m}}\mathrm{|}\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{S}}\mathrm{|}=l_{\mathrm{m}}\left({\left({\displaystyle \frac{\partial U_x}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial U_y}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial U_z}{\partial z}}\right)}^{\mathrm{2}}\right.\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial U_x}{\partial y}}+{\displaystyle \frac{\partial U_y}{\partial x}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial U_x}{\partial z}}+{\displaystyle \frac{\partial U_z}{\partial x}}\right)}^{\mathrm{2}}\\ & {\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial U_y}{\partial z}}+{\displaystyle \frac{\partial U_z}{\partial y}}\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}\end{array}$$

From the literature (Iungo et al., 2017), it appears that in the wake of a wind turbine, the dominating term (in cylindrical coordinates) is $\frac{\partial U}{\partial r}$. 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.

$$\begin{array}{}\text{(18)}& u_{\mathrm{0}}=l_{\mathrm{m}}\cdot \sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial U_x}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial U_x}{\partial z}}\right)}^{\mathrm{2}}}\end{array}$$

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

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{k}_{x,\mathrm{MF},\mathrm{am}}& ={\left({\displaystyle \frac{u_{\mathrm{0}}}{C_{\mathit{\mu }}^{\mathrm{1}/\mathrm{4}}}}\right)}^{\mathrm{2}}\\ & ={\displaystyle \frac{l_{\mathrm{m}}^{\mathrm{2}}}{\mathrm{2}{C}_{\mathit{\mu }}^{\mathrm{1}/\mathrm{2}}}}\cdot [{\left({\displaystyle \frac{\partial U_{x,\mathrm{MF}}}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial U_{x,\mathrm{MF}}}{\partial z}}\right)}^{\mathrm{2}}]\\ & ={\displaystyle \frac{l_{\mathrm{m}}^{\mathrm{2}}}{\mathrm{2}{C}_{\mathit{\mu }}^{\mathrm{1}/\mathrm{2}}}}\cdot [{\left(-U_{\mathrm{\infty }}\left(z\right){\displaystyle \frac{\partial \mathrm{\Delta }{U}_{\mathrm{MF}}}{\partial y}}\right)}^{\mathrm{2}}\\ & +{\left(-U_{\mathrm{\infty }}\left(z\right){\displaystyle \frac{\partial \mathrm{\Delta }{U}_{\mathrm{MF}}}{\partial z}}+\left(\mathrm{1}-\mathrm{\Delta }{U}_{\mathrm{MF}}\right){\displaystyle \frac{\partial U_{\mathrm{\infty }}\left(z\right)}{\partial z}}\right)}^{\mathrm{2}}].\end{array}\end{array}$$

In Eq. (19), the last term $(\mathrm{1}-\mathrm{\Delta }{U}_{\mathrm{\infty }})\frac{\partial U_{\mathrm{\infty }}\left(z\right)}{\partial 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 f_c in Eq. (7) did not lead to any analytical solution.

Figure 4 Profiles 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,

$$\begin{array}{}\text{(20)}& (\mathrm{IV}{)}_{\mathrm{am}}=k_{x,\mathrm{a},\mathrm{am}}=\max \left(k_{x,\mathrm{\infty }},f_{\mathrm{c}}**k_{x,\mathrm{MF},\mathrm{am}}\right),\end{array}$$

with

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{k}_{x,\mathrm{MF},\mathrm{am}}& ={\left(U_{\mathrm{\infty }}\left(z\right)l_{\mathrm{m}}^{*}\right)}^{\mathrm{2}}\left[{\left({\displaystyle \frac{\partial \mathrm{\Delta }{U}_{\mathrm{MF},\mathrm{am}}}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathrm{\Delta }{U}_{\mathrm{MF},\mathrm{am}}}{\partial z}}\right)}^{\mathrm{2}}\right]\\ & =K_{\mathrm{MF}}(x,z)\left[{\left({\displaystyle \frac{y}{{\mathit{\sigma }}_y^{\mathrm{2}}\left(x\right)}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{z}{{\mathit{\sigma }}_z^{\mathrm{2}}\left(x\right)}}\right)}^{\mathrm{2}}\right]\\ & \cdot \exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}\left(x\right)}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}\left(x\right)}}\right),\end{array}\end{array}$$

where $K_{\mathrm{MF}}=(U_{\mathrm{\infty }}C{l}_{\mathrm{m}}^{*}{)}^{\mathrm{2}}$ and $l_{\mathrm{m}}^{*}$ is the modified mixing length $l_{\mathrm{m}}^{*}=l_{\mathrm{m}}/\sqrt{\mathrm{2}}{C}_{\mathit{\mu }}^{\mathrm{1}/\mathrm{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: $l_{\mathrm{m}}^{*}/l_{\mathrm{m},\mathrm{\infty }}^{*}\left(z_{\mathrm{hub}}\right)$, where $l_{\mathrm{m},\mathrm{\infty }}^{*}$ is defined in Eq. (13) and $l_{\mathrm{m}}^{*}$ is computed as

$$\begin{array}{}\text{(22)}& l_{\mathrm{m}}^{*}={\displaystyle \frac{\sqrt{k_{x,\mathrm{MF}}}/U_{\mathrm{\infty }}}{\sqrt{{\left(\frac{\partial \mathrm{\Delta }{U}_x}{\partial y}\right)}^{\mathrm{2}}+{\left(\frac{\partial \mathrm{\Delta }{U}_x}{\partial z}\right)}^{\mathrm{2}}}}}.\end{array}$$

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 $l_{\mathrm{m}}^{*}/l_{\mathrm{m},\mathrm{\infty }}^{*}\left(z_{\mathrm{hub}}\right)$ 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 $\mathrm{\Delta }{k}_{x,\mathrm{MF},\mathrm{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 $k_{x,\mathrm{MF},\mathrm{am}}=\mathrm{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 $\partial U_x/\partial 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 f_c,am.

3.1.3 Wake meandering

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

$$\begin{array}{}\text{(23)}& f_{\mathrm{c},\mathrm{am}}(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{fy}\left(x\right){\mathit{\sigma }}_{fz}\left(x\right)}}\exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}\left(x\right)}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}\left(x\right)}}\right).\end{array}$$

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

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{f}_{\mathrm{c},\mathrm{am}}(y,z)& **\left(U_{\mathrm{\infty }}\left(z\right)\cdot g(y,z)\right)\\ & =U_{\mathrm{\infty }}\left(z\right)\cdot \left[f_{\mathrm{c},\mathrm{am}}(y,z)**g(y,z)\right].\end{array}\end{array}$$

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

$$\begin{array}{}\text{(25)}& \begin{array}{rl}\left(\mathrm{I}{)}_{\mathrm{am}}\right(y,z)& =f_{\mathrm{c},\mathrm{am}}(y,z)**\left[U_{\mathrm{\infty }}\left(z\right)\left(\mathrm{1}-\mathrm{\Delta }{U}_{\mathrm{MF},\mathrm{am}}(y,z)\right)\right]\\ & =U_{\mathrm{\infty }}\left(z\right)(\mathrm{1}-\int \int \mathrm{\Delta }{U}_{\mathrm{MF},\mathrm{am}}\left(y-y_{\mathrm{c}},z-z_{\mathrm{c}}\right)\\ & \cdot f_{\mathrm{c},\mathrm{am}}\left(y_{\mathrm{c}},z_{\mathrm{c}}\right)\mathrm{d}{y}_{\mathrm{c}}\mathrm{d}{z}_{\mathrm{c}})\\ & =U_{\mathrm{\infty }}\left(z\right)\left(\mathrm{1}-\mathrm{\Delta }{U}_{\mathrm{FF},\mathrm{am}}\right).\end{array}\end{array}$$

Since it has been shown in the companion paper that term (II) of Eq. (6) is negligible, we approximate that $U_{x,\mathrm{FF},\mathrm{am}}=(\mathrm{I}{)}_{\mathrm{am}}$. The velocity deficit in the FFOR ΔU_FF,am is thus the convolution product of two Gaussian functions. It is known that the convolution product of two normalized Gaussian functions of variance ${\mathit{\sigma }}_a^{\mathrm{2}}$ and ${\mathit{\sigma }}_b^{\mathrm{2}}$ is a normalized Gaussian function of variance ${\mathit{\sigma }}_a^{\mathrm{2}}+{\mathit{\sigma }}_b^{\mathrm{2}}$ (Teitelbaum, 2022). Equation (25) can be written as the product of two convolution products, leading to

$$\begin{array}{}\text{(26)}& \begin{array}{rl}\mathrm{\Delta }& U_{x,\mathrm{FF},\mathrm{am}}=\mathrm{2}C\mathit{\pi }{\mathit{\sigma }}_y{\mathit{\sigma }}_z[\int {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_y}}\exp \left(-{\displaystyle \frac{{\left(y-y_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{fy}}}\exp \left(-{\displaystyle \frac{y_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}\right)\mathrm{d}{y}_{\mathrm{c}}\cdot \int {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_z}}\exp \left(-{\displaystyle \frac{{\left(z-z_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{fz}}}\exp \left(-{\displaystyle \frac{z_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\mathrm{d}{z}_{\mathrm{c}}]\\ & =C\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}\exp (-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}).\end{array}\end{array}$$

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 ${\mathit{\sigma }}_{ty,tz}=\sqrt{{\mathit{\sigma }}_{y,z}^{\mathrm{2}}+{\mathit{\sigma }}_{fy,fz}^{\mathrm{2}}}$ and a maximum velocity deficit of

$$\begin{array}{}\text{(27)}& C_{\mathrm{FF}}=C\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}.\end{array}$$

To fulfil the conservation of momentum as in Eq. (2), one would need $C_{\mathrm{FF}}=\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}/(\mathrm{8}{\mathit{\sigma }}_{ty}{\mathit{\sigma }}_{tz}/D^{\mathrm{2}})}$, 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:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{U}_{x,\mathrm{FF},\mathrm{am}}(y,z)& =U_{\mathrm{\infty }}\left(z\right)\left(\mathrm{1}-C\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}\right.\\ & \left.\cdot \exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(29)}& k_{x,\mathrm{am}}=k_{x,\mathrm{m},\mathrm{am}}+k_{x,\mathrm{a},\mathrm{am}}.\end{array}$$

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 $U_{\mathrm{\infty }}^{\mathrm{2}}$ out of the convolution product, and Eq. (26) is reused to compute the right-hand side of term (III): ${\stackrel{\mathrm{‾}}{\widehat{\stackrel{\mathrm{‾}}{U_{\mathrm{MF}}}}}}^{\mathrm{2}}$. On the left-hand side, there is a convolution of the Gaussian function f_c,am with $\mathrm{\Delta }{U}_{x,\mathrm{MF},\mathrm{am}}^{\mathrm{2}}$, which is also a Gaussian function of widths $\sqrt{\mathrm{0.5}}{\mathit{\sigma }}_y$ and $\sqrt{\mathrm{0.5}}{\mathit{\sigma }}_z$. It is thus possible to use the fact that the convolution of two Gaussian functions is a Gaussian function (Teitelbaum, 2022).

$$\begin{array}{}\text{(30)}& \begin{array}{rl}(\mathrm{III}{)}_{\mathrm{am}}& =k_{x,\mathrm{m},\mathrm{am}}(y,z)=\left[f_{\mathrm{c},\mathrm{am}}**U_{x,\mathrm{MF},\mathrm{am}}^{\mathrm{2}}\right]-U_{x,\mathrm{FF},\mathrm{am}}^{\mathrm{2}}\\ & =U_{\mathrm{\infty }}^{\mathrm{2}}\left(z\right)\int \int {\left(\mathrm{1}-\mathrm{\Delta }{U}_{x,\mathrm{MF},\mathrm{am}}\left(y-y_{\mathrm{c}},z-z_{\mathrm{c}}\right)\right)}^{\mathrm{2}}\\ & \cdot f_{\mathrm{c},\mathrm{am}}(y_{\mathrm{c}},z_{\mathrm{c}})\mathrm{d}{y}_{\mathrm{c}}\mathrm{d}{z}_{\mathrm{c}}-U_{\mathrm{\infty }}^{\mathrm{2}}\left(z\right){\left(\mathrm{1}-\mathrm{\Delta }{U}_{x,\mathrm{FF},\mathrm{am}}\right)}^{\mathrm{2}}\\ & =U_{\mathrm{\infty }}^{\mathrm{2}}\left(z\right)\int \int \mathrm{\Delta }{U}_{x,\mathrm{MF},\mathrm{am}}^{\mathrm{2}}\left(y-y_{\mathrm{c}},z-z_{\mathrm{c}}\right)\\ & \cdot f_{\mathrm{c},\mathrm{am}}\left(y_{\mathrm{c}},z_{\mathrm{c}}\right)\mathrm{d}{y}_{\mathrm{c}}\mathrm{d}{z}_{\mathrm{c}}-U_{\mathrm{\infty }}^{\mathrm{2}}\left(z\right)\mathrm{\Delta }{U}_{x,\mathrm{FF},\mathrm{am}}^{\mathrm{2}}\\ & ={\left(C{U}_{\mathrm{\infty }}\left(z\right)\right)}^{\mathrm{2}}[\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}}\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}\\ & \cdot \exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)-{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}\\ & \cdot {\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)]\end{array}\end{array}$$

The shape of term (III) is thus not a double-Gaussian, as it may be interpreted in the literature (Stein and Kaltenbach, 2019; Ishihara and Qian, 2018), but rather a Gaussian of width $\sqrt{\mathrm{0.5}{\mathit{\sigma }}^{\mathrm{2}}+{\mathit{\sigma }}_f^{\mathrm{2}}}$ minus a thinner and less pronounced Gaussian of width $\sqrt{\mathrm{0.5}{\mathit{\sigma }}^{\mathrm{2}}+\mathrm{0.5}{\mathit{\sigma }}_f^{\mathrm{2}}}$. 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:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}(\mathrm{IV}{)}_{\mathrm{am}}& (y,z)=k_{x,\mathrm{a},\mathrm{am}}(y,z)=\max \left(k_{x,\mathrm{\infty }};k_{x,\mathrm{MF},\mathrm{am}}**f_{\mathrm{c},\mathrm{am}}\right)\\ & =\max (k_{x,\mathrm{\infty }};{\displaystyle \frac{K_{\mathrm{MF}}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{fy}{\mathit{\sigma }}_{fz}}}\int \int [{\left({\displaystyle \frac{y_{\mathrm{c}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{z_{\mathrm{c}}}{{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)}^{\mathrm{2}}]\\ & \cdot \exp \left(-{\displaystyle \frac{y_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}-{\displaystyle \frac{z_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)\exp \left(-{\displaystyle \frac{{\left(y-y_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{{\left(z-z_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\mathrm{d}{y}_{\mathrm{c}}\mathrm{d}{z}_{\mathrm{c}})\\ & =\max (k_{x,\mathrm{\infty }};{\displaystyle \frac{K_{\mathrm{MF}}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{fy}{\mathit{\sigma }}_{fz}}}[\int {\left({\displaystyle \frac{y_{\mathrm{c}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)}^{\mathrm{2}}\\ & \cdot \exp \left(-{\displaystyle \frac{y_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}-{\displaystyle \frac{{\left(y-y_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}\right)\mathrm{d}{y}_{\mathrm{c}}\int \exp \left(-{\displaystyle \frac{z_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}}}-{\displaystyle \frac{{\left(z-z_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\mathrm{d}{z}_{\mathrm{c}}\\ & +\int {\left({\displaystyle \frac{z_{\mathrm{c}}}{{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)}^{\mathrm{2}}\exp \left(-{\displaystyle \frac{z_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}}}-{\displaystyle \frac{{\left(z-z_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\mathrm{d}{z}_{\mathrm{c}}\\ & \cdot \int \exp \left(-{\displaystyle \frac{y_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}-{\displaystyle \frac{{\left(y-y_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}\right)\mathrm{d}{y}_{\mathrm{c}}\left]\right).\end{array}\end{array}$$

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:y\to y^{\mathrm{2}}\exp (-y^{\mathrm{2}}/{\mathit{\sigma }}_y^{\mathrm{2}})$ with a Gaussian function and the second of two Gaussian functions. The first convolution product has been solved with a computer algebra tool (Scherfgen, 2022), and the other has already been solved in Eq. (30). It gives

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\int & {\left({\displaystyle \frac{y_{\mathrm{c}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)}^{\mathrm{2}}\exp \left(-{\displaystyle \frac{y_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}}}-{\displaystyle \frac{{\left(y-y_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}\right)\mathrm{d}{y}_{\mathrm{c}}\\ & \cdot \int \exp \left(-{\displaystyle \frac{z_{\mathrm{c}}^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}}}-{\displaystyle \frac{{\left(z-z_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\mathrm{d}{z}_{\mathrm{c}}\\ & ={\displaystyle \frac{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{fy}\left({\mathit{\sigma }}_y^{\mathrm{2}}{y}^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{4}}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{4}}\right)}{{\mathit{\sigma }}_y{\left({\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}\right)}^{\mathrm{5}/\mathrm{2}}}}\\ & \cdot \exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}\right){\displaystyle \frac{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{fz}{\mathit{\sigma }}_z}{\sqrt{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}\exp \left(-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\\ & =\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{fy}{\mathit{\sigma }}_{fz}{\displaystyle \frac{{\mathit{\sigma }}_y{\mathit{\sigma }}_z}{\sqrt{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}\sqrt{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}{\displaystyle \frac{\left({\mathit{\sigma }}_y^{\mathrm{2}}{y}^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{4}}\right)}{{\mathit{\sigma }}_y^{\mathrm{2}}{\left({\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}\right)}^{\mathrm{2}}}}\\ & \cdot \exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right).\end{array}\end{array}$$

From Eq. (32), we need to add the same quantity with y←z and z←y, then factorize and simplify to deduce the model for $k_{x,\mathrm{a},\mathrm{am}}$ as follows:

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{k}_{x,\mathrm{a},\mathrm{am}}& =\max [k_{x,\mathrm{\infty }};\\ & K_{\mathrm{FF}}({\displaystyle \frac{y^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_y^{\mathrm{2}}{\mathit{\sigma }}_{fy}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{4}}}{{\mathit{\sigma }}_y^{\mathrm{2}}{\left({\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}\right)}^{\mathrm{2}}}}+{\displaystyle \frac{z^{\mathrm{2}}{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_z^{\mathrm{2}}{\mathit{\sigma }}_{fz}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{4}}}{{\mathit{\sigma }}_z^{\mathrm{2}}{\left({\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}\right)}^{\mathrm{2}}}})\\ & \cdot \exp (-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}})],\end{array}\end{array}$$

where

$$\begin{array}{}\text{(34)}& K_{\mathrm{FF}}={\displaystyle \frac{K_{\mathrm{MF}}}{\sqrt{\mathrm{1}+\mathrm{2}{\left({\mathit{\sigma }}_{fy}/{\mathit{\sigma }}_y\right)}^{\mathrm{2}}}\sqrt{\mathrm{1}+\mathrm{2}{\left({\mathit{\sigma }}_{fz}/{\mathit{\sigma }}_z\right)}^{\mathrm{2}}}}}.\end{array}$$

It can be noted that in the absence of meandering, i.e. for ${\mathit{\sigma }}_{fy}={\mathit{\sigma }}_{fz}=\mathrm{0}$, the model retrieves its MFOR form (Eq. 21). As for terms (I) and (III), the expression of $k_{x,\mathrm{a},\mathrm{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 $U_{\mathrm{\infty }}^{\mathrm{2}}$, 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 $l_{\mathrm{m}}^{*}$, 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 ΔU_MF through a non-linear least-squares method.

$$\begin{array}{}\text{(35)}& f\left(y,z,C_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}},\mathit{\sigma }\right)=C_{\mathrm{0}}+C\exp (-{\displaystyle \frac{{\left(y-y_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}-{\displaystyle \frac{{\left(z-z_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}})\end{array}$$

C is fixed as a function of σ (Eq. 2 with ${\mathit{\sigma }}_y={\mathit{\sigma }}_z=\mathit{\sigma }$), and the optimization is run on parameters $\mathit{\{}{C}_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}},\mathit{\sigma }\mathit{\}}$, where y₀, z₀ are the mean wake centres; σ is the wake width (the parameter of interest); and C₀ 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 5 Wake 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:

$$\begin{array}{}\text{(36)}& \mathit{\sigma }/D=(aI+b){\displaystyle \frac{x}{D}}+c\sqrt{\mathit{\beta }},\end{array}$$

where a, b and c are parameters to fit; I is the total turbulence intensity (Eq. 9); and $\mathit{\beta }=\mathrm{0.5}\left(\mathrm{1}+\sqrt{\mathrm{1}-C_t}\right)/\sqrt{\mathrm{1}-C_t}$ (Bastankhah and Porté-Agel, 2014). 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 3 Parameters for the wake width in the MFOR.

4.2 Modified mixing length

The modified mixing length $l_{\mathrm{m}}^{*}$ 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 $l_{\mathrm{m}}^{*}/l_{\mathrm{m},\mathrm{\infty }}^{*}$ is plotted in Fig. 6, and it shows an approximately linear behaviour with the downstream distance.

Figure 6 Normalized 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$:

$$\begin{array}{}\text{(37)}& l_{\mathrm{m}}^{*}=l_{\mathrm{m},\mathrm{\infty }}^{*}\left(d{\displaystyle \frac{x}{D}}+e\right),\end{array}$$

where $l_{\mathrm{m},\mathrm{\infty }}^{*}$ is deduced from Eqs. (13) and (14), in which u_* and L_MO 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 4 Parameters for the mixing length.

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, y_c(x,t) and z_c(x,t):

$$\begin{array}{}\text{(38)}& {\mathit{\sigma }}_{fy}\left(x\right)=\sqrt{\stackrel{\mathrm{‾}}{y_{\mathrm{c}}(x,t{{)}^{\prime }}^{\mathrm{2}}}};{\mathit{\sigma }}_{fz}\left(x\right)=\sqrt{\stackrel{\mathrm{‾}}{z_{\mathrm{c}}(x,t{{)}^{\prime }}^{\mathrm{2}}}}.\end{array}$$

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 I_z 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 7 Normalized 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:

$$\begin{array}{}\text{(39)}& {\mathit{\sigma }}_{fy,fz}(x{)}^{\mathrm{2}}=\mathrm{2}{k}_{y,z}\underset{\mathrm{0}}{\overset{x/U_{\mathrm{c}}}{\int }}\left({\displaystyle \frac{x}{U_{\mathrm{c}}}}-\mathit{\zeta }\right){\mathcal{A}}_{v,w}(\mathit{\zeta })d\mathit{\zeta },\end{array}$$

where U_c 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 U_c=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 (Ainslie, 1988), 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:

$$\begin{array}{}\text{(40)}& {\mathit{\sigma }}_{fy}={\displaystyle \frac{\sqrt{k_y^{\mathrm{D}}}}{U_{\mathrm{\infty }}}}{\displaystyle \frac{x}{D}}\end{array}$$

and similarly for σ_fz. In Eq. (40), $k_y^{\mathrm{D}}$ 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 $k_y^{\mathrm{D}}$ and $k_z^{\mathrm{D}}$ 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 d_disc, $k_{y,z}^{{\mathrm{d}}_{\mathrm{disc}}}$ to the total turbulence is computed for k_y and k_z for each case.

Figure 8 Ratio 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 $k_y^{{\mathrm{d}}_{\mathrm{disc}}}/k_y$ and $k_z^{d_{\mathrm{disc}}}/k_z$ with an analytical function:

$$\begin{array}{}\text{(41)}& \begin{array}{rl}& k_y^{\mathrm{d}}/k_y=\exp \left(-d_{\mathrm{disc}}/{\mathrm{\Gamma }}_y\right),\\ & k_z^{\mathrm{d}}/k_z=\exp \left(-d_{\mathrm{disc}}/{\mathrm{\Gamma }}_z\right).\end{array}\end{array}$$

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 5 Parameters for the large-scale turbulence function.

Γ 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 9 Results 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 $l_{\mathrm{m}}^{*}$ 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. 9–12.

• 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. 9–12.

• 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 $l_{\mathrm{m}}^{*}$ 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. 9–12.

Additionally, in Figs. 9–12 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=\mathrm{2}$, $x/D=\mathrm{5}$ and $x/D=\mathrm{8}$.

Figure 10 Results 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 Cathelain, 2020) 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 C_T and I_x 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 11 Results 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 12 Results 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 I_x and C_T values. However, the Meso-NH simulations show significant differences, in particular the fact that around x=5 D, 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 $l_{\mathrm{m}}^{*}$. When computed directly from the simulation, the values of $l_{\mathrm{m}}^{*}$ are very similar between the neutral and unstable cases, whereas the values of $l_{\mathrm{m},\mathrm{\infty }}^{*}$ are much greater in the unstable case (see Fig. 3). Therefore, the value of $l_{\mathrm{m}}^{*}$ 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 I_x or C_T) 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 13 Results of the axial turbulence analytical velocity model in the neutral case, for different values of parameter d in the calibration of $l_{\mathrm{m}}^{*}$.

6.1 Calibration improvement

In Figs. 9–12, 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 $l_{\mathrm{m}}^{*}$. Contrarily to σ and σ_f, which can be computed on a wake no matter what, our computation of the modified mixing length $l_{\mathrm{m}}^{*}$ 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 $\partial U_{\mathrm{\infty }}/\partial z$ is voluntarily omitted in Eq. (20).

On one hand, all of these assumptions make the measure of $l_{\mathrm{m}}^{*}$ 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 $l_{\mathrm{m}}^{*}$. Therefore, a small over- or underestimation of $l_{\mathrm{m}}^{*}$ 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 $l_{\mathrm{m}}^{*}$ 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 $U_{\mathrm{\infty }}^{\mathrm{2}}$ as a factor in k_m,am and k_a,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 6 Calibration parameters of the model.

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 $\partial U_{\mathrm{\infty }}/\partial 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

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{U}_{x,\mathrm{am}}(y,z)& =U_{\mathrm{\infty }}\left(z\right)\left(\mathrm{1}-C\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}\right.\\ & \left.\cdot \exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\right)\end{array}\end{array}$$

and for the turbulence

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{k}_{x,\mathrm{am}}& =\max [k_{x,\mathrm{\infty }},{\displaystyle \frac{{\left(C{U}_{\mathrm{\infty }}\left(z\right)l_{\mathrm{m}}^{*}\left(x\right)\right)}^{\mathrm{2}}}{\sqrt{\mathrm{1}+\mathrm{2}{\left({\mathit{\sigma }}_{fy}/{\mathit{\sigma }}_y\right)}^{\mathrm{2}}}\sqrt{\mathrm{1}+\mathrm{2}{\left({\mathit{\sigma }}_{fz}/{\mathit{\sigma }}_z\right)}^{\mathrm{2}}}}}\\ & \cdot ({\displaystyle \frac{y^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_y^{\mathrm{2}}{\mathit{\sigma }}_{fy}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{4}}}{{\mathit{\sigma }}_y^{\mathrm{2}}{\left({\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}\right)}^{\mathrm{2}}}}+{\displaystyle \frac{z^{\mathrm{2}}{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_z^{\mathrm{2}}{\mathit{\sigma }}_{fz}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{4}}}{{\mathit{\sigma }}_z^{\mathrm{2}}{\left({\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}\right)}^{\mathrm{2}}}})\\ & \cdot \exp (-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}})]+{\left(C{U}_{\mathrm{\infty }}\left(z\right)\right)}^{\mathrm{2}}\\ & \cdot [\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}}\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}}\exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\right)\\ & -{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\sigma }}_z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}}\exp (-{\displaystyle \frac{y^{\mathrm{2}}}{{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{fy}^{\mathrm{2}}}}-{\displaystyle \frac{z^{\mathrm{2}}}{{\mathit{\sigma }}_z^{\mathrm{2}}+{\mathit{\sigma }}_{fz}^{\mathrm{2}}}})],\end{array}\end{array}$$

where $C=\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}/(\mathrm{8}{\mathit{\sigma }}_y{\mathit{\sigma }}_z/D^{\mathrm{2}})}$, C_T is the thrust coefficient, D is the turbine diameter, and k_x,∞ 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 $l_{\mathrm{m}}^{*}$. 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 ϕ, U_c=0.8U_∞ and $l_{\mathrm{m},\mathrm{\infty }}^{*}$ 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 $l_{\mathrm{m}}^{*}$ 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 U_y and U_z. 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.