The neural network is supplied with data extracted from the high-fidelity LES which allows it to predict the rotor effective transverse velocity, w _WT(t). The corresponding ANN input consists of the instantaneous rotation speed of the wind turbine, Ω(t) as well as its projected yaw, γ(t), and global pitch angles, β(t). Finally, following Bertelè and Bottasso (2020), the reconstructed blade loads harmonics are also appended to this input state vector: m ̄ t m c t m s t = 1 3 1 3 1 3 2 3 cos ⁡ θ t 2 3 cos θ t − 2 π 3 2 3 cos θ t + 2 π 3 2 3 sin ⁡ θ t 2 3 sin θ t − 2 π 3 2 3 sin θ t + 2 π 3 × m 1 t m 2 t m 3 t (4) with θ(t), the instantaneous azimuthal position of blade 1. In this expression, a Multi-Blade Coordinate (MBC) transform derived by Coleman and Feingold (1958) is applied in order to map the edgewise, m _e,b and flapwise bending moments, m _f,b of the blades, b = 1, 2, 3 to their 1-P harmonics.

The resulting state vector fed to the neural network thus comprises 11 different fields: i t = Ω t , γ c / s t , β c / s t , m e / f , 0 t , m e / f , c t , m e / f , s t (5) where t is the temporal index and ϕ _c/s denotes the projected cos(ϕ) and sin(ϕ) components, respectively. Data is sampled at a frequency f _ann and the input and output state vectors are eventually concatenated into a single data base entry: S t = i t f ; w t f (6) where (⋅) ^f denotes that an Exponential Moving Average filter with time constant τ _ann has been applied (here both to the input and output data).

We finally introduce some information about the past history by concatenating several input filtered states into a single data base entry: S t h = i t − T a n n f , … , i t − 1 f , i t f ; w t f . (7) T _ann determines the number of successive time indices considered: the larger T _ann , the longer the neural network memory. We should however note that, in this configuration, each data base entry is considered as an independent data sample, this may very rapidly increase the complexity of the input data.

The different date base entries are then collected inside one large data base whose different fields, ϕ, are normalized using a minmax scaler: ϕ n o r m = ϕ − min ϕ max ϕ − min ϕ (8) which casts all the database values into the interval [0, 1] thereby facilitating the neural network convergence.

This data formatting procedure is then applied in order to aggregate a LES data base. The latter is splitted into three subsets: a training, a validation and a testing set. The training set provides never-seen data against which the performances of the trained neural network can be validated. The testing set on the other hand constitutes the bulk of the data base and is used for training purposes. Finally, the validation set provides a metric to avoid overfitting: the optimization procedure should be stopped when the loss metric evaluated on the validation set has reached its minimum. This procedure ensures that the calibrated neural network generalizes well to other data sets.

Different LES are performed in order to assemble the LES database. Their results are uniformly sampled at a f _ann = 0.25 Hz and then appended to the data base. The description of the simulation setup used for each data base subset (i.e., validation, training and testing sets) is presented hereunder.

• The wind farm layout from which the training and validation data set are extracted consists of 19 turbines with a random spacing ranging from 4D to 10D thereby accounting for both waked and unwaked conditions. 6,250 s long simulations are carried for each of the ABL precursors resulting in a total 89,000 data base entries. 20% of this data is used for validation purposes while the remaining 80% are fed to the ANN as part of its learning. Within both subsets, the data is shuffled so that the neural network does not overfit the last wind turbine of the data set. Indeed, if no shuffling was applied, the ANN would converge toward an optimum representative of the last wind turbine encountered as part of its training which may not be corresponding to the global optimum.

A Multi-Layer Perceptron (MLP) is used to learn the mapping between the input state vector consisting of the 11 original inputs, i t f , and the target effective transverse velocity output, w WT , t f . As illustrated on Figure 2, MLPs are a class of ANNs consisting of at least three fully-connected layers of nodes: an input, a hidden and an output layers. Both the hidden and outputs layers use non-linear activation functions thereby allowing this class of neural networks to capture non-linear behaviors.

The different layers are connected through weight matrices whose values are updated as part of the training procedure until convergence is reached. Learning is performed through the back-propagation of the reference loss function through the network during the training step. The complete training data set is first split into mini-batches each consisting of N _ann,b randomly picked data base samples. At each iteration, a new mini-batch is fed to the MLP and the network weights are updated using the standard Adam optimizer algorithm with learning rate lr _ann in order to minimize the target loss function, L _MSE: L MSE = ∑ t w ̂ WT , t − w W T , t f 2 (9) which corresponds to the Mean Square Error (MSE) with w ̂ WT , t and w WT , t f the estimated and reference rotor effective transverse velocities, respectively.

The MLP used here comprises N _ann,L hidden layers of N _ann,N neurons. The latter use Rectified Linear Unit (ReLU) activation functions preceded by a Normalization Layer in order to reduce the training time (Ba et al., 2016). Finally, the Adam optimization algorithm is used to update the layers weight matrices during the training procedure. The resulting neural network architecture is implemented and trained using PyTorch, an open source machine learning framework for deep neural networks (Paszke et al., 2019). A heuristic approach was adopted in order to compute the hyperparameters. Different configurations were tested and only the set of hyperparameters achieving the lowest MSE on the validation set was conserved. We should therefore note that a more exhaustive hyperparameters search may lead to improved performances of the ANN. However, it would inevitably inflate the computational cost. The neural network training and subsequent hyperparameters tuning are not discussed in more details here. The optimized network hyperparameters are: N _ann,L = 3, N _ann,N = 64, T _ann = 1, τ _ann = 8s, lr _ann = 3.5 × 10⁻⁵ and N _ann,b = 64.

The freestream flow field is represented by a series of freestream particles, F _i , shed at successive timesteps t _i , and described by a position, x F i ( t ) and a characteristic velocity, u F i ( t ) . The source state of a given particle is provided by the measurements of the wind turbine it emanates from: x F i t i = x WT (12a) u F i t i = α w u F i − 1 t i + 1 − α w u ̂ WT t i (12b) where α _w is the Exponential Moving Average parameter computed from its filtering timescale, τ _w .

The discrete information conveyed by the freestream particles is interpolated using a radial basis function based spatio-temporal interpolation scheme similar to Rott et al. (2020): u f x , t = ∑ i u F i p F i x , t ∑ i p F i x , t (13) where the interpolation weights, p F i , are computed as follows: p F i x , t = exp − r x , t 2 2 σ r 2 + ξ x , t 2 2 σ ξ 2 exp − t − t i 2 2 σ t 2 . (14) ξ(x, t) and r(x, t) are the local streamwise and radial curvilinear coordinates of the point considered in the particle frame of reference. The weight relative to a particle thus fades out as the particle gets further away from the point of evaluation for the field u _f (x, t). Furthermore, the information contained by a particle is considered less relevant as the particle progressively gets older, i.e., t ≫ t _i . This temporal decay is tuned by σ _t , the scaling constant of the temporal distance. The σ _r and σ _ξ parameters, on the other hand, govern how localized the information should be in the streamwise and radial directions: a low value means the particles have a small trust zone. Far away from the particles sources, the weights collapse and the mean average flow features are recovered.

We note that this framework allows different values of the spatio-temporal scheme weighting parameters to be selected depending on the nature of the phenomenon studied. Indeed, frozen turbulence demands the freestream velocity field to be advected as smoothly as possible whereas the coherence with the self-induced velocity of the wake requires a narrower interpolation kernel. A distinction is thus made for the interpolation of the advection velocities for the large-scale freestream field, u _f , and the local wake, u ′ f ; those are evaluated using the weighting parameters [σ _r , σ _ξ , σ _t ] and [σ _r , σ ξ ′ , σ t ], respectively, with σ ξ ′ < σ ξ .

Following Eq. 11, particle positions are transported by an effective freestream convection velocity u ̃ f : x F i t + Δ t = x F i t + u ̃ f x F i , t Δ t (15) while their characteristic velocity, u F i , is assumed to remain unchanged. The effective freestream convection velocity is obtained by subtracting the speed deficit from the freestream velocity field itself: u ̃ f x F i , t = u f x F i , t − C f π R W 2 ∫ S Δ u x , t d s (16) where C _f is a tuning constant and S denotes the circular area of radius R W = 2 D 2 centered at x F i and normal to the time-averaged flow streaklines defined by Eq. 15.

A wake particle, W _i , shed at time t _i , is described by a position, x W i ( t ) , an orientation, n W i ( t ) , a curvilinear coordinate along the wake centerline, ξ W i ( t ) and finally a thrust coefficient, C T , W i ( t ) , and a rotor-effective turbulence intensity, T I W i (t).

The position of the particle, x W i = x W i e ̂ x , WT + z W i e ̂ z , WT , is measured in the wind farm inertial frame and initially coincides with that of the wind turbine hub. Its remaining states are directly retrieved from the wind turbine state at shedding time. The subsequent source state of the wake particle is therefore similar to equation and reads: ξ W i t i = 0 (17a) x W i t i = x WT (17b) n W i t i = cos γ t i , − sin γ t i (17c) T I W i t i = α w T I W i − 1 t i + 1 − α w T ̂ I WT t i (17d) C T , W i t i = α w C T , W i − 1 t i + 1 − α w C ̂ T , WT t i (17e) where the orientation of the wake particle is computed based on the turbine yaw angle at shedding time, γ(t _i ).

Each wake particle, W _i , is associated to a relative speed deficit field: Δ u W i ( t , r W i ( x ) ) where r W i ( x ) and t denote the radial position of x relative to the wake particle and the time considered, respectively. The orientation of the particle, n W i , is introduced in order to allow the computation of the radial position of an arbitrary point, x, in the wake particle frame: r W i x = x − x W i − x − x W i ⋅ n W i n W i . (18)

The speed deficit field relative to the wake particle is assumed quasi-steady and then paramatrized based on the analytical expression proposed by Bastankhah and Porté-Agel (2014) unlike most DWM implementations (Larsen et al., 2007; Jonkman and Shaler, 2021) that rather rely on the tuning-intensive Ainslie viscosities wake model. Indeed, even though both models were shown to predict accurately the far-wake speed deficit (Zhan et al., 2020), the former overall demonstrates improved computational time and facilitated tuning due to its limited set of parameters.

The Bastankhah (Bastankhah and Porté-Agel, 2014) wake model is derived through the momentum conservation equation assuming a negligible pressure term. Consistent with the thin shear layer assumption, the speed deficit is then parametrized as a Gaussian whose characteristics are expressed as functions of the thrust coefficient and turbulent intensities: Δ u W i x , t = u ̂ WT t i 1 − 1 − C T , W i 8 σ W i t / D 2 × exp − 1 2 σ W i t / D 2 r W i x D 2 . (19) The width of the Gaussian, σ W i ( t ) , grows linearly with the material coordinate ξ W i ( t ) through a growth rate k that is itself assumed to be a linear function of the turbulent intensity, T I W i (Doekemeijer et al., 2019; Duc et al., 2019): σ W i t = k T I W i ξ W i t + ε W i D = a k + b k T I W i ξ W i t + ε W i D (20) with a _k and b _k tuning constants. The initial value of the width is σ _Wi,0 = ɛ _Wi D and determines the characteristic diameter of the streamtube at the outlet of the near-wake region: ε W i = ε 1 2 1 + 1 − C T , W i 1 − C T , W i (21) with ɛ = 0.2 (Bastankhah and Porté-Agel, 2014). The value of these constants has been extensively studied in the literature (e.g., Niayifar and Porté-Agel, 2015). Most of these studies have however targeted the development of time-averaged models where k captures the time-averaged effects of both wake expansion and wake meandering. The value used here should therefore be lower than the traditional values of the expansion constant as wake meandering is already accounted for through our Lagrangian framework.

The speed deficit estimated by Eq. 19 is however not valid in the near wake region. Other analytical field models introducing various near-field wake correction strategies such as the ones proposed by Blondel et al. (2020) or Keane (2021) were thus investigated. They were eventually dismissed as, in the near-wake, only the effective convective velocity scale is important. The exact shape of the speed deficit is not pertinent as no wind turbine is generally located within this wake region. Following Bastankhah and Porté-Agel (2016), the characteristic speed deficit within the near wake is thus recovered using the potential core theory. This theory is consistent with the present framework because it provides coherent estimates of the wake velocity while requiring minimal tuning even though it largely overlooks the near-wake physics. The wake potential core is described as a region of uniform wake speed deficit, Δ u W i N W t , r W i ( x ) , whose inlet is located at the wind turbine location. Due to its interactions with the ambient flow, it gradually gets smaller until ξ W i , 0 = ξ W i ( t ) where the Gaussian self similarity solution is recovered: Δ u W i N W x , t = u ̂ WT t i 1 − 1 − C T , W i for 2 r W i x D < 1 − ξ W i t ξ W i , 0 . (22)

For isolated wakes, the global speed deficit field, Δu(x, t), is obtained by projecting the local speed deficit, Δ u W j , along the wake particle orientation, n W j (Jonkman and Shaler, 2021): Δ u x , t = Δ u W j t , r W j x n W j . (23) W _j is the particle whose plane is the closest to x. In the context of multiple superposing wakes, an altered version of the standard root-sum-square weighting superposition strategy (Gebraad and van Wingerden, 2014; Duc et al., 2019; Jonkman and Shaler, 2021) is applied: Δ u x , t = ∑ k ∈ x , z R ∑ j ∈ S Δ u W j 2 t , r W j x × n W j , k | n W j , k | e ̂ k , W T . (24)

R ( ϕ ) = sign ( ϕ ) | ϕ | is the sign preserving square root operator and S denotes the set of the closest wake particles for each superposing wake. This procedure allows for an efficient merging of both streamwise and transverse components of the speed deficit.

The interpolated wake field can then be used to update the wake particles states using a Lagrangian advection scheme as in Eq. 15: x W i t + Δ t = x W i t + u ̃ w x W i , t Δ t . (25) The wake curvilinear coordinate, ξ W i , is updated accordingly: ξ W i t + Δ t = ξ W i t + u ̃ w x W i , t Δ t . (26) Consistent with Hill’s vortex theory and experimental investigations (Machefaux et al., 2014; Larsen et al., 2020), the wake advection velocity is assumed to be directly related to the maximum speed deficit observed. Following Eq. 24, the effective wake convection velocity, u ̃ W i ( t ) , is thus readily computed from the superposition of the local reconstructed freestream field and local speed deficit: u ̃ w t = u ′ f x W i , t − C w Δ u x W i , t (27) where C _w is a tuning constant.

The remaining particle states used to compute the local speed deficit, T I W i and C T , W i , are conserved as the particle moves downstream.

The identified constant for the wake convective velocity, C _w , is in line with the values reported in the literature. Larsen et al. (2020) found C _w = 0.4 by exploiting the wake/ring-vortex analogy along Hill’s vortex theory to approximate the wake self-induction. This theory is investigated in more details by Machefaux et al. (2014) who report slightly higher values (0.63 and 0.51) based on LES and experimental data. Keck et al. (2014), on the other hand, did not directly account for wake-induced velocity and adopted the wake transport velocity at 80% of free-stream velocity.

An initial guess for the wake expansion constants is retrieved from the high-fidelity data. The value of the precursor velocity field is first subtracted from the wind farm one thereby isolating the speed deficit from the ambient flow. The wake centerline is then extracted from the wind farm LES through the algorithm introduced by Coudou (2021). This allows to translate the speed deficit, initialy in the inertial LES frame, to the meandering frame of reference. A Gaussian speed deficit is finally fitted on the profiles obtained for various time windows and ambient conditions. Even so, the wake expansion constant obtained (ɛ = 0.2, a _k = 0.022 and b _k = 0.14) still overestimated wake recovery due to the inaccuracies introduced by the wake tracking algorithm and by the flow asymmetry. As expected, the calibrated constants are lower than the ones presented in the literature (Niayifar and Porté-Agel, 2015).

Following the frozen turbulence hypothesis, the freestream convective constant, C _f , is traditionally set to zero: this is equivalent to postulating a one-way coupling between the freestream velocity field and the wake. While a uniform convective velocity was indeed observed for freestream flow, the analysis of the reference LES rotor-averaged quantities indicates a reduction of the effective freestream convection speed for waked conditions (Figure 4). The cross-correlation between the histories of the rotor-scale flow quantities and the rotor measurements is evaluated at several downstream locations, ξ _f . The time offset, T _f , which achieves the maximum cross-correlation, and ξ _f are then combined to define the velocity scale, U f = ξ f T f . This analysis is repeated for different wind farm configurations and demonstrates a drop in U _f in the near wake while the freestream convection speed is progressively recovered as one considers positions further downstream.

A quantitative analysis of wake centerline characteristics is now performed for WT0 and WT3. Figure 9 shows the histories of their transverse location, z _C , evaluated at several downstream positions (6D, 9D and 12D). As differences in tracked wake centerline positions will likely occur depending on the size of the mask considered for the LES post-processing, different mask sizes are evaluated and used to define the so-called wake centerline envelope. Highly asymmetric or incoherent wakes are logically harder to track and hence tend to lead to larger differences in wake centerline position depending on the size of the mask selected. A wider wake envelope is therefore indicative of an asymmetric or incoherent wake. In both the high- and the low-TI case, less agreement is obtained between the different wake tracking masks for WT3 than for WT0 resulting in broader wake envelopes on average. The reference mask is then chosen quite wide, i.e., σ _x = σ _y = σ _z = D/2, in order to prevent the spurious numerical behaviors of the wake centerline detection.

The wake centerline computed by the model is in good overall agreement with the wake centroid extracted from the LES. The model appears to correctly advect the large wake features and to reproduce consistent wake deflection amplitudes. The accuracy of the wake centerline prediction is however limited by that of the underlying flow estimator. The transverse velocity estimations for WT0 clearly outperform those obtained for WT3 (Figure 5) which directly translates into a better tracking of the wake center downstream of WT0. This increased mismatch also reflects the higher complexity of the wake downstream of WT3 as indicated by the broader wake envelope. Moreover, the performances of the Lagrangian flow model degrade as the wavelength of the meandering mode considered decreases. As noted previously, this may be partially explained by the reduced turnover time of small-scale eddies. Small eddies quickly evolve as they travel downstream thereby making the information conveyed by the particles less relevant as they age. As a result, the signal cross-correlations, ρ z C , tend to decrease with the downstream distance.

The wake centerlines are further compared in terms of their statistics on the basis of the work of Foti et al. (2018). Each wake centerline is first low-pass filtered spatially with a cut-off length of 1D. The maxima and minima of the subsequent signals are then tracked and the statistics are evaluated. The amplitude, A, is defined as half the transverse distance between two consecutive maxima-minima while the wavelength, λ, is computed as the streamwise distance between two successive maxima or minima. Finally the Strouhal number is introduced as: S t = u ̃ W U A B L ⋅ D λ (28) where u ̃ W is the wake convection speed. These quantities are then binned, with respect to the downwind distance from the rotor. The resulting wake statistics are finally averaged over the turbines WT0 to WT4 and the results are plotted in Figure 10. The greyed areas (i.e., x D < 3 and 11 < x D ) denote areas that are not exploitable owing to the method used for the computation of the wake statistics. It is nonetheless plotted for the sake of transparency.

The method is first applied to the wake centerline time series extracted from the LES using a Gaussian mask with σ x = σ y = D 2 and σ _z = D. The LES results demonstrate good agreement with previous experimental investigations by Coudou (2021). However, while they report strong influence of the ABL characteristics on the wake meandering statistics, we do no distinguish such a trend in the current context: the wake statistics exhibit a similar behavior for the low- and high-TI configurations.

When comparing the Lagrangian flow model statistics to that of the LES, reasonable match is obtained regarding the predicted binned amplitudes. Both approaches result in a drop in amplitude beyond 8D and the binned amplitude probability distributions look similar.

The wake convection speed is slightly overestimated by the Lagrangian flow model even though both curves present similar trends. The wake recovery leads to higher u ̃ w in the far wake. We however note that wake convection speed behavior in the near wake ( x D < 4 ) is not consistent with that of LES. A more consistent behavior is obtained when the near wake correction is disabled. This however leads to a significant time offset between the LES and the flow model time responses.

We finally investigate the mean wavelength: this metric is roughly constant across all downstream distances but appears inflated by the model while the LES produces values in line with the ones obtained by Coudou et al. (2018) (i.e., 2.5 < λ D < 3.25 for isolated wind turbines). This in turn results into an underestimation of the Strouhal number. Wake lengths under 3D are mostly overlooked by the model. Improving the resolution (i.e., the particle shedding frequency) can substantially improve this match but does not translate into a better correspondence of the time series. If a wider Gaussian mask σ _x = σ _y = σ _z = D/2 is employed for the LES wake centerline extraction, better agreement is recovered regarding the binned wavelengths while similar agreement is conserved for the other statistics.

The rotor-effective streamwise velocity, u _RE , describes how the rotor-averaged velocity of a fictive wind turbine placed at some location downstream would evolve with time. The 2D flow field inferred by the Lagrangian Flow model is extrapolated outside the reference x − z plane assuming wake axi-symmetry and neglecting the flow shear while the reference u _RE value is simply extracted from the LES using a circular mask. Both curves are plotted on Figure 11.

The rotor effective velocity can be considered as the superposition of the freestream velocity and of the speed deficit. The former dictates the slow velocity fluctuations caused by the heterogeneous inflow and the latter, strongly influenced by the meandering phenomenon, leads to the marked drops in velocity observed within the wake. Specifically, wake meandering deviates the wake into or away from the fictious downstream rotors, thus successively decreasing and increasing the rotor effective wind speed. This phenomena is remarkably well captured downstream of WT0 for the low-TI configuration (and to some more limited extent for the high-TI one). In the low-TI case, the extreme wake centerline deflection occurring 9D behind WT0 at times 45T _C , 65T _C and 105T _C coincides with a sharp increases in u _RE . Similar peaks can be found in the high-TI case for times 40T _C , 45T _C and 105T _C . Still, most small-scale velocity fluctuations are overlooked by the rotor as they corresponds the eddies whose characteristic dimensions fall beyond the model resolution. The model performances regarding this metric are however not consistent across all wind turbines: turbines whose meandering is poorly captured are obviously less likely to produce accurate predictions u _RE . Nevertheless, in all cases, the initial drop in streamwise velocity as well as the slow dynamics of the deficit recovery are well captured by the Lagrangian flow model.