Calibration and validation of the Dynamic Wake Meandering model Part I: Bayesian estimation of model parameters using SpinnerLidar-derived wake characteristics

In this first part of two, we study the calibration of the Dynamic Wake Meandering (DWM) model using high spatial and temporal resolution SpinnerLidar measurements of the wake field collected at the Scaled Wind Technology Facility (SWiFT) located in Lubbock, Texas. We derive two-dimensional wake flow characteristics including wake deficit, wake turbulence and wake meandering from the lidar observations under different atmospheric stability conditions, inflow wind speeds and downstream distances up to five-rotor diameters. We then apply Bayesian inference to obtain a probabilistic calibration of 5 the DWM model, where the resulting joint distribution of parameters allows both for model implementation and uncertainty assessment. We validate the resulting fully-resolved wake field predictions against the lidar measurements and discuss the most critical sources of uncertainty. The results indicate that the DWM model can accurately predict the mean wind velocity and turbulence fields in the far wake region beyond four rotor diameters, as long as properly-calibrated parameters are used and wake meandering time series are accurately replicated. We demonstrate that the current DWM-model parameters in the IEC 10 standard lead to conservative wake deficit predictions. Finally, we provide practical recommendations for reliable calibration procedures.


Problem statement
As described above, there is no consensus for the values of the DWM model parameters when studying load predictions at any given site. Also, and perhaps most importantly, we do not know the sources of uncertainty observed in previous studies that used 95 the model (Larsen et al., 2013;Churchfield et al., 2015;Reinwardt et al., 2018), which need to be addressed to provide reliable load predictions. The common practice has been to derive optimized sets of model parameters based on limited synthetic or experimental data. This has lead to an unknown confidence in the overall model prediction ability; incorrect calibration of the model parameters may impact significantly the model performance and lead to suboptimal wind turbine designs. To address this issue, we estimate uncertainties in the calibration parameters of the DWM model by applying Bayesian inference (Box and 100 Tiao, 1973), which consists in updating any related prior information on model parameters by incorporating new knowledge obtained from wake flow characteristics derived through lidar measurements. Further, the Bayesian calibration provides a systematic approach to include various types of uncertainty such as physical variability as well as measurement and modeling errors. This paper is the first part of a two-part study dedicated on improving and validating the calibration of DWM model parameters using lidar-derived data. This first part has a four-fold primary purpose: 105 1. Derive wake flow features such as the two-dimensional velocity deficit and wake-added turbulence profiles, as well as time series of the wake meandering in both lateral and vertical directions from the SpinnerLidar measurements under different inflow wind speeds and atmospheric stability conditions. 2. Calibrate the DWM model-based wake deficit and wake-added turbulence predictions using the SpinnerLidar-derived wake flow features and the Bayesian inference framework. 3. Propagate modeling uncertainties in fully-resolved wake flow fields for robust predictions that take into account the calibrated uncertainties. 4. Conduct a sensitivity analysis to determine the most significant sources of uncertainty in simulated wake fields that are typically inputs to aeroelastic simulations.
This study contributes to the ongoing discussion regarding the accuracy of power and load predictions of wind turbines 115 operating under wake situations (Conti et al., 2020c,b), by quantifying uncertainties in wake simulations performed with the DWM model under a variety of inflow wind conditions. The outcomes of this study are useful for improving currently adopted wake simulation procedures for load analysis in the IEC standards, as well as to provide practical recommendations for wake model calibration studies based on measurements from nacelle-mounted lidars.
The work is organized as follows. Section 2 describes the DWM model. The SWiFT layout and relative wind site conditions 120 are described in Sect. 3. In Sect. 4, we present the wind field retrieval assumptions used to derive wake features from Spin-nerLidar measurements. The Bayesian calibration of the DWM model is performed in Sect. 5. We carry out the validation of the wind turbine wake simulations and conduct a sensitivity analysis to investigate the most influential parameters in Sect. 6.
Finally, the last two sections are dedicated to the discussions and conclusions.
The DWM model resolves three main wake features: the quasi-steady velocity deficit, the wake-added turbulence and the wake meandering. Each model component is described separately in the following subsections.

Quasi-steady velocity deficit
The quasi-steady velocity deficit component describes the wake expansion and recovery caused partly by the recovery of the rotor pressure field and partly by turbulence diffusion moving farther downstream of the rotor (Larsen et al., 2013). The wake 130 deficit is formulated in the meandering frame of reference (MFoR), which is a coordinate system with origin in the center of symmetry of the deficit.
In the far-wake region, i.e., distances larger than two rotor diameters (Sanderse, 2015), the deficit evolution is assumed to be governed by turbulent mixing and is described by the thin shear layer approximation of the rotational symmetric N-S equations with the pressure term disregarded (Madsen et al., 2010). To account for the neglected pressure gradient effects, an 135 initial wake deficit is analytically formulated based on the turbine's axial induction derived from blade element momentum (BEM) theory (Madsen et al., 2010). The turbulence closure of the N-S equations is obtained by means of an eddy viscosity term, and the momentum equation is solved numerically using a finite difference scheme with the artificial initial deficit as boundary condition (Madsen et al., 2010). Here, we use the numerical scheme of the standalone DWM model (Liew et al., 2020;Larsen et al., 2020). We refer to the generalized definition of the non-dimensional eddy viscosity term by Keck et al. 140 (2012), who considered two major drivers to the turbulence mixing: the ambient turbulence (T I amb ) and turbulence induced by the wake shear layer: where ν T is the eddy viscosity, U amb is the ambient wind speed at hub height, and R is the rotor radius. The first term to the right-hand-side of Eq. (1) describes the contribution of the ambient turbulence and the second the self-generated turbulence by 145 the wake shear layer. Madsen et al. (2010) proposed the instantaneous wake radius R w (x), wherex is the downstream distance normalized by R, and the maximum velocity difference (U amb − U min ), where U min is the minimum wind speed in the wake, as the turbulent length and velocity scales that govern turbulent mixing due to the wake shear layer.
Based on classical mixing length theory, Keck et al. (2012) defined the turbulence stresses to be proportional to the local velocity gradient ∂U (x, r)/∂r, which provides a two-dimensional eddy viscosity formulation that is function of the axial 150 and radial coordinates,x and r, respectively. The max operator is included to avoid underestimating the turbulent stresses at locations where the velocity gradient of the deficit approaches zero. Both terms in Eq.
(1) include a filter function (F 1 (x) and F 2 (x)) and a model constant (k 1 and k 2 ). The filter functions are required to model the turbulence development behind the rotor and have values in the range 0-1 depending on the downstream distance only (Keck et al., 2012). F 1 accounts for the delay of the ambient turbulence entrainment into the wake and is assumed to 'activate' ambient turbulence effects at downstream 155 distances where the pressure has recovered (≈ 2D, where D is the rotor diameter (Sanderse, 2015)). F 2 compensates for the initial non-equilibrium between the mean velocity field and the turbulent energy content created due to the rapid change in mean flow gradients close to the rotor. We refer to Eqs. (17) and (18) in Keck et al. (2015) for the mathematical formulation of F 1 and F 2 . k 1 and k 2 are calibration parameters that govern the turbulence mixing and presumably do not change with wind turbine design and ambient conditions. 160

Wake turbulence
The wake turbulence is composed of three turbulence sources and can be defined as (Vermeer et al., 2003): where T I m denotes the turbulence induced by the meandering of the wake deficit and T I add is the wake-added turbulence. T I m is commonly denoted as the apparent turbulence (Madsen et al., 2005), as the stochastic meandering of the wake deficit induces 165 additional velocity fluctuations into time series taken at fixed locations in the wake. This term is considered the main source of added turbulence in the far-wake (Madsen et al., 2010), while its spatial distribution can be computed by the convolution of the wake deficit in the MFoR, and the probability distribution function (PDF) of the wake meandering in the lateral and vertical directions (Keck et al., 2014a). T I add accounts for the shear-and mechanical-generated turbulence due to blade tip and root trailing vortices. The inhomogeneity of the wake-added turbulence is modeled by scaling the local turbulence using the factor 170 k mt (Madsen et al., 2010) as: where U def,M F oR is the velocity deficit in the MFoR, and k m1 and k m2 are constants calibrated based on CFD results (Madsen et al., 2010). The wake-added turbulence derived from Eq. (3) is presumed to meander together with the wake deficit, thus being displaced by the large-scale eddies in the atmosphere.

Meandering model
Here, the meandering model is confined to a single wake scenario, whereas multiple wake dynamics are described in Machefaux (2015). Wakes are considered to act as passive tracers driven by the large-scale atmospheric turbulence structures. The wake field is modeled by considering a cascade of consecutive wake deficits that are displaced by the large-scale lateral and vertical velocity fluctuations, i.e., the wake transport velocities (v c and w c ), corresponding to the lateral (y) and the vertical axis 180 (z), respectively. Adopting Taylor's hypothesis, the downstream advection of these deficits is assumed to be controlled by the mean wind speed of the ambient wind field. Larsen et al. (2008) estimated v c and w c by low-pass filtering atmospheric turbulence fluctuations. They defined a filtering cut-off frequency f cut,of f = U amb /(2D) thus excluding contributions from smaller eddies to the meandering dynamics. This assumption was verified using full scale lidar-based measurements collected behind an operating turbine . The wake displacements are computed as: wheret = x/U amb defines the time for an air particle to move from the rotor to the downstream distance in the wake region.
An appropriate choice of the transport velocity of the wake advection lies between the ambient wind speed and the centre velocity of the wake deficit (Keck et al., 2014b;Machefaux et al., 2015). The contribution from the yaw misalignment, which can redirect wakes in the lateral direction, is accounted for by h yaw (x,t) = x tan(θ(t)), where θ(t) is the yaw offset at the 190 specific time (Machefaux, 2015;Vollmer et al., 2016). The contribution of the rotor tilt is considered by h tilt (x,t).

The SWiFT facility
The SWiFT facility is a research site located in Lubbock, Texas, operated by Sandia National Laboratories (Herges et al., 2017). The site includes three Vestas V27 wind turbines, two meteorological towers, and a SpinnerLidar (Peña et al., 2018) mounted on the nacelle of one of the turbines and looking backwards. The entire site is on a fiber optic data acquisition 195 and control network that synchronizes recordings from masts, turbines, and the SpinnerLidar (Herges et al., 2017(Herges et al., , 2018. The measurement campaign took place between 2016 and 2017 with the main objective of characterizing wake fields and investigating wake steering control strategies (Herges et al., 2017).
Figure 1 provides an overview of the test layout together with the notation used along the manuscript. In this study, we analyze data collected at the meteorological mast (METa1), the turbine (WTGa1), and the SpinnerLidar mounted on the nacelle 200 of the WTGa1. The METa1 (hereafter referred as the mast) is 60-m tall and instrumented with sonic anemometers at 10, 18, 32, 45, and 58 m, sampling at 100 Hz. Other instruments installed on the mast are reported in Herges et al. (2017). The mast is placed 2.5 D south of WTGa1 in compliance with the IEC standard guidelines (IEC, 2015(IEC, , 2017. As southerly winds are prevalent at the site (see Fig. 1-right), this layout allows to retrieve concurrent incoming wind conditions from METa1, wake measurements behind the WTGa1 performed by the SpinnerLidar, and power and load measurements on the waked-WTGa2 205 installed 5D downstream. The WTGa1 and WTGa2 are variable-speed and pitch-regulated turbines with hub height of 32.1 m, D = 27 m, and a maximum power output of 192 kW (Herges et al., 2018). The supervisory control and data acquisition (SCADA) is available for both turbines providing records of the rotor speed, pitch and yaw angles, and power production, among others, at 50 Hz.

210
The SpinnerLidar is a research Doppler wind lidar developed at DTU based on a continuous-wave (CW) laser system (Peña et al., 2018). Hereafter, SpinnerLidar and lidar denote the same system. The SpinnerLidar has been mounted either in the spinner or on top of the nacelle of a wind turbine (Angelou and Sjöholm, 2015;Peña et al., 2018). The SpinnerLidar scans the rotor wake at high temporal and spatial resolution so that wake features can be derived. For the SWiFT campaign, the SpinnerLidar scanned continuously in a rose-pattern every 2 s (see Fig. 2), and the system internally subdivided the rose into 215 984 sections. The accumulated Doppler-shifted spectra at each of the sections was also recorded (Herges et al., 2017).
Once a scan was completed, the SpinnerLidar refocused at a different range and this process took about 2 s (Herges et al., 2018). For the SWiFT campaign, few scanning strategies were adopted that are described below: -Strategy I: the SpinnerLidar scanned seven downstream distances: 1, 1.5, 2, 2.5, 3, 4, and 5 D. A full cycle (i.e., from 1 to 5 D) took 30-42 s. This dataset is suitable for investigating the wake deficit evolution and recovery behind the rotor; 220 however, the frequency is too low to properly derive turbulence estimates or meandering dynamics.
-Strategy II: the SpinnerLidar scanned at the fixed distance of 2.5 D ensuring both high spatial and temporal resolution. ≈ 298 rosette scans were generated within a 10-min period. This dataset is suited for turbulence and meandering investigations.
-Strategy III: the SpinnerLidar scanned at the fixed distance of 5 D behind the rotor, generating about 298 scans each 225 10-min. During this period, power and load measurements were recorded on WTGa2. This dataset is suitable for load validation analysis. Since it provides a description of the wake flow field, including velocity deficits, turbulence and meandering at a distance that corresponds to typical spacings in wind farms, it is a valuable dataset for validating fullyresolved wake flow predictions as long as induction effects are accounted for.

Site conditions 230
For extended periods of the campaign, WTGa1 operated under large yaw misalignment, as wake steering strategies were being investigated (Herges et al., 2017). To consider periods where WTGa1 is nearly aligned with the mean inflow, we filtered out 10-min periods characterized by an average yaw offset larger than ±10 • compared to the free-stream wind direction (Conti et al., 2020a). Further, we focus the analysis on periods for which the free-stream wind direction is within 90 • -270 • (thus South winds, see Fig. 1 -right). This leads to about 850 available 10-min periods.  α is computed from the sonic measurements at 18 and 45 m. As shown, the site is characterized by a wide range of turbulence and shear conditions, which are consequence of the varying atmospheric stability Conti et al., 2020a).
Further, relatively low wind speeds are recorded (3-10 m/s); thus WTGa1 operates below rated power as seen in Fig. 3 (c).

240
Because of this range of operating conditions, high rotor thrust coefficients that induce strong wake deficits characterize this dataset.

Atmospheric stability
Here, we investigate the variability of the wake flow characteristics under varying stability and inflow wind speed conditions.
We classify each 10-min sonic-derived statistic into atmospheric stability classes defined by ranges of the dimensionless sta-245 bility parameter (z/L), where L is the Obukhov length (Monin and Obukhov, 1954) computed from the sonic measurements as: where u * is the friction velocity, k = 0.4 is the von Kármán constant, g is the acceleration due to gravity, T is the mean surfacelayer temperature, the vertical velocity component is denoted by w, and Θ is the potential temperature (which we approximate (−0.2 < z/L < 0.2), and stable (0.2 < z/L < 2) atmospheric conditions. We use the measurements at the 18-m sonic to derive the stability within each 10-min period.

Data statistics 255
The statistics of the inflow wind parameters are presented in separate Tables 1, 2 and 3, according to the relative SpinnerLidar scanning strategy. Table 1 presents data collected during strategy I. There is a sufficient amount of 10-min periods to characterize the variability of the wake deficit with respect to atmospheric stability, inflow wind speeds and downstream distances. The For Strategy II, represented in Table 2, less 10-min values are found; however this is sufficient to characterize wake turbulence and meandering under different stability conditions. For Strategy III, represented in Table 3, the dataset is characterized by stable conditions mainly, as the records correspond to night hours within three consecutive nights in July 2017 mainly.

Lidar measurements processing 265
As lidars only measure the line-of-sight (LOS) velocity (v los ), assumptions are needed to reconstruct the three-dimensional wind field u = (u, v, w), where u is the longitudinal, v the lateral and w the vertical velocity component. If we neglect any probe volume averaging along the beam, v los depends on the unit directional vector n = (cos φ cos θ, cos φ sin θ, sin φ), which describes the scanning geometry through the elevation (φ) and azimuth (θ) angles, and the wind field u, v los (φ, θ) = u cos(φ) cos(θ) + v cos(φ) sin(θ) + w sin(φ).
270 Table 1. Dataset from Strategy I. The data are classified according to wind speed bins of 1 m/s and three atmospheric stability classes; stable (s), near-neutral (nn) and unstable (u). The number of 10-min samples is also indicated. α is the power-law shear exponent; T I amb is the turbulence intensity defined as the standard deviation of horizontal wind speed divided by the mean wind speed. The wind speed and turbulence parameters are obtained from sonic observations at 32 m height.  Considering the small elevation angles and the typical low values of w, we assume w = 0 (Doubrawa et al., , 2020Debnath et al., 2019). Following the approach of Doubrawa et al. (2020), we can combine the u-and v-velocity components into a total horizontal wind vector, U , and Eq. (6) becomes whereθ 0 is the yaw offset and the overbar indicates a smoothed signal, as we apply a moving average operator with a 15-s 275 window to the yaw misalignment to account for any temporal delay from the spatial distances among the mast, turbine's nacelle and SpinnerLidar measurements (Conti et al., 2020a). With Eq. (7), we can reconstruct horizontal wind velocity measures at each individual scanned point within the rosette pattern. Further, we linearly interpolate the reconstructed wind speeds across the rosette pattern into a two-dimensional regular grid with a 2-m resolution, which is sufficient to characterize the spatial characteristics of the wind field in wakes (Fuertes et al., 2018;Conti et al., 2020a).
280 Table 3. Similar as Table 1, but for Strategy III

Lidar-estimated wake deficit
To perform comparisons with predicted velocity deficits from the DWM model, we aim at isolating the contribution of the wake deficit from that of the vertical wind shear in lidar measurements. As defined in Trujillo et al. (2011), the quasi-instantaneous wake deficit profile can be obtained by subtracting the mean vertical shear profile (U amb (z)) from the quasi-instantaneous wake recording as: where U (x, y, z) is estimated from lidar measurements using Eq. (7), and U amb (z) is the relative 10-min average inflow vertical wind speed profile measured at the mast. The deficit is then normalized with respect to the ambient wind speed profile. The v los measurements, and indeed the reconstructed U wind velocities, are defined on a coordinate system that is either attached to the nacelle (nacelle frame of reference, NFoR), which rotates with the yawing of the turbine, or to the ground (fixed frame of 290 reference, FFoR). To perform direct comparisons with the DWM model predictions, the lidar-estimated deficits obtained from Eq. (8) need to be computed in the MFoR. Here, this is performed by tracking the wake center position through the method of Trujillo et al. (2011), where a bivariate Gaussian shape is fitted to the velocity deficit flow field and the wake center is the geometric centroid of the Gaussian:

295
where (µ y , µ z ) define the wake center location, (σ wy , σ wz ) are width parameters of the wake profile in the y and z directions, respectively, (y i , z i ) denote the spatial locations of the lidar measurements, and A is a scaling parameter. Each scanned point of the quasi-instantaneous wake recording can be translated into the MFoR using the estimated µ y and µ z from Eq. (9) (Reinwardt et al., 2020). Therefore, we can compute the multiple wake recordings within a 10-min period in the MFoR, and subsequently compute flow statistics such as the ensemble-average deficit profile as well as the spatial distribution of the wake turbulence 300 in the MFoR. To ensure a high-quality fit, we reject scans where the estimated wake center location is within ≈ 10% of the https://doi.org/10.5194/wes-2020-135 Preprint. Discussion started: 4 February 2021 c Author(s) 2021. CC BY 4.0 License.
lateral bounds of the scanning area and at more than 0.75 D from the hub height in the vertical direction (Conti et al., 2020a;Doubrawa et al., 2020). Figure 4 illustrates ensemble-average measured deficit profiles in the MFoR at 2, 3, 4 and 5 D behind the rotor obtained from all 10-min periods characterized by an incoming wind speed of 7 m/s, and under varying stability regimes (see Table 1 for 305 reference). We can clearly observe the impact of the atmospheric stability and in particular of the associated turbulence levels on the wake recovery behind the rotor. A strong and well-defined symmetric wake deficit shape is seen under stable conditions (top row), whereas the deficits recover faster moving downstream as the atmosphere becomes more unstable (bottom row).

Lidar-estimated wake turbulence
Turbulence measures derived from lidar radial velocity measurements are 'filtered' because of their relatively large probe 310 volume (Peña et al., 2017) and so they are generally lower than those obtained from sonic observations. Nevertheless, if the Doppler spectrum of the v los is available, we can potentially circumvent the averaging effects and estimate the unfiltered variance of v los (Peña et al., 2017;Mann et al., 2010). Mann et al. (2010) assumes that the ensemble-averaged Doppler spectrum over a time period S(v los ) is related to the probability distribution of the v los at the focus distance, and can be computed as: where ϕ(s) is the spatial averaging function of the lidar that depends on the position along the beam s, and p(v los |s) denotes the PDF of v los at the location s. If we assume that the PDF of v los is independent of s, (i.e., there is no velocity gradient along the beam), then Eq. (10) reduces to S(v los ) = p(v los ). As a result, the v los statistics (i.e., mean and variance) can be computed from the first and second central moments of p(v los ) as: where µ v los and σ 2 v los denote the mean and unfiltered variance of v los , respectively. Following the procedure of Peña et al. (2019), we compute the ensemble-averaged Doppler spectrum within 10-min periods by thresholding the noise-flattened spectra with a value of 1.2 and correcting them by subtracting the background spectrum. We accumulate the LOS Doppler spectra onto the regular grid of the scanned area and estimate µ v los and σ 2 v los for each grid cell using Eq. (11). As discussed in Herges 325 and Keyantuo (2019), invalid measurements occur due to the boresight and ground return, as well as the return from the rotating rotor of WTGa2, if in operation. These invalid observations appear as very high return signal in the Doppler spectrum in proximity of low wind speeds (i.e., at approximately 1 m/s) and are removed. The filtering effects due to the probe volume can be quantified by computing the ratio between filtered and unfiltered LOS variances across the rosette pattern; we find ratios in the range 0.8-0.9 at 2.5 D, which vary according to stability conditions (not shown).

330
Examples of 10-min ensemble-averaged Doppler spectra obtained at three fixed locations across the scanned area, a wake center, a wake edge, and a wake-free position, are shown in Fig. 5 for an incoming wind speed of 7 m/s and low turbulence.
A narrow spectrum with a single-peak distribution centered at about 7 m/s for the wake-free location (green) is seen, whereas spectrum broadening effects induced by small-scale generated turbulence are noticeable for the positions within the wake. The wake center (red) shows a wider spectrum with a peak at a significantly lower wind speed than the incoming flow, whereas the 335 wake edge (cyan) shows a double-peak distribution that may be partially due to the inhomogeneity of the wind field along the beam (Herges and Keyantuo, 2019) and also due to the meandering occurring within the analyzed 10-min period. To characterize the spatial distribution of the wake turbulence within the scanned area, we derive σ 2 U estimates directly by applying the variance operator to Eq. (7):

340
where σ 2 U is the variance of the horizontal wind speed, and as shown, covariance terms are neglected. As the LOS is almost never aligned with the u-velocity component across the rosette, except at the center of the pattern, σ 2 v los can be 'contaminated' by the variances and covariances of the other velocity components (Peña et al., 2017). Therefore, the relation in Eq. (12) can lead to inaccurate estimations of the longitudinal velocity variances. Peña et al. (2019) estimated the contamination of different components on the LOS variances for the SpinnerLidar, and showed that the ratio of the unfiltered LOS velocity variance to 345 the variance of the longitudinal velocity component is generally lower than one across the scanned area, except at the center where the ratio is one, and within an area above the center where it can be higher than unity. Although the adopted retrieval assumption in Eq. (12) introduces uncertainties in the turbulence measures, we can account for the expected errors in the Bayesian inference framework. Figure 6 illustrates the spatial distribution of the unfiltered σ 2 U computed in the MFoR, normalized with the u-velocity 350 variance of the ambient wind field measured at the 32-m sonic (σ 2 u,amb ). Under stable conditions, and for a downstream distance of 2.5 D, we can observe an enhancement in turbulence levels in proximity of the rotor tips, especially in the upper part of the rotor (see Fig. 6 (a)). The observed added turbulence is caused by the breakdown of the rotor tip vortices. These features are no longer noticed as the atmosphere becomes more unstable, where a more uniform and less prominent distribution of the turbulence is found (see Fig. 6 (b) and (c)).  Figure 6. Two-dimensional spatial distribution of the horizontal wind velocity variance σ 2 U computed in the MFoR and normalized with the u-velocity variance of the ambient wind field (σ 2 u,amb ) for an incoming wind speed between 7 and 8 m/s under different stability conditions.

Calibration of the DWM model
The calibration of the wake deficit and wake-added turbulence components are conducted in the MFoR using a Bayesian inference framework. We describe the Bayesian model in Sect. 5.1, provide calibration results for the wake deficit in Sect. 5.2, and for the wake-added turbulence in Sect. 5.3. We investigate wake meandering dynamics separately in Sect. 5.4.

Bayesian inference formulation 360
The basis of the Bayesian inference is to estimate the probability distribution of the model parameters based on available observations. Let θ m = {k 1 , k 2 , ..., k m1 , k m2 } be a set of model parameters to be estimated using lidar-derived wake features (i.e., wake deficit and wake-added turbulence profiles in the MFoR) denoted by y d = {y d1 , y d2 , ..., y dn }, where n is the number of available observations. We consider that the experimental data and the model predictions satisfy the prediction error equation: 365 whereĝ(θ m , X m ) denotes the DWM model predictions obtained from a particular set of model parameters (θ m ) and a set of observable variables (X m ). Here X m = {T I amb , U amb , α, C T , x, y, z} includes the inflow wind conditions measured at the mast (T I amb , U amb , α), the rotor thrust coefficient of the turbine (C T ), which is derived from the BEM model (Madsen et al., 2010), and the spatial locations of the scanning pattern (x, y, z). = y + m denotes a random prediction error composed of two terms: the measurement error y and the model prediction error m . The former is described by a zero mean normal 370 distribution with standard deviation σ y , which is determined from field observations. The latter is assumed to have zero mean, which implies unbiased model predictions, and a standard deviation σ m to be determined by the Bayesian estimation along with the model parameters. To facilitate statistical inference, we assume X m as deterministic inputs (i.e., free of uncertainty), and that the model error m is independent on the set of input variables X m , and described by a normal distribution. This implies that the model predictions are normally distributed for a given X m , which is a reasonable choice for wake deficit 375 profiles. The Bayesian approach for model calibration deals with updating the combined parameter set (θ m , σ m ), given a set of observations (y d , X m ) by applying the Bayes theorem: where defined as: The outcome of the calibration is a joint probability distribution of the inferred model parameters. From this joint PDF, we can estimate the posterior PDF of any wake feature simulated by the DWM model, i.e., the wake deficit/wake-added turbulence profiles in the MFoR, or the fully-resolved wakes in the FFoR, among others, which we denote by q: The posterior distribution of the wake feature q in Eq. (16) can be solved numerically using sampling methods (e.g., Monte Carlo simulations) so its first and second moment can be estimated.

Wake deficit parameter estimation
We use lidar-derived wake deficit profiles in the MFoR collected during Strategy I and Strategy II, and employ the Bayesian 400 model to infer uncertainty in the k 1 and k 2 parameters of the eddy viscosity term in Eq. (1). These parameters were found to be the most sensitive on the resulting wake deficit predictions (Keck et al., 2012). The prediction model of Eq. (13) is constructed as following. The experimental data (y d ) comprises two-dimensional ensemble-average lidar-estimated deficit profiles binned according to downstream distances (2, 3, 4 and 5 D), atmospheric stability (i.e., stable, near-neutral and unstable), and wind speed bins of 1 m/s in the range 3-9 m/s. The set of observable variables comprises X m = {T I amb , U amb , α, C T , x, y, z},

405
where the inflow parameters (T I amb , U amb , α) are provided in Tables 1 and 2; the rotor thrust coefficient C T is derived from the BEM model implemented in the aeroelastic code HAWC2 (Larsen and Hansen, 2007) and based on the aerodynamics and airfoil inputs of the SWiFT turbine (Doubrawa et al., 2020) (C T =0.84 that is nearly constant for wind speeds below 9 m/s); and x, y and z refer to the spatial coordinates of the deficits resolved in the MFoR. The uncertainties in measured deficit profiles are computed as y d (r) = σ(r)/ √ n, where σ(r) is the standard deviation of all 10-min deficits within the analyzed case at the 410 radial position r, and n is the number of 10-min periods (also referred to as samples in Table 1).
We select uniform prior distributions on model parameters k 1,prior ∼ U(0.001, 0.2) and k 2,prior ∼ U(0.001, 0.2) on intervals that consider physical constraints, ensuring convergence of results, and covering previous calibrations reported in the literature.
Thus, we employ a MCMC algorithm to sample from the posterior PDFs of the calibration parameters using the Bayesian framework. The inferred joint and marginal posterior distributions of the model parameters (k 1 and k 2 ) are shown in Fig.   415 7, together with point values from earlier studies (Madsen et al., 2010;Keck et al., 2012;Larsen et al., 2013;IEC, 2019;Reinwardt et al., 2020). As shown, the lidar-based wake deficits are informative and we obtain well-defined posteriors that follow a normal distribution with k 1 ∼ N (0.081, 0.017) and k 2 ∼ N (0.015, 0.003). The negative correlation between k 1 and k 2 seen in Fig. 7 indicates the interdependence of the physical induced-effects, as both parameters contribute to turbulence diffusion. It is found that the posterior means of the informed parameters k 1 and k 2 differ from those recommended in the IEC 420 standard (IEC, 2019). Generally, low k 1 and k 2 values attenuate the degree of turbulence mixing in the wake, which lead to strong velocity deficits persisting at farther downstream distances. We provide the statistical properties (mean, variance and coefficient of variation) of the inferred parameters in Table 4 and correlation measures in Table 5.

Wake deficit predictions in the MFoR
We propagate the uncertainties of k 1 and k 2 to predict wake deficits in the MFoR, and compare them with the ensemble-425 average lidar-derived profiles in Fig. 8. First, we observe that the lidar-estimated deficits exhibit a faster wake recovery as the atmosphere becomes more unstable compared to stable regimes. This effect is mainly caused by the enhanced turbulence mixing occurring under unstable conditions, as they are characterized by ambient turbulence levels 2-3 times higher than those of the stable cases (see Table 1). The lidar-observed maximum deficit varies between 30% and 60% within the first five rotor diameters, depending on the inflow turbulence conditions. Similar behaviors were reported in recent lidar measurement deficits are approximately Gaussian under stable to near-neutral conditions, whereas the Gaussian shape is lost under more unstable conditions. This may result from errors in the wake tracking procedure due to the larger meandering amplitudes, and also due to the presence of large-scale turbulence structures in the inflow (Conti et al., 2020a). The rotor thrust is another factor governing the variability of the wake recovery (Zhan et al., 2020b), however, its influence is secondary for the here-analyzed The DWM-predicted deficit profiles with parameters specified by their posterior distributions are in good agreement with the lidar observations for distances beyond 4D (see Fig. 8). For these distances, the turbulence mixing effects dictated by the ambient and self-generated wake turbulence on the deficit recovery are fairly well-captured by the inferred parameters. The nominal model predictions generally fit the observations, whereas an overlap between the measurements and the region of 440 modeling uncertainty is found. The largest deviations between predicted and measured deficits are found at shorter distances (2-3D) and mostly under unstable conditions. These deviations are mainly due to both the model inadequacy to simultaneously fit all the experimental measurements and experimental uncertainties. Similar findings were reported in Keck et al. (2015) and Machefaux et al. (2016), who showed the inaccuracy of deficit predictions at short distances.
It can be observed that the uncertainties in k 1 and k 2 parameters primarily influence the depth of the wake (i.e.,the maximum 445 deficit), while their sensitivity decrease significantly with the outer radial distance. It is also noticed that the uncertainty of the deficit predictions increases for high ambient turbulence and far downstream distances. This is because both k 1 is proportional to T I amb in Eq. (1) and the increasing sensitivity of the model parameters to wake recovering. To provide a measure of the uncertainty of wake deficit predictions, we compute the coefficient of variation COV = σ/µ, where σ is the standard deviation and µ is the mean value of the maximum deficit, obtained by propagating the PDFs of k 1 and k 2 , as in Eq. (16). We find 450 COV = 3% under stable conditions (T amb =0.07) that increases to 6% under unstable conditions (T amb =0.14) for an incoming wind speed of 7 m/s. This result confirms that the uncertainty of wake deficit predictions increases for higher turbulence, but it also shows that uncertainties in k 1 and k 2 parameters do not lead to uncertainty of the same magnitude on deficit predictions resolved in the MFoR (for reference COV k1 = 21% and COV k2 = 19% as reported in Table 4).
We provide comparisons between measured and predicted wake deficit profiles using the calibration from this work as well 455 as those reported in early studies in Fig. 9. For this particular analysis, we analyze predictions at 5D behind the rotor for an inflow wind speed of 7 m/s under stable, near-neutral and unstable regimes. The main discrepancy among the models is the relative sensitivity of the wake recovery to the ambient turbulence. This is primarily governed by k 1 ; however, in the eddy   Table 1 for details on the inflow wind conditions). while the shaded areas indicate its 95% confidence interval.

Improved wake-added turbulence formulation
The wake-added turbulence model (Eq. 3) assumes that turbulent structures (i.e., tip and root vortices) are unaffected by atmo-470 spheric turbulence. The rotor-induced vortices are rapidly disrupted under high turbulence conditions causing the breakdown within the first 2D (Madsen et al., 2005). However, vortices can persist and extend at farther distances under low to moderate turbulence combined with stable stratification conditions (Ivanell et al., 2009;Subramanian et al., 2018;Conti et al., 2020a).
Thus, the wake-added turbulence profiles can exhibit a more pronounced double-peak feature in the proximity of the rotor tips or a more uniform distribution depending on the atmospheric turbulence conditions (this effect is also seen in Fig. 6). Eq. (3) 475 also assumes radially symmetric wake-added turbulence profiles. However, the inflow vertical wind shear re-distributes the wake turbulence. Enhanced turbulence levels are actually observed in the proximity of the upper tip of the rotor blade (Vermeer et al., 2003;Chamorro and Porte-Agel, 2009;Conti et al., 2020a). Figure 6 (a) shows this effect.
Due to these two assumptions, we propose an improved semi-empirical formulation of the wake-added turbulence scaling factor k * mt , which produces wake profiles in better agreement with the lidar observations. This is achieved by relating both the 480 depth and the velocity gradient terms to the ambient turbulence, and by including the effect of the inflow vertical wind shear on the vertical velocity deficit gradient as: where U def (y, z) * = (U (z)U def (y, z))/max (U (z)U def (y, z)), and k * m1 , k * q1 , k * m2 , k * q2 are parameters to be determined using Bayesian inference. Figure 10 illustrates the two-dimensional profiles of the depth and velocity deficit gradient terms of Eqs.

485
(3) and (17). As illustrated, by including the term U def (y, z) * , we obtain a turbulence field that mimics qualitatively well the observed enhanced turbulence within the upper wake region. 1 U def (r)

Estimation of wake-added turbulence parameters
The calibration parameters in Eq. (17) are inferred based on the lidar-estimated wake-added turbulence profiles in the MFoR collected during Strategy II. For this particular dataset, the SpinnerLidar scans at a fixed distance of 2.5D ensuring about 298 490 scans for every 10-min period. The inflow characteristics are reported in Table 2. As the Doppler LOS velocity spectrum is available, we derive unfiltered LOS variances as described in Sect. 4.2, and subsequently the turbulence intensity as the ratio of the standard deviation to the mean of the horizontal wind speed. From Eq.
(2), we can isolate T I add (wake-added turbulence term) by firstly resolving the wake recordings in the MFoR, which eliminates the contribution of T I m , and then by subtracting T I amb that is derived from the 32-m sonic observations in front of the rotor.

495
As a first step, we fit the 'original' analytical formulation of k mt from Eq. (3) to each individual lidar-estimated wake-added turbulence profile to estimate the values of k m1 and k m2 using a simple least-squares optimization algorithm. Figure 11 a,d show the relation between the estimated optimal parameters (in markers) and the ambient turbulence. It is shown that k m1 increases and k m2 decreases almost linearly for increasing turbulence intensity. This indicates the strong effect of the deficit gradient term (proportional to k m2 in Eq. (3)) under low turbulence conditions, which both amplifies the double-peak feature 500 of the wake turbulence profile at the rotor tips and mimics the wake vortices. As the ambient turbulence increases, k m1 and lower k m2 become larger, which indicates a more uniform distribution of the wake turbulence. These effects are observed in Fig. 11 b,c (see markers), which shows the spanwise distribution of the lidar-derived wake-added turbulence for two 10-min periods with relatively low and high values of ambient turbulence intensity, 7% and 12%, respectively. We also compare the measured and predicted vertical distribution of T I add in Fig. 11 e,f, which shows that slightly improved predictions can be 505 obtained using Eq. (17), i.e., considering the effects of the atmospheric shear on wake turbulence.
The resulting posterior PDFs of the parameters follow a normal distribution (not shown) with statistical properties tabulated in Table 4. We provide the mean values of the inferred k * m1 , k * q1 , k * m2 , k * q2 values in the form of a linear regression model that relates the depth and deficit gradient terms of Eq. (17) to the ambient turbulence in Fig. 11 a,d. The shaded area represents the 515 95% confidence interval, which is obtained by propagating the PDFs of the turbulence-related and wake deficit parameters. As shown, the predictions are characterized by a relatively high degree of uncertainty. This is because, e.g. there are few available observations to characterize turbulence, the measurements uncertainties are relatively high, and the modeling simplifications, such as the semi-empirical k * mt factor. We provide the correlations among all the inferred parameters obtained from the joint PDF in Table 5. As expected, the wake deficit parameters k 1 and k 2 are negatively correlated to k * m1 and k * q1 , as k * m1 and k * q1 520 are proportional to the wake deficit term in Eq. (17), which is the main driver to the intensity of the wake-added turbulence.

Wake meandering
Here, we investigate the relationship between the inflow turbulence fluctuations and the lidar-tracked wake positions to characterize the large-scale eddies responsible for the meandering. The analysis is carried out by comparing the spectra of the lidar-tracked meandering time series, which are derived by means of the tracking algorithm in Eq. (9), with that simulated from 525 the meandering model in Eq. (4), where v c and w c are obtained from the sonic observations at 32 m. Yaw misalignment from SCADA is also included. We conduct the analysis using the dataset collected at 2.5 D behind the rotor during Strategy II, and classify all the available 10-min periods according to stability to derive ensemble-average spectra from multiple observations with similar inflow conditions. Note that the measured inflow wind speeds are below rated, and turbulence levels from all 10 min periods within each stability class are similar as those reported in Table 2.

530
Results of the spectral analysis for both lateral and vertical meandering are shown in Fig. 12. Here, the ensemble-average spectra from the SpinnerLidar observations are compared to those from the meandering model without low-pass filtering the incoming turbulence fluctuations (denoted as DWM wf ). The spectra are normalized with their relative variances and plotted as function of the commonly used Strouhal number St = f D/U ∞ , where f denotes frequency and U ∞ is the aggregated wind speed from the ensemble-average statistics. As shown, the slope of the lidar-based spectra (red lines) matches that of DWM wf 535 (blue lines) up to St = 0.3-0.5, which corresponds to three and two rotor diameters, respectively. For St > 0.5, the energy content of the lidar-estimated spectra remarkably decreases compared to that of DWM wf . These observations indicate that large-scale turbulent structures (>2D) are dominant in the wake meandering (Trujillo et al., 2011;Heisel et al., 2018). When compared to the stable case, the spectra under unstable conditions show a higher energy content at large turbulence scales, or equivalently low Strouhal number. lateral wake-added turbulence profiles resolved in the MFoR at hub height and obtained at 2.5D behind the rotor are shown in (b) and (c), for TI amb = 7% and 12%, respectively, whereas the relative vertical profiles are shown in (e) and (f). The error bars indicate the measurements uncertainty, whereas the shaded areas that of the model predictions relative to the 95% confidence interval. Figure 12 shows that a stochastic description of the large-scale eddy size might be appropriate. Thus, we describe the largescale eddies responsible for the wake meandering by introducing the stochastic variable D m , which is normally distributed with mean equal to 2.5D (it corresponds to the wake diameter at 5D behind the rotor) and a standard deviation of 0.3D based on the observations in Fig. 12. The resulting 95% confidence intervals are shown as shaded areas in Fig. 12. The uncertainty in D m is found negligible when computing wake meandering time series; this is shown in Appendix A.

6 Validation of the DWM model
The validation is performed by resolving wake fields in the FFoR, thus the simulated wakes include the combined effects of the velocity deficit, added turbulence, and wake meandering dynamics in both lateral and vertical directions. This analysis is carried out using data from Strategy III, i.e., at a fixed distance of 5D behind the rotor, ensuring a sufficient amount of scans to derive unfiltered turbulence estimates as well as wake meandering time series within a 10-min period. The analyzed dataset The 95% confidence interval in the large-scale eddies definition is shown by the grey area (see text for more details).
is primarily characterized by stable conditions as seen in Table 3, i.e., low turbulence intensities (6-8%) and strong vertical shears (α = 0.18-0.38). The wind speed is mainly lower than 9 m/s, thus WTGa1 operates below rated power.

Correction for rotor induction effects
The SpinnerLidar measurements collected during Strategy III were taken in the induction zone of the WTGa2. For induction correction, we employ the two-dimensional induction model of Troldborg and Meyer Forsting (2017), which accounts for both 555 longitudinal and radial variation of the induced wind velocity: where a 0 is the induction factor at the rotor center area, a 0 = 0.5(1 − √ 1 − γ a C T ), γ a = 1.1 (Troldborg and Meyer Forsting, 2017), ξ x = x/R is the distance upfront the rotor normalized by the rotor radius, ρ a = y 2 + z 2 /R denotes the radial distance from the rotor center axis, and a = ρ a / λ a (η a + ξ 2 x ) being β a = √ 2, α a = 8/9, λ a = 0.587, η a = 1.32 (Troldborg and factor in Eq. (18). The estimated induction factors indicate that the wind speed can be reduced at hub height by up to 12% upstream the WTGa2 and below rated power (not shown).

Uncertainty propagation of simulated wake fields in the FFoR
We derive the two-dimensional spatial distribution of the mean wind speed in the wake region (U F F oR ) by the convolution 565 between the wake deficit in the MFoR (U def,M F oR ), and the PDF of the meandering path (f m ) (Keck et al., 2015): where (y m , z m ) denote the spatial coordinates of the wake meandering time series, and (y m, , z m, ) are measures of their relative uncertainties. We introduce these errors to account for incorrect wake tracking positions that can arise due to the adopted wake tracking algorithm. y m, and z m, are assumed to be uncorrelated and to follow a normal distribution with zero mean 570 and standard deviation such as the 95% percentile corresponds to approximately 4 m, which is twice the resolution adopted to interpolate SpinnerLidar measurements onto the regular grid (see Sect. 4). Note that the atmospheric shear profile U amb (z) is superposed after the wake deficit calculation (Madsen et al., 2010). Similarly, the two-dimensional spatial distribution of the u-velocity variance (σ 2 u F F oR ) can be computed as: (19), and σ 2 u amb is the variance of the ambient u−velocity component. Alternatively, U F F oR and σ 2 u F F oR can be equivalently derived by superposing the wake deficit and the wake-added turbulence factor k * mt on stochastic turbulence fields with constrained meandering path, and subsequently by computing the first-and second-order statistics of the synthetic wind fields. However, the analytical forms in Eqs. (19) and (20) can be easily used to propagate the posterior PDFs of the calibration parameters to predict profiles of U F F oR and σ u F F oR with relative 580 uncertainties. The inflow parameters (U amb , α and σ 2 u amb ) in Eqs. (19) and (20) are required inputs for the DWM model and are estimated from mast measurements.
We compute U F F oR and σ u F F oR from Eqs. (19) and (20) by constraining the meandering path (f m ) either using the lidartracked meandering time series, which we denote as DWM * , or by using the meandering model in Eq. (4) with low-pass filtered v c -and w c -velocity fluctuations obtained from the 32-m sonic and the yaw misalignment of WTGa1 obtained from SCADA, 585 which we denote as DWM * * . In the latter model, the low-pass filtered frequency is defined as function of the stochastic Table 4. Mean (µ), standard deviation (σ) and coefficient of variation (COV= σ/µ) estimated from the posterior PDFs of model parameters.
The values of k1, k2, σ def , km1, kq1, km2, kq2 and σ add are determined using Bayesian inference. Dm denotes the spatial size of the large-scale eddies governing wake meandering dynamics, and ym, and zm, denote the wake tracking position errors expressed in meters.

Modules
wake deficit wake-added turbulence wake meandering  variable associated to the size of the large-scale eddies (D m ) as discussed in Sect. 5.4. Figure 13 shows the comparison between observed and predicted (with both DWM * and DWM * * models) profiles of the mean wind speed and its standard deviation obtained from two different 10-min periods.
A good agreement between measurements and predictions is found. The vertical profile of the mean wind speed exhibits a 590 single-peak shape resulting from the combined effects of the inflow vertical shear (modeled by a power-law) and the wakeinduced Gaussian-like deficit shape. The wake turbulence in the lateral direction exhibits a double-peak shape with larger values near the locations associated with strong velocity gradients that are further enhanced by the wake meandering. Enhanced turbulence levels in proximity of the upper wake region are observed from lidar measurements. The deviations between DWM * -and DWM * * -model predictions are more pronounced for the simulated turbulence than for the wind speed fields, and 595 are exclusively due to differences in the meandering representations.
The wake simulations uncertainties shown in Fig. 13 are determined by propagating uncertainties in model parameters (Table   4), and by accounting for relative correlations (Table 5) using Monte Carlo simulations. The uncertainties of lidar-measured wind speeds account for volume averaging effects and for errors introduced by the retrieval assumptions (e.g., w sin(φ) = 0 m/s) . The uncertainties of lidar-derived turbulence (here defined as the standard deviation of the wind 600 speed) account for errors introduced by neglecting cross-contamination effects in Eq. (12) . To evaluate the performance of the DWM model, we calculate two flow metrics that are relevant in aeroelastic simulations Murcia Leon et al., 2018) and compare them to relative measured quantities: the rotor effective wind speed (U ef f ) defined as the weighted sum of wind speeds across the rotor area, and the effective wake turbulence (σ u,ef f ) that is derived as the weighted sum of turbulence estimates across the rotor. Figure 14 shows the one-to-one comparison between 605 measured and DWM * -model predicted flow metrics for all 10-min periods. Similar statistics are obtained with the DWM * *model (not shown). We find a slope that deviates < 1% from unity and R 2 = 0.95 for U ef f and a bias of 4% and R 2 = 0.93 for σ u,ef f . The observed scatter is explained by the large measurement uncertainties, by the uncertainties of inflow wind parameters, and by those of the model predictions. The latter is estimated by propagating uncertainties of model parameters provided in Table 4, which result in a COV of 1% for U ef f , and of 3% for σ u,ef f . These findings indicate that the inferred 610 uncertainties in model parameters do not propagate into an error of the same magnitude on the fully-resolved wakes that are eventually input to aeroelastic simulations.
Further, a sensitivity analysis indicates that the variations in U ef f and σ u,ef f are mostly explained by the uncertainties in k 1 and k 2 and the wake tracking position y m, , z m, . This is shown in Appendix B. The contribution of the wake-added turbulence from Eq. (17) on the total wake turbulence σ u,ef f in the FFoR is marginal and accounts for approximately 2-7%. 615 The turbulence induced by the meandering of the wake deficit is thus the major source of added turbulence (Madsen et al., 2010

Discussion
Atmospheric stability significantly impacts wake evolution and recovery (Iungo et al., 2013;Zhan et al., 2020b). SpinnerLidar measurements of the wake show that the wake recovers faster under unstable compared to stable conditions primarily due to the 620 high turbulence levels of the former. The wake effects can be predicted accurately by the DWM model when using appropriate calibration parameters. Further, accurate reproduction of the wake meandering dynamics in both lateral and vertical directions is key for accurate wake simulations.
We recommend conducting DWM model calibrations by resolving the wake flow features such as the wake deficit and wakeadded turbulence profiles in the MFoR. This requires high spatial and temporal resolution scanning strategies, which cover the 625 two-dimensional plane upfront the rotor to track both the horizontal and vertical wake displacements, and to reconstruct the spatial distribution of the wind and turbulence fields. When using power production for such calibrations (the IEC standard (IEC, 2019) values are based on this approach), we are unable to distinguish between uncertainties from inaccurate wake deficit predictions and those from erroneous wake meandering representations. The latter plays an important role in the accuracy of the fully-resolved wake fields, which are inputs to aeroelastic simulations.
Datasets that include observations of the wake deficit profiles under varying stability conditions, inflow wind speeds, and downstream distances are imperative. Here, we find that the velocity deficit's recovery rate for increasing turbulence (see Fig.   9) is the main difference among calibrations.
As wind turbines are typically spaced 5D and beyond, it is recommended to focus future measurement campaigns to include those regions. To this extent, Bayesian inference is a valuable approach for updating the PDFs of model parameters estimated within this work by directly including data from future observations while retaining information from the earlier observations. Power and load validation analyses using the proposed calibration parameters at multiple sites can further increase the confidence in our calibration methodology.
Characterizing wake turbulence using lidars is challenging due to the limited sampling frequency and probe volume effects (Peña et al., 2017). We demonstrate the usefulness of Doppler radial velocity spectra to compute unfiltered LOS variances in 640 wake conditions. In addition to the enhanced turbulence intensity, the wake turbulence is characterized both by being highly isotropic and reduced turbulence length scale compared to the ambient turbulence (Madsen et al., 2005). These turbulence characteristics were not investigated in this work.

Conclusions
We analyzed high spatial and temporal resolution SpinnerLidar measurements of the wake field collected at the SWiFT facility 645 and derived wake features such as the wake deficit, wake-added turbulence, and wake meandering under varying atmospheric stability conditions, inflow wind speeds, and downstream distances. The SpinnerLidar-estimated wake characteristics computed in the MFoR were used to determine uncertainties in the DWM model parameters using Bayesian inference. The uncertainties in model parameters were propagated to predict fully-resolved wake flow fields in the FFoR. This approach allowed us quantifying uncertainties in the DWM-simulated wake fields and to investigate the sensitivity of the DWM model parameters on flow 650 features that primarily affect power and load predictions.
The SpinnerLidar-derived wake deficit profiles revealed the strong impact of atmospheric stability on wake evolution. In particular, we observed the faster recovery of the deficit under unstable compared to stable regimes, as higher turbulence intensities characterized the former. These effects were accurately reproduced by the eddy viscosity term of the DWM model with the inferred parameters for distances beyond 4D. Our results indicate that the currently adopted parameters in the IEC 655 standard lead to conservative velocity deficit predictions (up to 18% for moderate to high ambient turbulence T I amb ≥ 12%) at distances up to 5D behind the rotor.
We proposed and verified an improved semi-empirical formulation of the wake-added turbulence model that captured the effects of the atmospheric shear and the ambient turbulence on the spatial re-distribution of the wake turbulence observed at 2.5 D. We also demonstrated that the wake meandering is the major source of added turbulence in the wake region.

660
The underlying hypothesis of the DWM model, i.e., wakes are advected passively by the large-eddies in the incoming wind field, was verified by means of the SpinnerLidar-tracked meandering time series. The spectral analysis indicated that large-eddies associated with sizes larger than 2D are the main responsible for the wake meandering; however, the large-eddies 'definition' had only marginal effects on the predicted wake fields. Accurate tracking of the wake center position was the most influential factor in simulating wake flow fields accurately. We expect that it also plays a central role in the accuracy of power 665 and load predictions.
In the Part II of this work, we will quantify uncertainties in power and load predictions based on the proposed calibration at two different sites, the SWiFT facility and the Nørrekaer Enge wind farm in Denmark (Peña et al., 2017;Dimitrov, 2019;Conti et al., 2020b).
Appendix A: Comparisons of wake meandering time series