Modeling dynamic wind direction changes in large eddy simulations of wind farms

The wind direction in atmospheric boundary layers changes continuously due to meso-scale weather phenomena. Developing accurate simulations of these changes is essential for understanding their effect on the performance of large wind farms. Our study introduces a new technique to model dynamic wind direction changes obtained from meso-scale simulations or ﬁ eld measurements in micro-scale large eddy simulations. We propose a method in which the simulation domain is treated as a non-inertial rotating reference frame. The primary bene ﬁ t of our approach is that it is straightforward to imple-ment and reproduces desired wind direction changes excellently. We veri ﬁ ed our approach in neutral atmospheric boundary layers and show that the observed boundary-layer characteristics for dynamic wind directions agree very well with those observed for constant mean wind directions when the wind direction is changed slowly such that the ﬂ ow is quasi-stationary. Further, we show that atmospheric measurements of the wind direction can be reproduced by our method. To underline the importance of the method, we conclude with a representative scenario, which shows that dynamic wind direction changes can affect the performance of large wind farms.


Introduction
With the increasing size of wind farms, there is a growing need to understand how wind farm performance is affected by changes in atmospheric conditions. The effect of atmospheric phenomena with scales that are much larger than the typical size of wind farms is still terra incognita [1,2] and needs further exploration. In particular, changing wind directions can strongly influence the performance of wind farms as the wakes from upwind turbines can greatly affect the power production of downstream turbines, and this effect depends strongly on the wind direction [3e5]. While the effect of different mean wind directions on wind farm performance is well explored, the effect of dynamic wind direction changes, which originate from meso-scale weather phenomena, on wind farm performance is not well-understood and needs further investigation [6e9].
Meso-scale simulations in which dynamic large-scale wind direction changes are accounted for, are usually restricted to horizontal resolutions larger than the turbine diameter [10e14]. Different parameterizations have been developed to represent wind farms in these meso-scale models. Commonly used models include the use of increased surface roughness to parametrize the effect of a wind farm [15e17], or a more detailed approach in which momentum is extracted, and turbulent kinetic energy is added at rotor height as pro-posed by Fitch et al. [18]. However, the horizontal resolution in these models is extremely coarse, due to which the interaction between the individual turbines cannot be investigated [19,20].
These interactions are commonly studied in micro-scale large eddy simulations (LES) of wind farms [8,9,15,21]. However, a vast majority of these studies focus on small scale turbulence and consider cases in which the flow is forced to approach steady-state conditions [15,22,23]. Although these simulations provide great insights into the steady-state interaction of wind farms and atmospheric boundary layers (ABLs), they ignore the influence of largescale effects such as the influence of dynamic wind direction changes on wind farm performance.
To simulate more realistic inflow conditions meso-scale forcings have to be included in the LES [1,24e26,]. This can be achieved by nesting the LES within a meso-scale simulation domain, e.g. by coupling LES to meso-scale models like the Weather Research and Forecasting model (WRF) [24,27e34]. As meso-scale simulations do not resolve the turbulent structures up to the same scale as LES, this affects the LES modeling itself [35e38]. An alternative approach [6,7] is to represent the effect of meso-scale effects in the microscale simulation so that the LES modeling itself is not affected. The benefit of this approach is that it allows a direct comparison with LES results in which meso-scale effects are not included. This allows one to study the influence of these large-scale effects, such as dynamic changes in the wind direction, independently.
To model the effect of dynamic wind direction changes obtained from field measurements or meso-scale simulations in LES, Munters et al. [7] proposed to use a concurrent precursor method [39] in which they rotate the horizontally periodic precursor domain. Chatterjee et al. [6] modified the method of Munters et al. [7] and proposed to rotate the inflow velocity vector instead of the precursor domain. They used data from LIDAR scans to model the effect of dynamic wind direction changes on the operation of the Alpha Ventus wind farm. Both studies revealed that dynamic-wind directions can significantly impact wind farm performance.
However, rotating the precursor simulation requires a sequence of geometrical interpolations and significant MPI communication.
Here, we propose a more straightforward method to simulate dynamic wind direction changes in LES. Our approach is inspired by the use of a proportional-integral-derivative (PID) controller, which is used in wind farm simulations to keep the wind angle constant at a particular height [40e42]. We therefore treat the simulation domain as a non-inertial rotating reference frame, which is an attractive approach as it only requires small changes to the governing equations and is straightforward to implement. A major advantage of our approach is that it avoids the geometrical interpolations and associated computational overhead that are required in the previously considered methods [6,7]. Besides, the methods discussed above [6,7] require a concurrent precursor inflow technique, i.e. an additional concurrent simulation from which the inflow data for the wind farm simulation is sampled. This condition makes it impossible to perform simulations with periodic boundary conditions. Such simulations are, for example, used to perform simulations of infinite wind turbine arrays that are considered in the development of analytical wind farm models [15,43]. In contrast, we will demonstrate that our method can be applied with and without a precursor method, allowing simulations of both finite and infinite wind farms.
The remainder of this paper is organized as follows: a description of the governing equations used to model dynamic wind direction changes in LES is outlined in section 2. A validation of the approach for neutral ABLs is presented in section 3, and the method is applied to a representative scenario in section 4 in which dynamic wind direction changes affect the wind farm performance. The conclusions will be presented in section 5.

Rotation of the mean wind direction in LES
The simulations are performed using an updated version of the LES code developed by Albertson and Parlange [44,15]. The updated code has been successfully used to study neutral, stable, and unstable ABLs [45,39], as well as the flow dynamics in extended wind farms [46e48]. The governing equations are the filtered incompressible continuity and momentum equations. The aim is to include dynamic wind direction changes qðtÞ, which can be obtained by meso-scale simulations or field measurement data, in the LES. For this purpose, the reference frame is rotated with an angular velocity u ¼ 0:5v t qðtÞ and corresponding non-inertial forces are added to equation (2). Here, the factor 0.5 can be explained as follows: half of the Coriolis acceleration arises due to the relative velocity and half due to the turning of the frame of reference [49]. As a consequence, the wind direction relative to a fixed axis is changed by an angle qðtÞ in the time where the frame of reference rotates by an angle equal to 0:5qðtÞ (see, e.g., Persson [49]). The resulting equations are: Here, the tilde represents spatial filtering with a spectral cut-off filter at the LES grid scale D andũ i represents the filtered velocity field components. t ij ¼ g u i u j Àũ iũj is the trace-less part of the subgrid scale (SGS) stress tensor and it is modeled with a standard Smagorinsky model [50] using a constant Smagorinsky coefficient C s ¼ 0:16 [51]. The trace of the SGS stress tensor is absorbed into the filtered modified pressurep þ ¼p=r À p ∞ =r À t kk =3, note that p * is defined below. The force f i is added for modeling the effects of the wind turbines, which are parameterized using an actuator disk approach [15,52]. Since the simulations are performed at very high Reynolds numbers we neglect viscous stresses [53], which is a common practice in LES of ABLs. The wall shear stress at the ground is modeled using the Monin-Obukhov similarity theory [54,55]. We use a surface roughness of z 0 ¼ 0:1m. The boundary conditions at the top of the domain are zero vertical velocity and zero shear stress. The mean-pressure gradient v i p ∞ =r defines a reference friction velocity u * ¼ À ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Hp ∞ =r p . Length and time scales are nondimensionalized with the domain height H and H=u * , respectively. To present the results in dimensional form, we assume that the incoming wind velocity at z ¼ 150m is 10 m/s, which is representative for the typical value observed in offshore conditions, see the Dutch Offshore Wind Atlas [56]. From this we obtain that u * z0:5m=s [4], which means that each non-dimensional time unit corresponds to t dim ¼ H=u * ¼ 2000sðz0:56hÞ.
While the continuity equation (1) is rotational invariant, the momentum equation (2) is modified by adding the Coriolis force 2uũ j ε ij3 . Besides, the direction of the mean-pressure gradient, which is driving the flow, is aligned with the desired wind direction qðtÞ. In a non-inertial, rotating reference frame, in addition to the Coriolis force, the centrifugal force Àu 2 x i ðd i1 þd i2 Þ and the Euler force ε ij3 x j du=dt are introduced [57]. The centrifugal force can be rewritten as a conservative force À u 2 Once rewritten as a conservative force, the centrifugal force is combined with the pressure term: The Euler force is neglected here, as a periodic domain is considered and the distance to the axis of rotation, which is required for its calculation, is unknown. Time integration is performed using a second-order accurate Adams-Bashforth scheme. Derivatives in the vertical direction are calculated using a second-order central finite difference scheme. In streamwise and spanwise directions a pseudo-spectral method is applied. Thus, doubly periodic boundary conditions are considered in the horizontal directions, which implies that an infinite wind farm is considered [15]. To model finite size wind farms we employ a concurrent precursor inflow method [39]. In this approach we sample flow data from a periodic turbulent ABL simulation performed in a precursor domain. The sampled data is introduced as inflow velocity into a fringe region of the wind farm simulation domain. To ensure a smooth transition between the velocity in the wind farm domain u i;WF and the inflow velocity sampled from the precursor simulation u i;Pre a symmetric weighing function wðxÞ is applied in the fringe region: u i;Fringe ðx;y;z;tÞ¼wðxÞu i;Pre ðx;y;z;tÞþð1ÀwðxÞÞu i;WF ðx¼L start ;y;z;tÞ (3) where: The parameter L s sets the starting point of the fringe region. Here, we select L s ¼ L x À 0:2L x such that the length of the fringe Fig. 1 shows the corresponding weighing function and Fig. 7 shows the fringe regions in the wind farm simulation domain. However, before we employ our method to a simulation of a representative wind farm, we first test its performance in a neutral ABL in section 3.

Dynamic wind direction changes in LES of neutral ABL
In section 3.1 we validate the proposed method to model dynamic wind direction changes in LES for flows that are quasistationary. In section 3.2 we show that our method can be used to reproduce the dynamic wind direction changes obtained from atmospheric field measurement data.

Validation of the approach
We perform LES of a neutral ABL to validate the described method to incorporate dynamic wind direction changes. The selected size is L x ¼ L y ¼ 10km and L z ¼ 1km, L y , and L z are the domain length in the spanwise, and vertical direction, respectively.
The simulations are performed on a grid with 256 Â 256 Â 48 nodes. We perform three simulations in which the wind direction is changed linearly with time, i.e. qðtÞ ¼ q 0 t and consider different rotation speeds q 0 ¼ ½9 + =h; 27 + =h; 54 + =h. In addition, the following combination of sines and cosines: is considered to assess the performance for more or less random wind direction changes. These simulations are compared with a reference simulation with a constant wind direction. Fig. 2 confirms that the horizontal wind direction at hub-height CqDðtÞ ¼ tan À1 ðCvDðtÞ =CuDðtÞÞ follows this imposed wind direction excellently. It is worth noting that this excellent overlap would not be achieved by only varying the driving pressure gradient, because a large phase-lag (up to several hours [7]) is visible between the pressure gradient and the mean wind direction when the additional Coriolis force is neglected [4].
qðtÞ ¼ À1 + t þ ½15 + À 3 + sinð0:7t=5Þsinð2t=5Þ À 11 + cosðtÞsinðt=10Þ À15 + sinð2t=15Þ þ cosð1:5t=5Þ þ sinðt=2Þ À 3 + cosðt=4:5Þ; A visualization of the horizontal velocity magnitude u h ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Fig. 3. In the top row, the horizontal wind direction is rotated from The visualizations reveal streamwise-elongated coherent structures typically observed in neutral ABL simulations [21]. For the 0 + wind direction (left panel in Fig. 3), these structures are oriented parallel to the x-axis. When the mean wind direction is rotated the large-scale flow structures orient themselves with the mean flow direction, and we do not observe any unusual stretching of the flow structures due to the rotation. We also compare the time and horizontally averaged turbulent statistics obtained from the reference simulation to the simulation results in which the mean wind direction is dynamically rotated to validate the proposed method. The mean velocity u h for the different rotation rates is depicted in Fig. 4(a). For all cases, the velocity profiles agree well with the reference result. We observe only small differences in the velocity at the top of the domain. This difference could be caused by slight variations in the large-scale structures when the flow is rotated. When the wind direction changes, the domain length in the flow direction continuously changes, see Fig. 3, and the high and low-velocity streaks tend to adjust themselves to this. In addition to this inevitable domain effect, the neglected Euler force might also cause these small differences when comparing the statistics with the constant mean wind direction case. Fig. 4(b) shows that the time and horizontally averaged variance of the horizontal velocity for the different rotation rates agrees well with the reference case. All the cases show the same trend, exhibiting a maximum close to the ground. The maximum indicates the height up to which the influence of the surface friction dominates. At this height, the skewness of the horizontal velocity, displayed in Fig. 4(c), turns from positive, super-Gaussian to sub-Gaussian at higher positions. While there are small quantitative differences for higher rotation rates, the qualitative trend is consistent. This consistency is also present for the flatness presented in Fig. 4(d). The flatness increases with height above the surface layer, corresponding to an increase in rare but extreme deviations from the mean velocity. This qualitative trend is the same for all rotation rates. It is worth mentioning here that the higher-order statistics, such as skewness and flatness, provide a stricter validation of the proposed method than lower-order statistics. The vertical velocity is zero at the bottom and top of the for a constant mean wind direction (q ¼ 0 + ), linearly rotating wind direction (qðtÞ ¼ 9 + =h,t; 27 + =h,t; 54 + = h, t) and a randomly varying wind direction (see Fig. 2). b) Profiles of the variance, c) skewness and d) flatness of the horizontal velocity magnitude as function of z. domain and shows a maximum at zz0:2km, see Fig. 5(a). The figure shows that the variance obtained from the simulation with changing wind direction agree very well with the reference case. Furthermore, the higher-order statistics such as skewness and flatness ( Fig. 5(b)e(c)) only vary slightly with increasing rotation rate.
Overall, the quantitative and qualitative results agree for all rotation rates, which indicates that the flow characteristics remain the same when the wind direction is rotated slowly. This validates that the non-inertial rotating reference frame method has been implemented correctly.

Comparison with field measurement data
To assess the ability of our method to represent dynamic wind direction changes from meso-scale weather phenomena into micro-scale LES, we compare the simulation results to field measurement data. In Fig. 6(a) we reproduce the wind direction measurements taken from a wind vane at a height of 87 m on the M5 meteorological mast at National Wind Technology Center [58], of National Renewable Energy Laboratory (NREL), for a 7 h period on the 1st of February, 2018 from 6 a.m. to 1 p.m. Fig. 6(b) shows the corresponding power spectrum of the wind direction changes.
We performed LES of a neutral ABL in a domain of L x ¼ L y ¼ 5km; and L z ¼ 1km using a 128 Â 128 Â 48 and a 256 Â 256 Â 96 grid. Fig. 6 shows that LES with a constant mean wind direction captures the high-frequency wind direction changes fairly accurately, especially considering that we perform neutral ABL simulations, instead of matching the atmospheric conditions of the observational data. The figure shows that the higher resolution simulation captures the high-frequency wind direction changes better. However, the low-frequency wind direction changes are not represented in the LES that considers a constant wind direction.
The low-frequency wind direction changes can be included in the LES using the low-pass filtered field measurement data as input to the LES. Here, the low-pass filter's cutoff frequency is chosen as 0.0008 Hz. In agreement with the results in section 3.1, we find in Fig. 6(a) that the mean wind direction of the LES perfectly follows the desired wind direction. More importantly, Fig. 6(b) shows that the spectrum of wind direction changes obtained from LES now accurately represents the entire frequency range. As intended, our approach models the low-frequency wind direction changes obtained from field observations, or a meso-scale simulation, while it does not affect the high-frequency range, which we assume to be represented accurately by our micro-scale LES. In the next section, we will apply our method to a representative scenario to demonstrate that these low-frequency wind direction changes can significantly affect the performance of extended wind farms.

Case description
We perform LES of a symmetric wind farm with 6Â 6 turbines in a neutral ABL with the same surface roughness considered previously. Both the wind farm and the precursor domains have a size of L x ¼ L y ¼ 7:5km and L z ¼ 1km. The last 1.5 km of each horizontal direction is used as a fringe region. The wind farm layout is presented in Fig. 7. The simulations are performed on a uniform grid with 384 Â 384 Â 64 nodes and are used to demonstrate that dynamic wind direction changes can significantly affect wind farm performance.
To demonstrate this we consider the following representative scenario with sinusoidal wind direction changes: qðtÞ ¼ 20 + sinð2pt =T q Þ with T q ¼ 0:6h and T q ¼ 2:2h and 13 additional reference simulations in which a constant mean wind direction is considered, see Fig. 7. At each time step of the simulation, the wind turbines are rotated perpendicular to the local incoming wind direction to ensure that there is no yaw misalignment. We note that this instantaneous rotation of the disks is idealistic as wind turbines adjust their orientation with respect to the incoming wind direction with a time delay [59]. This can lead to additional yaw effects, which are not included in our representative scenario. The turbines have a thrust coefficient of C T ¼ 3=4, a diameter of D ¼ 150m, a hub-height of z hub ¼ 150m, and the distance to neighboring turbines is four turbine diameters in both horizontal directions. In the following, the total power production of the wind farm P tot is normalized by taking into account the velocity at hubheight and the power obtained for the reference case (P ref ) with constant mean wind direction of q ¼ 0 + . Fig. 8(a) and (b) present the time-variation of the imposed and the measured wind direction at hub-height for the two cases under consideration. Circular markers indicate when the wind direction is Fig. 8(c)e(d) depict the time-variation of the total power production of the wind farm. The circular markers denote the power at the time instant when q ¼ 0 + . Furthermore, the power production of two constant mean wind direction cases q ¼ 0 + and q ¼ 15 + are given as a reference. For time-varying wind directions, we observe that the power is mostly higher than for the case with a constant mean wind direction of 0 + at both slow and fast rotation rates, see Fig. 8(c)e(d). This is due to the strong wake-effects for the Fig. 8. a) & b) Imposed and measured wind directions at hub-height with T q ¼ 2:2h and T q ¼ 0:6h, respectively. Circular markers represent the reference wind direction

Hysteresis effects in wind farm power production
Corresponding normalized wind farm power production for T q ¼ 2:2h and T q ¼ 0:6h, respectively. q ¼ 0 + wind direction. Surprisingly, for fast wind direction changes (case T q ¼ 0:6h, see Fig. 8(d)), the minimum power is not reached for q ¼ 0 + . Instead, the power minima are reached after the wind direction passed q ¼ 0 + . A comparison with the constant mean wind direction q ¼ 15 + reference case reveals that the wind farm power production can be a bit higher due to the effect of the dynamic wind direction changes. Fig. 9(a)-(b) display the wind farm power production as function of the wind direction q. These results are obtained by binning the time-varying power production based on the instantaneous wind direction. A solid black line displays results that are binned over the entire simulation. Additionally, results are divided up into time periods during which dq=dt > 0 (red, dash-dotted line) and dq= dt < 0 (blue, dashed line), respectively. Each circle represents simulation results obtained with constant mean wind directions. Fig. 9(a) shows that the wind farm power production agrees well with the results obtained from the constant mean wind direction cases when the wind direction changes slowly ðT q ¼ 2:2hÞ Due to the symmetric layout of the wind farm, the power production is symmetric around q ¼ 0 + . The power production is lowest for q ¼ 0 + when the wind is aligned with the farm layout. The maximum inter-turbine spacing and thus maximum power production is reached at qz15 + : The wind farm power production is nearly independent of the sign of the wind direction change dq= dt. However, Fig. 9(b) shows that the wind farm power production depends on the sign of the wind direction change dq= dt for faster rotation rates. For dq=dt > 0 the power production agrees well with the values obtained for constant mean wind directions with negative q. In contrast, the power production is lower than for the constant mean wind direction cases for positive q and the minimum power production is observed for qz3 + . An exception is found between q ¼ 12 + and q ¼ 20 + ; where the power production is higher than for the constant mean wind direction cases. Due to the symmetric farm layout, a similar pattern is found for dq= dt < 0 with the symmetry axis positioned at q ¼ 0 + : These hysteresis effects can be explained by examining the development of the wake between the turbine rows, as displayed in Fig. 10 and the corresponding movie (see the supplementary materials). The figure shows flow snapshots of the horizontal velocity magnitude normalized by the inflow-velocity for different wind directions.
Each column represents one wind direction between q ¼ À9 + and q ¼ 9 + . The top row displays snapshots of the horizontal velocity magnitude from simulations performed with a constant mean wind direction. In the lower rows, instantaneous flow fields for the same range of wind directions are shown for the simulation in which the wind direction varies with a period of T q ¼ 0:6h. The difference between the middle and bottom rows is the direction in which the wind direction changes, i.e. from dq=dt > 0 in the middle row and dq=dt < 0 in the bottom row.
For dq=dt > 0 the minimum power production is lower than for the 0 + reference case and observed at qz3 + . When dq=dt < 0 the minimum is positioned at qz À 3 + . We note that also Munters et al. [7] found that the flow angle for which the minimum power production is observed can change when the wind direction changes dynamically. The observed effects for dq=dt > 0 can be explained by the low-velocity zones between the turbine rows, which originate from earlier time steps (see Fig. 10, middle row). Due to the fast rotation rate, the low-velocity zones between the turbine rows at x > 3 km did not mix with the incoming high-velocity inflow yet. Therefore, turbines at x > 3 km cannot entrain energy from the sides, which is possible for the constant mean q ¼ 0 + wind direction case for which high-velocity wind-speed channels are formed between the turbines, see Fig. 10. Besides, we speculate that the dynamic wind direction changes may influence the vertical kinetic energy flux that brings down high-velocity wind from above the wind farm to the hub-height plane. In previous work it was namely observed that wakes recover faster when their inter turbine distance is smaller, which leads to a relatively strong wake recovery for the aligned configuration [60]. Unfortunately, as the vertical kinetic energy flux cannot be conditionally sampled on the wind direction, we cannot verify this hypothesis at the moment. In the range of wind directions selected in this study ðq ¼ ½ À 20 + ; 20 + Þ; we observe that for dq=dt > 0 the power production at q ¼ 15 + is higher than for the corresponding mean wind direction case. The movie shows that the turbines at x > 3 km and y > 1:5 km benefit from the high-velocity wind speed zones between the turbines (see the supplementary movie). However, when dq=dt < 0 or when the mean wind direction is static at q ¼ 15 + , the turbines in this region of the wind farm are continuously in the wakes created by upstream turbines. The representative scenario considered above is presented to show that dynamic wind direction changes can significantly affect the performance of large wind farm. As we have shown, this can lead to the hysteresis effect in the power production. We emphasize that these effects will be wind farm-specific and will depend on wind farm design parameters and atmospheric flow conditions. For example, we would expect that the hysteresis effect will be more pronounced when the turbine thrust coefficient is higher. Besides, we expect hysteresis effects to increases with wind farm size in which longer length and time scales corresponding to the dynamic wind direction changes are more important than in smaller farms. Further studies are required to assess such effects in more detail.

Conclusions
We have presented a new technique to incorporate dynamic wind direction changes in LES of ABLs. The time evolution of the wind direction can be obtained from meso-scale simulations or field measurements. Our method is advantageous compared to previously considered methods [6,7] as our approach only requires small changes to the governing equations and is easy to implement. Besides, our approach can be applied to simulations of infinite wind farms, which is not possible with previously considered methods [6,7].
We performed neutral ABL simulations in which we varied the rotation rate of the wind direction to validate our approach. We find an excellent agreement between the imposed and simulated wind direction. We showed that the mean and higher-order flow statistics in the simulations with varying wind direction agrees very well with results obtained from a simulation with a constant mean wind direction when the flow direction is rotated slowly. Comparisons to measurement data demonstrate that our method produces a similar power spectrum of wind directions. This confirms that the non-inertial rotating reference frame is a good technique to model dynamic wind direction changes in LES.
Subsequently, we applied our method to a representative scenario to demonstrate some potential effects of dynamic wind direction changes on wind farm performance. We performed simulations for various wind directions and cases in which a sinusoidal wind direction variation (q ¼ 20 + sinð2pt =T q Þ with a time periods of T q ¼ 0:6h and T q ¼ 2:2h) is enforced. In agreement with previous studies [6,7], we show that dynamic wind direction changes can significantly affect the performance of wind farms. The presented demonstration case shows that dynamic wind direction changes can positively and negatively affect the wind farm power production.
However, we emphasize that further studies are required to better understand these effects. The observed hysteresis effect can, for example, depend on the wind farm design, atmospheric conditions, yaw misalignment with respect to the incoming flow direction, and the turbine thrust coefficient. In this work, we considered a neutral boundary layer situation, but we emphasize that there are no restrictions in extending the presented approach to stable and unstable boundary layer simulations. The present work focuses on modeling wind direction changes, but we note that other meso-scale phenomena like wind shear and temperature variations require further studies. CRediT authorship contribution statement

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.