An analytical solution for wind deficit decay behind a wind energy converter using momentum flux conservation validated by UAS data

The wind deficit behind a wind energy converter (WEC) results from a complex interaction of forces. Kinetic energy is removed from the atmosphere, but coherent turbulent structures prevent a swift compensation of momentum within the wake behind the WEC. A detailed description of the wake is beneficial in meso-scale weather forecast (e.g. WRF models) and numerical simulations of wind wake deficits. Especially in the near to intermediate wake (0.5− 5 rotor diameter D), the dominating processes characterising the wake formation change along the wake. The conservation equation of momentum 5 is used as a starting point to map the most important processes assuming the WEC operates at maximum efficiency in a neutral stratified boundary layer. The wake is divided into three different regions to accommodate the changing impact of atmospheric turbulence and the shear created by the WEC onto the wake. A differential equation that depicts the variable momentum transport into the wind deficit along the wake is derived and solved analytically. Additionally, a numerical solution (Euler method) of the simplified momentum conservation equation is shown to provide a quality control and error estimate to 10 the analytical model. The analytical solution is compared to in-situ wake measurements behind an Enercon E-112 converter, located in the Jade Wind Park near the North Sea coast in Germany, captured by the MASC-3 UAS (unmanned aircraft system) of the University of Tübingen. The obtained UAS data cover the distance from 0.5−5 D at hub height behind the nacelle. The analytical and numerical model are found to be in good agreement with the data of the three measurement flights behind the WEC. 15


Introduction
Wind energy converters (WECs) remove kinetic energy from the mean flow in the atmosphere which results in a wind deficit in the wake behind each converter. The interaction of these individual wakes with the atmosphere and even with neighbouring converters needs to be understood, for example in order to increase wind park efficiency. With increasing sizes of WECs and 20 wind farms the wind deficit created by the wind parks and individual converters affects the environment increasingly. The WIPAFF (Wind PArk Far Field) project determined the length of single wind park wakes and found their dimensions large enough to influence local wind conditions Bärfuss et al., 2019). Therefore, an influence on local weather (and neighbouring wind farms) is likely. WRF (weather research and forecasting model) simulations by Siedersleben et al. (2018) display the wind deficit using wind farm parametrisation (WFP). Emeis (2010) introduced a wind deficit model for wind farms, concluding that the wind deficit decreases because of momentum inflow from higher altitudes. In his Lagrangian approach, an air parcel is described, travelling from the end of the wind farm downstream in the wake for several kilometres. transports on a smaller scale in their respective wake area. A key parameter is to resolve the wind deficit as a main driver of shear stress in the wake. In numerical simulations wind deficit models based on the thrust coefficient of the WEC are common practice. For example Magnusson and Smedman (1999) describe the wind deficit as a time dependent exponential function. Bastankhah and Porté-Agel (2014) also propose a wind deficit model using the thrust coefficient of the energy converter. They also introduce the wake widening as an additional parameter to obtain a more precise prediction in the far wake. While the model by Bastankhah and Porté-Agel (2014) provides a spatial distribution for the wind deficit, it lacks atmospheric boundary conditions (turbulence, thermal stratification etc.) that have an influence on the wake development. Bastankhah and Porté-Agel (2017) refined their analytical model in a wind-tunnel, however, atmospheric conditions are hard to simulate in wind-tunnel 15 experiments. Usually, wind-tunnel results suffer from scaling errors when applied to the real world. Yet, for the simulated conditions (thermally neutral stratification and along the dimensions of the wind-tunnel) the model from Bastankhah and Porté-Agel (2017) works fine.
Hirth and Schroeder (2013) conducted a field measurement, using a mobile Doppler radar system to span a high resolution grid of radar measurements behind a WEC, studying the wind deficit influence of a neighbouring converter, in order to quantify 20 the power losses in wind farms.
In the present study in-situ wake measurements are used to display the wind deficit behind an Enercon E-112 converter in the near to mid wake (0.5 to 5 rotor diameters D). In-situ measurements have the advantage that the measured wind is already a superimposition of turbulence created by the WEC and the free stream turbulence. Consequently a model derived from UAS in-situ data needs to accommodate turbulence and shear terms in some way. Disadvantages of real-world experiments 25 are, of course, additional boundary conditions that are hard to quantify, e.g. is vegetation, additional WECs or obstacles of any kind (dyke, buildings) creating random, unrelated (to the investigated problem) turbulence that adds uncertainties to the measurement. Despite these issues, in the following, a steady state spatial analytical solution for a wind deficit is presented, based on the momentum flux conservation equation.
2 Theory and methods 30

Analytical solution of the conservation of momentum equation
The starting point for the presented analytical wind deficit model is the equation for conservation of momentum in the mean flow using Einstein summation notation (Stull, 1988). After Reynolds averaging and the Boussinesq approximation have been https://doi.org/10.5194/wes-2020-92 Preprint. Discussion started: 6 July 2020 c Author(s) 2020. CC BY 4.0 License. applied, the conservation of momentum can be written as: Here, i, j = 1, 2, 3 for all three directions in space, g is the gravitational acceleration, p the static pressure, f c the Coriolis parameter, ν the kinematic viscosity and ρ the density of air.

5
Term I represents storage of mean momentum. Term II describes advection of mean momentum by the mean wind. Term III allows gravity to act in the vertical direction only. Term IV describes the influence of the Coriolis force. Term V describes the mean pressure-gradient force. Term VI represents the influence of viscous stress on the mean motions assuming incompressibility. Term VII represents the influence of Reynolds' stress on the mean motions.
Term II can be rewritten for an incompressible flow (and thus divergence-free) (Etling, 2008): A one dimensional, horizontal steady state wind field with no buoyancy is assumed (no change in wind direction and wind speed at each location and time). No meso-scalic changes in the static pressure occur in the model. Thus, term I, III and term 10 V can be cancelled from the equation. In this study the wake is observed up to 5 D ≈ 600 m, as a consequence the Coriolis force (term IV) can be neglected. Term VI represents the viscous stress and observations in the atmosphere indicate that the molecular diffusion terms are several order of magnitudes smaller compared to the other terms and can be neglected (Stull, 1988). Equation 1 can now be simplified for u = u r , the reduced horizontal wind speed in the wake along the x direction, with the remaining terms II and VII, as: With the x-axis rotated into the main wind direction v becomes zero. Therefore term B can be cancelled from Eq. 3. For further simplification, only vertical momentum transport is assumed and lateral momentum transport is ignored, due to the absence of change in the lateral shear stress (term D and E can be cancelled). In a homogeneous unstable to neutral boundary layer subsidence can be neglected and therefore w = 0; term C can be cancelled (Stull, 1988). 20 ∂u 2 Equation 1 is now reduced to a much simpler differential equation in one dimension by cancelling insignificant parts, considering the order of magnitude of each term. The resulting Eq. 4 states that a change in the wind speed in x direction corresponds to The Reynolds shear stress can be expressed using a momentum transfer coefficient K m : In analogy to the derivation of a wind deficit model for wind parks by Emeis (2010) a vertical space ∆z is defined above hub 5 height z within which the wind speed is reduced to u r and where momentum is transported into the wake from above. Above ∆z + z an undisturbed flow velocity u 0 is assumed which varies only little with height (see Fig. 1). Thus the momentum flux within ∆z can be approximated by Eq. 5.
For the calculations presented later in this study, K m needs to be resolved: With u * the shear stress velocity, altitude (hub height) z and κ the van Kárman constant (Foken, 2006). Regarding the temperature profile, wind conditions and turbulence intensity (s.a Sec. 3.3), a stability parameter ζ = z/L of approximately 0 to 1 can be concluded, using Businger et al. (1971). Joffre et al. (2001) studied the variability of the stable and unstable atmospheric boundary layer height and documented a dependence of the shear stress velocity u * on the stability parameter ζ. From there, a value of u * ≈ 0.3 can be concluded for the respective atmospheric stability, present at the time of the field measurements.
Slight differences in u * solely shift the solution along the y axis. In Sec.5 a short sensitivity evaluation is presented. Now, that the Reynolds shear stress is parametrised, Eq. 4 can be re-written by making the occurring partials to differences and ∆u = u 0 − u r the difference between the undisturbed wind speed u 0 and the reduces wind speed u r in the wake, cf. Eq.
18 in Emeis (2010). After inserting Eq. 5 into Eq. 4 the differential equation Eq. 4 becomes: Rearranging and moving ∆x to the other side integrating on both sides yields: is a non-homogeneous non-linear differential equation (DE) of first order. The solutions of the corresponding homogeneous DE are of the form: which actually is vividly realistic: the reduced windspeed u r increases by the square-root of the horizontal distance x -quite similar to the vertical development of turbulent boundary layer or an internal boundary layer in heterogeneous terrain (Garratt, 1987(Garratt, , 1994Hanna, 1987). However, a particular solution of the non-homogeneous DE 7 could not be found. In order to find an approximated solution, u r is treated as a constant on the right-hand side of Eq. 8.1 and the simplified DE is solved in the following. The obtained result agrees sufficiently well with the numerical solution introduced later (s. Sec. 2.2).
A short assessment of the order of magnitude asures that the introduced error remains small. It is assumed that the change of the velocity gradient along x is small (with the velocity gradient in the order of u0−ur(0) 1000 m ≈ 7 1000 s −1 ). This complies with Taylor's hypothesis of frozen turbulence, similarly implemented in the model by Emeis (2010Emeis ( , 2017. This will introduce a small error, as afore mentioned, will be shown later in the results, in a comparison with a numerical solution of Eq. 8.1. Using the described simplification, the analytical solution is much more convenient to solve.
Rearranging again and multiplying with u r rewrite as a quadratic equation This quadratic equation has two solutions: While only the positive solution is physically relevant, resulting in positive wind speeds, c can now be determined using initial boundary conditions, i.e. for x = 0. These conditions can be measured or a theoretical value can be taken, e.g. u r (0) ≈ 0.3 u 0 using Betz law (Betz, 2013). The equation for the reduced wind speed along the wake can then be rewritten as: When examining a wind deficit created by a single converter, a closer look on the deficit formation is necessary. The different wake regions can be divided in three areas that, in numerical modelling, need different approaches (s.a. shear starts to drain turbulence and the wind deficit is beginning to be reduced (Medici and Alfredsson, 2006;Frandsen, 2007).
Note that the radial shear stress is represented by the vertical shear in this 2-dimensional approach (Eq. 4). In the intermediate wake, the wind deficit starts to decay and turbulence and momentum flux gradients begin to decrease. In the far wake (10 to 20 D), the wind deficit dissipates. The above mentioned wake regions are not rigid and depend on the mean wind speed, thermal stratification, surface roughness etc.

15
In the presented approach, the wind deficit of a single turbine is described by the vertical momentum flux. And since the near wake and the intermediate wake are of interest, a closer look onto the turbulence distribution inside these parts of a wake, is necessary. Frandsen (2007) describes how the wind deficit is reduced by ambient turbulence acting from the wake boundaries (s.a. green arrows in Fig. 2). The wake can now be divided in turbulence created by the wind deficit (s.a. purple volume in Fig.   2) and a mix of atmospheric and wake turbulence that has 'forgotten its origin' in the outer regions of the wake. Hence, in a 20 wake, the area inhabiting turbulence created by the wind deficit would be rotational symmetric volume that slowly becomes smaller as the distance to the nacelle increases (approximately up to 10 D).
Considering that the deficit decay depends on the distance x, the model will need to satisfy for this condition. Therefore, the frequency α (see Eq. 7) is modified to increase with increasing distance from the nacelle in the wake. In practice, this is implemented by an approximation of the radius of the core wake, that is still untouched by the free-stream turbulence (s.a.

25
Eq. 11) and is shown by the purple volume in Fig 2. For the one-dimensional approach in this study, this radius r(x) of this rotational symmetric volume can be relabelled to the height ∆z(x): For this study, with the available data, the distance at which the free-stream turbulence caves into the wake is approximately two rotor diameters (2 D = 4 R) (Medici and Alfredsson, 2006;Frandsen, 2007). Conclusively, a new function ∆z(x) is needed. prevents atmospheric turbulence (green arrows) to enter the wake up to about 2 D. After this point, the helical tip vortex structure begins to decay and atmospheric turbulence (green arrows) can enter the wake. The area inhabiting turbulence 'knowing its origin' from the WEC thins out with increasing distance to the nacelle (purple volume).
the kink in its analytical solution (Fig. 3). In the real world this boundary would be fluent (cf. the numerical solution in Fig. 8).
The 'dynamic alpha' analytical solution is computed in a semi-numerical way, this means that α(x) is calculated using Eq. 12 for each distance and then inserted into Eq. 10. A sole analytical solution is possible inserting Eq. 14 into Eq. 8.1. However, for this proof of concept with the Euler method as validation (s.a. Sec. 2.2), the more convenient and easy to implement method of the semi-numerical solution is used.

Euler method
The Euler method is a simple way to solve differential equations of first order numerically (Atkinson et al., 2009;Faragó, 2013). Thus, the this method solves the presented Eq 8.1 numerically. : The dash denotes a continuously differentiable derivation of a u(x). Discretising each calculation step with: x n = x 0 + h n, with n = 1, 2, 3, ...
The differential form is approximated by a discrete expression, with a discretisation step h = 0.1 m in the presented study.
Reducing h further does not change the result, but increases the computation time. Thus, the step size h = 0.1 m is used for 5 the numerical solutions. Both, the forward and backward method, are used to validate the analytical solution and to obtain a numerical solution for the difference Eq. 8.1, a short derivation of the specific Euler solutions is shown below.
Inserting Eq. 8.1 in Eq. 17 with ∆x = h and u = u r the Euler forward solution is: Re-writing Eq. 18 for the backward solution method yields: Multiplying both sides of the equation with u(x n+1 ) Rearranging to a quadratic equation gives: The resulting specific Euler backward solution is the positive version of the solution of the quadratic Eq. 20.4: 3 Data acquisition and data availability  Fig. 4). With its set-up it is possible to calculate the 3-D wind vector at any point in space up to 30 Hz. A more detailed description of the UAS and its instrumentation can be found in Rautenberg et al. (2019b) or Wildmann et al. (2013Wildmann et al. ( , 2014a.

Experimental site
The underlying in-situ data were captured behind an Enercon E-112 converter that is part of the Jade wind park north of  But as it turned out, in flight 1 the legs were too short to cover the whole wake downstream, thus in flight 1 the number of measurement legs is reduced to the ones in the near wake. Due to high wind speeds the UAS needed a larger turning circle as usual and the flight pattern needed adjustment (s.a green tracks in Fig. 5).

Atmospheric conditions
The measurement flights took place November 15 2019 at 15:30 LT in the Jade wind park with an easterly wind direction. An 10 average wind velocity of 10.5 m s −1 in hub height was measured. The turbulence intensity at this altitude was about 5 − 8 %. Figure 6 shows a vertical profile of the virtual potential temperature in the inflow that suggests near neutral conditions above 20 m a.g.l.

Results
In this section the MASC measurements are compared to two analytical solutions (one with constant and one with dynamic α as 15 seen in Eq. 12) and the numerical solutions (Eq. 17 and Eq. 21) of the simplified momentum flux conservation equation. It can be shown that Euler forward and backward method are identical. Both introduce a small error that increases with increasing ∆x. However, the analytical solution Eq. 10 is not far off from the numerical solution. In the presented data the maximum deviation between the analytical and numerical solution is ∆ max ≈ 5 % at 5 D. An additionally occurring parameters in the https://doi.org/10.5194/wes-2020-92 Preprint. Discussion started: 6 July 2020 c Author(s) 2020. CC BY 4.0 License. equations is the hub height of the E-112 converter. This hub height is set to z = 125 m. In the α = const. calculations, ∆z(x) is set to half a rotor diameter throughout the calculations.

Solutions using a constant α
The assumption of a constant α is the simplest approach, but neglects the measurements of Medici and Alfredsson (2006) and the theory by Frandsen (2007) who studied i.a. turbulence entrainment in the near wake. Figure 7 shows the analytical model

Solutions using a dynamic α
To represent the MASC data and to do justice to the entraining turbulence accelerating the wind deficit decrease, Eq. 10 and Eq. 8.1 need a slight modification. Now, Eq. 14 is used, thus, α defined in Eq. 12, is now a function of x. As described 5 above, ∆z corresponds to the radius in the wake of 'pure wake turbulence', still knowing its origin. With this parameter decreasing downwind of the wake, α increases along the x axis. This essentially creates a high turbulent region in which the turbulent momentum flux from aloft increases the decay of the wind deficit and conversely increases the mean horizontal wind. Figure 8 shows the result of an increasing α using ∆z(x) defined by Eq. 14 and shown in Fig. 3. The analytical and the numerical solution now follow the measured data. The measured data scatter about ±5 % around the computed solutions. The most important feature is the kink at x = 2 D, where both solutions are adjusted by considering the entraining free-stream turbulence, following the progressively decreasing wind deficit.

Far wake behaviour of the model
In the sections above the behaviour of the models is shown up to 5 D from the WEC. The far wake, roughly beginning at 8 − 10 D and beyond, is characterised by 'uniform' turbulence, meaning that atmospheric turbulence and turbulence created by 5 the WEC are undistinguishable from one-another. The additionally created turbulence by the WEC has almost dissipated and the wake stream is reintegrated into the surrounding flow. Figures 9 and 10 show the behaviour of the models with constant and dynamic α for the far wake. While the constant-α model underestimates the wake behaviour the dynamic-α approach follows the measured data up to 5 D and paints a reasonable picture of the wind deficit decay. Eventually, the numerical solution runs into saturation at u r /u 0 = 1.0 at about 12 D.

Discussion
The presented analytical and numerical solution acknowledge the near-wake turbulence behaviour. As shown above, the adjusted model (using a dynamic α) follows the measured data. In this first approach a rather simple hyperbolic function (Eq. 14) has been used to describe the WEC turbulence length that is not continuously differentiable, but there might be a smoother function that simulates the same behaviour.

5
The presented model can further be understood as a simple solution of the turbulent momentum conservation equation for a specific initial condition. The resulting equations are the solution of the simplified momentum flux conservation for the atmospheric flow that has been brought out of its equilibrium. The presented solutions do not represent the shear stress that is created by the WEC itself (solely by its presence in the flow). The solutions simulate the surrounding atmosphere refilling a 'velocity gap', created by a horizontal wind depression in flow direction (u r (0) = 0.3 u 0 ). The shear stress created by the 10 presence of the WEC, is automatically dealt with by the behaviour (steepness, trend) of ∆z(x).
The shear velocity u * is usually measured at the ground. In this field campaign no ground station was available. Thus, a value for u * had to be derived from the thermal stratification of the lower atmosphere. Although the MOST studies are widely accepted, a derivation of u * from thermal stability always contains uncertainties. Therefore, the influence of the shear velocity u * has been examined by changing the input value by ±50%. The resulting analytical solution then deviates by a maximum of approximately ±5% at 5 D in the resulting wind reduction, from the solution shown in this study. Thus, for example a u * of 0.45 still results in a good agreement with the data, while still be representative for a neutral to slightly unstable atmosphere (Businger et al., 1971).

6 Conclusions
The derived analytical solution shows good agreement with the in-situ UAS measurements. It is shown that in the wake, considering the free-stream turbulence distribution, and its inability to affect the wind deficit unevenly in the wake, is important for the development of the wind deficit. Although the lack of data for the far wake, the dynamic-α model connects the near wake with the intermediate and far wake turbulence and wind deficit. In this study a literature value for u * is used. In a future 10 measurement campaign, the shear stress velocity could be derived from two measurement legs in front of the WEC in different altitudes, taking advantage of the logarithmic wind profile. A direct measurement of the Reynolds stress u w at hub height could also be feasible. Also the wake development farther down the wake would be interesting to cover in future measurements, e.g. up to 10 − 15 D. in the wake is reached at ca. 10 − 12 D. Also the numerical model reaches its saturation of ur/u0 = 1.0 at around 12 D.
Regarding the far wake behaviour of the dynamic-α approach, it could be shown, that the analytical and numerical solution can predict an internal boundary layer. One of the main features of this approach (e.g. for an internal boundary layer model) is the capability to reach a certain saturation level. This can be a great advantage in modelling wind transitions from land to sea (or vice versa) or thermal internal boundary layers (TIBL) in general.
The model presented in this case can also be applied to larger scale problems. For example to describe the wind deficit of a 5 wind farm. Another interesting case is the additional consideration of the lateral turbulent momentum flux that must occur, as soon as the wind deficit disturbs the idealised neutral conditions of the free atmospheric flow (term E in Eq. 3).
Data availability. The measured UAS data can be provided by the authors (Moritz Mauz, Jens Bange, Andreas Platis) by ftp. Due to the amount and complexity of the data, we advise a brief introduction by one of the authors.