An analytical solution for wind deficit decay behind a wind energy converter using momentum 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 diameters D), the dominating processes characterising the wake formation change along the wake. The conservation equation of momentum is 5 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 the 10 analytical model. The analytical solution is compared to conventional WEC wake models and 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 (Stevens, 2016). With increasing sizes 20 of WECs and 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 the WEC using thrust coefficient based calculations. The model does not use WEC parameters as e.g. the thrust coefficient.
Instead it is assuming that the WEC is running in optimal condition and that Betz (1920) theory is applicable. The only WEC parameter going into the model is the rotor diameter D of the WEC.
In the present study in-situ UAS (unmanned aircraft system) based 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) and to measure the initial conditions 5 in the surrounding atmospheric flow. The MASC UAS platform was successfully used in previous measurement campaigns Rautenberg et al., 2019;Braam et al., 2016;Wildmann et al., 2015Wildmann et al., , 2014a. The obtained data is compared to the WEC wake models from Frandsen et al. (2006) or Porté-Agel (2014, 2017).
In the context of this study, 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 10 to accommodate turbulence and shear terms in some way. Disadvantages of real-world experiments are, of course, additional boundary conditions that are hard to quantify, e.g. vegetation, additional WECs or obstacles of any kind (dyke, buildings) creating random, unrelated (to the investigated phenomenon) turbulence that adds uncertainties to the wake measurement.
2 Theory and methods 2.1 Conventional engineering wake models 15 As a comparison two established analytical models are used with airborne measured data. The first one is the Frandsen et al. (2006) wake model where residual wind velocity u r normalised with the free stream velocity u 0 at hub height and at distance x in the wake of a WEC with rotor diameter D is given by: K is in the order of 10 k, with k the wake growth rate: and h the hub height of the WEC and z 0 the roughness height. The model uses a parameter β that relates to the thrust coefficient The second analytical model is a modification of the Frandsen model by Bastankhah and Porté-Agel (2014) derived from wind 25 tunnel measurements. The normalised residual wind speed in their model is described by: where σ is the lateral wake width. In this model a linear wake growth k is assumed.

The MaST model: a likely 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 5 flow using Einstein summation notation (Stull, 1988). After Reynolds averaging and the Boussinesq approximation have been 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. 10 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.
A one dimensional, horizontal steady-state wind field (i = 1) with no buoyancy is assumed (no local change in wind direction and wind speed at each location and time). No meso-scalic changes in the static pressure gradient occur in the model. Thus, term I, III and term V can be cancelled from the equation. In this study the wake is observed up to 5 D ≈ 600 m only, thus the 15 Coriolis force (term IV) can be neglected. Term VI represents the viscous stress and observations in the atmosphere indicate that the molecular diffusion is several order of magnitudes smaller compared to the other terms and can be neglected (Stull, 1988). Equation 6 can now be simplified. The mean reduced horizontal wind speed u in the wake at hub height in x-direction along the wake centre line for an incompressible flow, with the remaining terms II and VII, can be approximated as: Applying the product rule to each term on the left hand side of Eq. 7 yields: The left side of Eq. 7 can now be rewritten as: Since a steady-state incompressible flow is considered, the part in brackets in Eq. 11 is zero (divergence-free flow). The left side of Eq. 7 can now be written using Eq. 11.

5
Now, we focus only on the wake. Therefore we denote the residual wind speed in the wake along the centre line u = u r .
Finally Eq. 12 is the momentum conservation equation along x in differential form for an incompressible steady-state flow (Lecheler, 2011).
The downstream velocity in the wake increases with distance (∂u/∂x > 0). Considering the general continuity equation we 10 have: If furthermore a circular symmetric wake is assumed and comparable flow conditions below, left and right to the wake exist due to the initial flow conditions surrounding the wake with du 0 /dz ≈ 0 (s.a Fig. 1), then at the wake centre line the following is true: The left hand side (term A, B and C) of Eq. 12 at the centre line can now be simplified using Eq. 14: Note that due to the circular symmetry and the convergence towards the wake centre line the lateral velocity components must be zero there (v| y=0 = w| y=0 = 0). If this were not the case the whole wake would be moved away from the rotor axis of the 20 WEC and the centre line were no longer identical to the rotational axis of the WEC.
Coming now to the right hand side of Eq. 12 or the divergence of the turbulent fluxes of the longitudinal velocity perturbation u . Recalling the rotational symmetry of the wake and the initial conditions with du 0 /dz ≈ 0 the lateral flux perturbations (term E in Eq. 12) are identical to the vertical perturbations (term F). Concerning term D it is common, to neglect this longitudinal change (i.e. in x-direction) of the longitudinal turbulent flux of longitudinal momentum in comparison to ∂u 2 r /∂x. Equation 12 can now be rewritten: Equation 6 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. 16 states that a change in the wind speed in x direction corresponds 5 to a change of vertical momentum (in-)flux. This equation is the Eulerian analogy to Eq. 17 in Emeis (2010) using a Lagrangian approach.
To further simplify the approach, an empirical relation (Gradient method) is applied, assuming a first order distribution of the Reynolds shear stress along z. Therefore, the Reynolds shear stress can be expressed using a momentum transfer coefficient K m and can be brought to a finite difference form at the wake centre line at height h: In order to obtain the flux-divergence at height h for all x, assuming a symmetric wake flow above and below the centre line, the resulting form is: With these bulk simplifications and going from ∂ to d, Eq. 16 is now of the form: With ∆u = u 0 − u r the difference between the undisturbed wind speed u 0 and the reduced wind speed u r in the wake. The newly introduced parameter α is the wind deficit decay rate. For now, it will be treated as a constant. An additional solution with a changing rate α along x, will also be introduced and discussed. Due to the rotational symmetry ∆z will be renamed to the turbulence thickness d that describes the thickness of a volume that consists solely of turbulence created by the WEC. Through 20 this volume of thickness d the wind speed is reduced to u r and momentum is transported into the wake. Above z = h + d an undisturbed flow velocity u 0 is assumed which varies only little with height (see Fig. 1). C is a parameter derived from the model geometry, describing the momentum influx over the circular symmetric surface. This means for C = 8 momentum influx from all sides is equally pronounced. In reality this is not the case and will be discussed later. K m is the momentum transfer coefficient (s.a. Sec. 2.3).
To get an analytical solution, Eq. 19 needs now to be solved: Rearranging and moving dx to the other side integrating over du r on the left and over dx on the right yields: Equation 19 is the difference form of a non-homogeneous non-linear differential equation ( which actually is vividly realistic: the reduced wind speed 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 19 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. 20.1 and the simplified DE is solved in the following. The obtained result agrees sufficiently well with the numerical solution introduced later (s. App. A).
Additionally, a short assessment of the order of magnitude assures that the introduced error by treating u r as a constant for the integration step, is given. 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, a comparison with a numerical solution of Eq. 20.1 will be shown in the results. Using the described simplification, the analytical solution is much more convenient to solve, and the introduced error is small.
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. The initial conditions can be measured or a theoretical value can be taken, e.g.
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. Fig. 1): the near, 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 combined vertical and horizontal shear in this approach (Eq. 16). In the intermediate wake, the wind deficit starts to decay and turbulence and momentum flux gradients begin to decrease. In the far 10 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.
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 15 (s.a. green arrows in Fig. 1). The wake can now be divided in turbulence created by the wind deficit (s.a. purple volume in Fig.   1) and a mix of atmospheric and wake turbulence that has 'forgotten its origin' in the outer regions of the wake. Hence, in a wake, the area inhabiting turbulence created by the wind deficit, is a rotational symmetric volume that slowly becomes smaller, as the distance to the nacelle increases (approximately up to 10 D).
The model has to consider that the deficit decay depends on the distance x. Therefore, the wind deficit decay rate α (see 20 Eq. 19) is modified to increase with increasing distance from the nacelle in the wake. In practice, this is implemented by an approximation of the thickness of the core wake, that is still untouched by the free-stream turbulence (s.a. Eq. 22) and is shown by the purple volume in Fig 1. Equation 22 is a hyperbolic function. In a first approach this asymptotic function has been chosen to represent the remaining WEC turbulence along x. This is an educated guess, since the remaining turbulence (in this model approach) needs to decay from R to 0 in some manner. Then, for the one-dimensional approach in this study, the thickness of 25 this rotational symmetric volume is identical to a 'turbulence thickness' d(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 (Medici and Alfredsson, 2006;Frandsen, 2007 in its analytical solution (Fig. 2). In the real world this boundary would be fluent (cf. the numerical solution in Fig. 7). The 'dynamic alpha' analytical solution is computed in a semi-numerical way, this means that α(x) is calculated using Eq. 23 for each distance and then inserted into Eq. 21. A sole analytical solution is possible inserting Eq. 25 into Eq. 20.1. However, for this proof of concept with the Euler method as validation (s.a. App. A), the more convenient and easy to implement method of Figure 1. Sketch of a WEC wake, the incoming wind profile with u0 at hub height, and the wake turbulence development (not true to scale).
The coherent helical tip vortex structures (beige) 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 semi-numerical solution is used.

The momentum transfer coefficient K m
For the calculations presented later in this study, K m needs to be resolved: With u * the friction velocity at the surface, altitude z and κ the von Kármán constant Eq. 26 is valid in the surface layer (Emeis, 2010;Foken, 2006). For the sake of this derivation and the assumption that the error is small, it is assumed that Eq. 26 can be applied (s.a. Sec. 5). The friction velocity is computed from a vertical profile that was measured in front of the WEC in the inflow. Details are given in Sec. 3.3. In Sec.5 the sensitivity of the model towards u * is discussed.

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  Fig. 4). The E-112 was operating near its rated conditions.

Atmospheric conditions
The measurement flights took place November 15, 2019, at 15:30 LT in the Jade wind park with an easterly wind direction. 20 An average wind velocity of 11.0 m s −1 with a deviation of ±1.5 − 2 m s −1 at hub height behind the WEC was measured with gusts up to ≈ 13 m s −1 . The turbulence intensity (TI) at this altitude was about 5 − 8 %. It is computed by calculating TI = u hor /u hor from a vertical profile flown in the inflow of the WEC directly before Flight 1. With u hor being the horizontal wind. The averaging window is derived from the computed integral length scale (e.g. around hub height). Figure 5 shows a vertical profile of the virtual potential temperature in the inflow that suggests near neutral conditions above 20 m a.g.l. The vertical profile of the horizontal wind in front of the WEC is shown by black dots. The black solid line is the moving average 5 of the horizontal wind along 500 data points. While the temperature profile suggests a neutral atmospheric stratification the wind profile of the first 50 m resembles a profile that fits to a stable stratification. However, we assume this anomaly is caused by the vicinity to the dyke that reshapes the wind profile in the lower altitudes.
The friction velocity u * is calculated from the wind data up to 20 m a.g.l resulting in: 10 The measured friction velocity is used for a logarithmic fit of the wind data in front of the WEC using: u(z) = u * κ ln z z 0 , with κ = 0.4, u * = 0.29 m s −1 , z 0 = 0.0001 m and z the altitude over ground.
The resulting fit suggests a wind profile that is dominated by marine conditions (z 0 = 0.0001 m), which is realistic due to the vicinity to the North Sea and the easterly wind direction.

Data treatment and derived parameters
For the normalised residual horizontal wind speed u r /u 0 , both parameters have to be retrieved from the UAS measurements.
Therefore, for each flight leg two regions have to be defined: one for the undisturbed flow velocity and one that represents the decrease in wind speed due to wind energy conversion. To facilitate identifying these regions the measured data is smoothed out by a moving average of 200 data points (2 s). The undisturbed horizontal wind speed is then derived from the measurement 5 for each flight leg individually that is not influenced by the wake. The reduced wind velocity u r is derived from the minimum in the wake that results from the wind energy conversion. This method is used, since especially in the near wake, nacelle effects influence and increase the residual wind speed along the wake centre line. Figure 6 shows two exemplary wake measurements of the horizontal wind u hor at distances x/D = 2 and x/D = 5. The residual and free stream wind velocity are indicated by red and blue lines. Some measurements did not allow for a useful data evaluation, e.g. the aforementioned too short flight paths, unfinished turning manoeuvres of the UAS, data lags due to flow angles out of specification for the five-hole probe. Hence, all flights can lag data points at various distances in the evaluation. Some flight legs had to be discarded.
The model by Bastankhah and Porté-Agel (2014) allows for a measured wake growth rate k. This parameter is derived from the measurements shown in Fig. 6 and calculated to: The thrust coefficient C T of the E-112 WEC is derived from the operational conditions at the time of the measurement flights.
Therefore the tip-speed ratio (TSR) λ is calculated using a rotational velocity of 12 revolutions per minute (rpm), the rotor radius R and an average wind velocity of 11 m s −1 .

λ = ΩR
u 0 = 2π 12 · 57 m 60 s · 11 m s −1 ≈ 6.5 (30) 10 The TSR is now used to estimate the thrust coefficient using an C T -TSR chart of the very similar NREL 5 MW WEC (Al-Solihat and Nahon, 2018) and a pitch angle β = 1 • yielding a thrust coefficient C T ≈ 0.7. The same method has been applied previously in Mauz et al. (2019). All parameters that are used in this study are shown in Tab. 1. Note that the friction velocity used in the MaST calculations has been lowered to fit the data better, to counteract the dyke effect and to represent a typical value found in neutral atmospheric conditions.

4 Results
In this section the MASC measurements are compared to two conventional wake models and two analytical solutions (one with constant and one with dynamic α as seen in Eq. 23) and the numerical solutions (Eq. A3 and Eq. A7) of the simplified momentum flux conservation equation. It can be shown that Euler forward and backward method are identical. Both introduce 20 a small error that increases with increasing ∆x. However, the analytical solution Eq. 21 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 parameter in the equations is the hub height of the E-112 converter. This hub height is set to h = 125 m.
In the α = const. calculations, d(x) is set to half a rotor diameter throughout the calculations.
The analytical solutions are based on Eq. 21 which is a simplified solution of the differential equation 19. However, the constant C, appearing in Eq. 19 and in Eq. 23 needs adjustment when applying the model to real-world conditions. The constant C = 8 as derived above is only true, when considering that the momentum influx from all directions is equally pronounced (we 5 recall the rotational symmetry of the wake in the derivation of the model). Yet, in a real-world scenario the momentum influx from below is limited by the vicinity of the ground and usually lower wind velocities, in contrast to the constant downward influx from aloft from the free stream. A similar argument can be made for the lateral momentum influx that is not as quickly replenished from the sides as the vertical momentum influx that is constantly filling up the wake. Due to these reasons a value of C < 8 is more realistic. Following these arguments C has been halved due to the almost non-existent momentum influx from 10 below and halved again due to the almost non-existent lateral replenishment of momentum flux outside of the wake, yielding C = 2 which is used in all calculations.

The conventional wake models
In Fig. 7 the Frandsen et al. (2006) and the Bastankhah and Porté-Agel (2014) model (with measured k) do not represent the measured data very well. Both models overestimate the wind deficit decay in the near wake but might represent the far wake 15 from 5 − 6 D onwards adequately. However, there is no data to support this region sufficiently.
The model by Bastankhah and Porté-Agel (2014) can be fitted by lowering the wake growth rate k and discarding the measured value. Then the model describes the wake from 3 D onwards quite well, which is also the minimum distance that Bastankhah and Porté-Agel (2014, 2017) claim for their model.

MaST solutions using a constant α
The assumption of a constant α is the simplest approach, but neglects the measurements of Lignarolo et al. (2014); Medici and Alfredsson (2006) and the theory by Frandsen (2007)

MaST solutions using a dynamic α
To represent the MASC data and to do justice to the entraining turbulence, accelerating the wind deficit decay, Eq. 21 and Eq.
20.1 need a slight modification. Now, Eq. 25 is used, thus, α defined in Eq. 23, is now a function of x. As described above, d 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

Far wake behaviour of the MaST model
In the sections above, the behaviour of the models is described up to 5 D from the WEC. For the far wake there are no data available, due to legal limitations and flight regulations at the measurement site. In general, the far wake of a single WEC, roughly beginning at 8 − 10 D and beyond, is characterised by 'uniform' turbulence, meaning that atmospheric turbulence and 25 turbulence created by the WEC are theoretically undistinguishable from one-another. The additionally created turbulence by the WEC has almost dissipated and the wake stream is reintegrated into the surrounding flow. Figure 7 shows the behaviour of all models including the MaST solutions for constant and dynamic α up to the far wake. While the constant-α MaST model underestimates the wake behaviour, the dynamic-α approach follows the measured data up to 5 D. Eventually, the numerical solution runs into saturation at u r /u 0 = 0.95 at about 10 D.

Discussion
The presented analytical and numerical solution for the MaST model shown in Sec. 4 acknowledge the near-wake turbulence behaviour. As shown above, the adjusted model (which uses a changing α along the wake) follows the measured in-situ data.
The far wake behaviour of the dynamic-α model paints a reasonable picture of the wind deficit decay, indicating a wake length of roughly 10 D which is a common wake length in neutral conditions (McKay et al., 2012). This solution facilitates, in a first 5 approach, a rather simple hyperbolic function (Eq. 25) to simulate the decrease of WEC turbulence along x. A linear function (or decrease) has been discarded, because this would require a distinct intersection point at some distance x/D. Consequently, a distance would have had to been defined where ∆z = 0. Such an approach would have contradicted the idea of a dynamic α. As of now, the resulting function that is used to describe the remaining untouched WEC turbulence (with unit length) is not continuously differentiable. However, in this semi-numerical implementation this is not necessary. Also there most certainly 10 exists a smoother function that simulates the same behaviour which could be a task for a future evaluation.
The presented model can further be understood as a simple solution of the momentum conservation equation for a specific initial condition. The resulting equations are the solution of the simplified momentum conservation equation for the atmospheric flow that has been brought out of its equilibrium. The presented solutions also represent the shear stress that is created by the WEC itself (solely by its presence in the flow). We recall the C parameter in Eq. 19.
The likely solution of the simplified Navier-Stokes equation simulates the surrounding atmosphere refilling a 'velocity gap', 5 created by a horizontal wind depression in flow direction (u r (0) = 0.3 u 0 ). In the mathematical derivation of the model a rotational symmetric wake has been assumed, which is realistic. But also a rotational symmetric momentum influx has been assumed, which is helping the mathematical derivation, but is not too realistic, and leads to C = 8 in Eq. 19. To prevent an over-pronunciation of momentum influx C = 2 has been used throughout this study. The results show that even a value of 1 < C < 2 could provide a better fit. 10 To implement a shear or friction velocity Eq. 26 has been used. This equation is only valid in the surface layer (or Prandtl layer) and could be applied up to the inversion height z i . The assumption was made that the WEC was still in the surface layer in order to derive the model. The authors are aware that the validity of the equation is stretched in this work. As a consequence K m might be over estimated (by using a too large, out of scope value for z), thus, ∆z is under estimated (e.q. in Eq. 23). The application of Eq. 26 could be circumvent by using a proper ground measurement of u * , e.g. by using an EC (eddy-covariance) 15 station. The friction velocity was only available from a vertical profile, flown by the UAS in front of the WEC. The influence of the shear velocity u * has been examined by changing the input value by ±50%. The resulting analytical solutions then change their slope which then results in different wake lengths. The calculated wake lengths differ around ±3 − 4 D. Longer wakes are calculated for stable atmospheric conditions and smaller wake lengths for unstable conditions. For the sensitivity calculation the distance where atmospheric turbulence could act on the wake deficit was kept at 2 D. In reality this distance would be a 20 function of the tip-speed ratio λ (tip-speed of blade normalised by the horizontal wind speed) and the atmospheric thermal stability.
The second parameter influencing the wind deficit decay rate α is the separation height ∆z or in this study later called the turbulence thickness height d. In Emeis (2010Emeis ( , 2017 ∆z is the distance from hub height to the layer above the wind park wake which introduces vertical turbulent momentum flux into the wake, resulting in a wind deficit decay. In this study, only a single 25 turbine is considered. Instead of speculating about the dimensions of the separation height, an unquestionable boundary, the boundary between the wake centre line (deficit) and the undisturbed flow, is chosen to be represented by d. The distance d has an additional meaning in this study: it represents the radius of the WEC turbulence, measured from the wake centre. In the above paragraph K m and u * are discussed. Together with d those parameters control the wind deficit decay rate α = C·K m /d 2 .
For fitting purposes α can be treated as a black box, thus, it does not matter, if K m is overestimated and d underestimated. 30 Consequently, the model could be fed with measured values at hub height (e.g. for K m ) to study the separation height ∆z in the future. But in the present study, empirical relations and most probable values for parameters derived from secondary results (virtual potential temperature profile) were used to quantify α.
The wind vector calculation relies on a constant relative true air speed (TAS) of the UAS of 18.5 m s −1 . When leaving and entering a wake, the TAS does change quickly. The error in wind velocity calculation is around 10% (Rautenberg et al., https://doi.org/10.5194/wes-2021-21 Preprint. Discussion started: 15 March 2021 c Author(s) 2021. CC BY 4.0 License. 2019a). Since in this study the residual wind speed in the wake is a normalised parameter and the undisturbed free stream velocity suffers the same error, the error can be assumed as insignificant.
The conventional WEC wake models, as shown also in this study, primarily use operational parameters like the thrust coefficient, pitch angle, and wake growth rate etc. to determine the momentum conversion of the WEC and the resulting initial wind deficit. From thereon, usually an empirical exponential decay or square-root function is used to describe the wind deficit.

5
The atmosphere and its properties are not considered in these models. The analytical solution proposed in this study tries to provide a more profound model based on the conservation of momentum in the atmosphere. Eventually, both model approaches could be combined in the future.

Conclusions
The derived analytical solution shows good agreement with the in-situ UAS measurements and plausible behaviour in the far 10 wake. It is shown that in the wake, that the free-stream turbulence is unable to affect the wind deficit in the near wake. In the mid wake atmospheric turbulence and momentum affect 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 future UAS measurements the wake development farther down the wake would be interesting to cover as well, e.g. up to 10 − 15 D.

15
Regarding the far wake behaviour of the dynamic-α approach, it could be shown, that the analytical and numerical solution can predict a capped 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 20 wind farm by altering Eq. 12, which is investigated in an additional study.
In the future the presented model could be combined with conventional existing momentum sink models to get a description of operating WEC and its surrounding flow.
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.
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. 20.1, a short derivation of the specific Euler solutions is shown below.
Inserting Eq. 20.1 in Eq. A3 with ∆x = h and u = u r the Euler forward solution is: While Eq. A3 is called the Euler forward or explicit method, also an Euler backward or implicit method (Eq. A5) is calculated, 10 where the solution is only available implicitly, hence the name. While the Euler forward method is straight forward and simple to solve, for the Euler backward solution Eq. A5 has to be solved for u(x n+1 ). Note the appearance of u(x n+1 ) on both sides of the equation.
u(x n+1 ) ≈ u(x n ) + h f (u(x n+1 )) (A5) Re-writing Eq. A4 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. A6.4: Author contributions. MM evaluated the UAS data and derived the formalism. BvK supported the methods with his expertise and his ABL knowledge. AP added to the discussion and understanding of the topic. SE and JB shaped the manuscript and added valuable input to form