Assessing Obukhov Length and Friction Velocity from Floating Lidar Observations: A Data Screening and Sensitivity Computation Approach

: This work presents a parametric-solver algorithm for estimating atmospheric stability and friction velocity from ﬂoating Doppler wind lidar (FDWL) observations close to the mast of IJmuiden in the North Sea. The focus of the study was two-fold: (i) to examine the sensitivity of the computational algorithm to the retrieved variables and derived stability classes (the latter through confusion-matrix theory), and (ii) to present data screening procedures for FDWLs and ﬁxed reference instrumentation. The performance of the stability estimation algorithm was assessed with reference to wind speed and temperature observations from the mast. A ﬁxed-to-mast Doppler wind lidar (DWL) was also available, which provides a reference for wind-speed observations free from sea-motion perturbations. When comparing FDWL- and mast-derived mean wind speeds, the obtained determination coefﬁcient was as high as that of the ﬁxed-to-mast DWL against the mast ( ρ 2 = 0.996) with a root mean square error (RMSE) of 0.25 m/s. From the 82-day measurement campaign at IJmuiden (10,833 10 min records), the parametric algorithm showed that the atmosphere was neutral (31% of the cases), stable (28%), or near-neutral stable (19%) during most of the campaign. These ﬁgures satisfactorily agree with values estimated from the mast measurements (31%, 27%, and 19%, respectively).


Introduction
Over the past decades, wind energy (WE) has achieved an important position in the global energy market due to its technical improvements and relevant advantages in terms of environmental impact [1,2]. In particular, interest has been rising in offshore WE due to the strong homogeneous winds over the ocean [3]. Nowadays, offshore WE is one of the core technologies of the European roadmap to becoming carbon neutral [4]. However, offshore wind farms still have very high production and deployment costs compared to onshore facilities [5].
In order to reach the European goal of 150 GW of installed offshore WE capacity by 2030, efforts toward commercial competitiveness are being dedicated by the WE industry and research institutions [6]. Traditionally, the wind resource of potential wind farm areas In the present work, we constructed a parametric solver algorithm along with a basic data screening procedure that, relying on FDWL measurements only and MOST, enables the estimation of the Obukhov length and friction velocity. Additionally, a sensitivity study was conducted for the retrieved variables and estimated atmospheric stability classes. In contrast to the approaches by Holtslag et al. [37] and Beljaars et al. [35], the parametric algorithm assesses the atmospheric stability without needing temperature inputs. Complementary to the method proposed by Basu [38], the proposed algorithm can be extended to an unlimited number of measurement heights.
The work is structured as follows: Section 2 describes the 3-month measurement campaign carried out at IJmuiden as well as the instrumental setup. Section 3.1 introduces the notation used in this work, Section 3.2 revisits MOST and atmospheric stability classifications. Section 3.3 describes the parametric solver algorithm. Sections 3.4 and 3.5 outline the derivation of the Obukhov length and friction velocity. Section 3.6 describes the data screening methods used to enhance measurement precision and data quality for both FDWL and the mast. Section 4 provides a discussion of the results and, finally, Section 5 provides our concluding remarks.

Materials
The validation campaign of the EOLOS TM test FDWL at the IJmuiden test site (52.848 N, 3.436 E) in the North Sea was run from April to June 2015. The aim of the measurement campaign was to validate the EOLOS TM FDWL precommercial buoy against the reference IJmuiden mast. Figure 1 shows the instrumental set-up of the campaign, depicting the instruments location and measurement heights. The IJmuiden metmast was installed 85 km far from the coast of The Netherlands, where the water depth is about 28 meters. Its top was 92 m above the Lowest Astronomical Tide (LAT), and its structure consisted of three booms pointing 46.5 • , 166.5 • , and 286.5 • , clockwise from cardinal north. Hereafter, all height values will be assumed as height above LAT. Three Thies TM First Class Advanced cup anemometers were installed at 27 m and 58.5 m in height (one on each boom), as well as three Metek TM USA-1 sonic anemometers at 85 m. Wind direction observations were recorded at 26.2 m, 57.7 m and 87 m by Thies TM First Class windvanes (one on each boom). Horizontal wind speed (HWS) and direction at ten different measurement heights (from 90 to 315 m) were measured by a ZephIR TM 300 DWL installed on a platform at 20.88 m in height (hereafter called the mast-DWL). This lidar measured 50 Line of Sights (LoSs) at equally spaced azimuth angles (7.2 • azimuth step between LoSs) at a sampling rate of 50 Hz along a conical scan with elevation of 30 • [39,40]. Two Vaisala TM HMP155D were installed at 21 and 90 m, which provided temperature and humidity observations every 0.25 s. Air pressure measurements were recorded at 21 m and 90 m by two Vaisala TM PTB210 every 0.25 s as well. A TRIAXYS TM wave-and-current buoy located near the metmast provided measurements of average wave height, wave period, current direction, and water temperature. The wave-buoy-derived observations were stored as hourly averaged measurements. Relevant IJmuiden site instruments used here are listed in Table 1. Table 1. Main specifications of the instruments from the IJmuiden test site. Detailed information about the sensors arrangement and specification can be found in [40].

Sensor
Parameter Unit Sampling Rate Height (1)  The EOLOS TM lidar buoy is a pre-commercial device specially designed to host a vertically pointed ZephIR TM 300 lidar. This 3.77-m wide, 3-ton device was conceived to flexibly and reliably assess the wind resources in order to meet wind-energy industry requirements. It hosts multiple sensors to measure the buoy's attitude as well as the main atmospheric parameters. The FDWL was configured to measure the wind vector at 25, 56, and 83 m to match the anemometers' heights of the mast. Additionally, the FDWL recorded measurements at a height of 38 m for calibration purposes.

Wind Notation Conventions
Wind measurements from the instruments above provided the three instantaneous wind vector components (u, v, and w) corresponding to the components in the horizontal x and y, and vertical z directions, respectively. Fluctuations u , v , and w are defined as the difference between the actual instantaneous velocity components and their respective mean velocities, e.g., u = u − U, where U is the mean value of u. Although we often use a Cartesian coordinate system so that the x, y, and z axes point north, east, and south, here, we align x along the wind direction. We skip usage of U to denote the mean wind component in order to avoid confusion with the vector notation − → U , which is used to denote a specific set of HWS measurements in subsequent sections. We retain, however, usage of the over-bar, (.), to denote mean over time (10 min) for all other variables.

Surface-Layer theory
The vertical wind gradient for a neutral, homogeneous, and stationary flow can expressed as [41] ∂U ∂z where u * is the friction velocity and l is the mixing length [42]. In the surface layer, which corresponds to the lowest 5-10% of the boundary layer, the magnitudes of stress and turbulent fluxes generally vary by less than 10% with height [43]. Thus, within the surface layer, two assumptions are commonly made: (i) the friction velocity variation with height is negligible, and (ii) the mixing length increases with height as l(z) = κz [44], where κ ≈ 0.4 is the von Kármán constant. The logarithmic wind profile for neutral atmospheric conditions is the integral of Equation (1): where z 0 is the roughness length, which is the height at which the mean wind speed becomes zero. Over sea, the roughness length is commonly expressed using the Charnock's relation, where g = 9.81 m/s 2 is the gravitational acceleration and α c is the Charnock's parameter, α c ≈ 0.012 [45]. According to MOST, the influence of atmospheric stability over the wind speed gradient within the surface layer is expressed by the dimensionless wind shear, φ m , as The dimensionless wind shear takes different functional forms depending on the atmospheric stability conditions. It is usually described by the experiment-based Businger-Dyer functions ( [46,47]): where L is the Obukhov length, z L is the dimensionless stability parameter, and β = 6.0 and γ = 19.3 are empirical constants [48], later validated by Holtstag et al. [33] for the IJmuiden site.
Thus, under non-neutral conditions the diabatic wind profile can be expressed as where Ψ m z L is the stability correction function [32,43] where x = 1 − γz L 1/4 . Although the Businger-Dyer functions were empirically derived from data collected over land, their applicability has been successfully tested and validated for the estimation of the offshore wind profile [33,49]. The Obukhov length can be used to classify atmospheric stability conditions, classified into stable (L > 0), neutral (L → ±∞), or unstable (L < 0). With a view to Section 4 two atmospheric stability classifications are listed in Tables 2 and 3 according to the definition classes for L of van Wijk et al. [50] and Gryning et al. [51].
The behavior of the stability correction Ψ m z L as a function of the Obukhov length L is depicted in Figure 2.  Table 3. Labels: vu, very unstable; u, unstable; nnu, near-neutral unstable; n, neutral; nns, near-neutral stable; s, stable; and vs, very stable.  [50].

Atmospheric Stability
Obukhov Length Range (m) Finally, we introduce the wind-speed ratio between HWSs at two heights, z 1 and z 2 , z 2 > z 1 , is therefore, The wind-speed ratio is a proxy of wind shear that will be used to analyze the dependence of wind shear on the dimensionless stability, z L , further.

Parametric Wind Model Estimation
We estimate atmospheric stability based on FDWL HWS profiles measured at a discrete number of heights. We propose solving Equations (6) and (7) for model variables L, u * , and z 0 . The algorithm's sensitivity to these three variables and consequent dimensionality reduction of the problem (two-variable problem, L and u * ) is provided in Section 4.2.
The optimization problem can be formulated as where the function U(z, L , u * , z 0 ) is the parametric wind profile model formulated by Equation (6) and piece-wise by Equation (7). Below, functions are assimilated into vectors and treated indistinctly. A constrained nonlinear least squares (NLSQ) method is used to solve the model parameters, L, u * , and z 0 , minimizing the error norm between the model vector, U(z), and the FDWL-measurement vector, U FDW L . The diabatic correction function Ψ m ( z L ) (Equation (7), Figure 2) becomes ill-posed for Obukhov lengths close to zero (L = 0) and exhibits a nearly flat behavior for large values of the Obukhov length (|L| → ∞). Additionally, Ψ m ( z L ) is a piece-wise function. To facilitate NLSQ convergence, we define upper and lower search bounds for each of the parameters to be optimized as The limiting bounds for the roughness length z 0 in offshore environments follow [43,52]. The friction velocity bounds are in agreement with those found in different experimental campaigns in the literature [53,54] and represent roughly a factor 1.5 over the range of practical friction velocities retrieved by the sonic anemometers (Section 3.5, Equation (14)). The search bounds for the Obukhov length, L, are, e.g., a factor two and four over the practical neutral stability thresholds given by Van Wijk et al. [50] (Table 2) and Gryning et al. [51] (Table 3), respectively. Larger values do not tend to substantially improve the algorithm search performance (see Section 4.2).
Optimization starting vector. Two start points, one positive, L + ∈ [1,2000], and one negative, L − ∈ [−2000, −1], are set to L + = 500 m and L − = −500 m, respectively, as initial guesses for the Obukhov length. This allows us to avoid the asymptotic discontinuity of Ψ m ( z L ) at L = 0 ( Figure 2). The initial guesses for the friction velocity, u * , and roughness length, z 0 , are set to 0.7 m/s and 5 × 10 −3 m, respectively, which are the mean values of these variables in the above search ranges. From these two starting vectors, (L + , u * , z 0 ) and (L − , u * , z 0 ), a positive and a negative search is started by the NLSQ algorithm and two candidate solution vectors are obtained. The vector with the smallest error norm is chosen as the solution.

Reference Measurements: Atmospheric Stability
Reference metmast-derived Obukhov length and friction velocity are needed in order to assess the performance of the proposed algorithm in Section 3.3 above. In the present and following subsection (Section 3.5), different methods to estimate such reference parameters are presented.
Conventionally, the Obukhov length is directly computed from turbulence fluxes measured (from wind and temperature observations at similar heights), e.g., with a sonic anemometer [43]. However, the sonic anemometer at 85 m available only recorded velocity observations. Alternatively, the method by Grachev and Fairall [55] to estimate the Obukhov length, L, provides a way out when air and sea temperature observations are available, as is the case. They found empirical dependency between the bulk-Richardson number, Ri, and the dimensionless stability parameter. Their method has been used in recent studies [33,56].
The bulk Richardson number is an approximation of the gradient Richardson number in which actual local gradients in an atmospheric layer are approximately computed by measurements at a couple (or a series of) discrete heights. The bulk Richardson number is defined as [43] where is the difference between the mean virtual potential temperature at the top layer height, z top , and bottom layer height, z bottom and ∆z = z top − z bottom . Here, z top = 21 m above LAT and z bottom = 0 m above LAT. θ v is the mean potential temperature in the layer and U is the mean wind observed at z top (at this point we note that U is formally defined as In practice, U is the 10 min averaged wind observed at 27 m. The potential temperature is estimated from the expression [43,57] θ where T (K) is the temperature, P 0 = 1000 hPa is the reference pressure, P is the air pressure, R ≈ 287 J/(K·kg) is the gas constant of air, C p ≈ 1004 J/(K·kg) is the specific heat capacity at a constant pressure for air, and r is the mixing ratio (unsaturated air).
In this study, we computed the mean potential temperature θ v at 21 m from the mast temperature, pressure, and humidity data at 21 m. θ v at sea-level was computed from the wave-buoy-measured 60 min mean water temperature [58], which was interpolated down to 10 min resolution, according to wind-energy standards [11]. The pressure and relative humidity at sea level were calculated by extrapolating the mast pressure and humidity vertical profiles from the sensors at 21 m and 90 m ( Figure 1) down to 0 m. Pressure and relative humidity at 0 m as well as water-temperature at 10 min resolution were computed via shape-preserving piecewise cubic interpolation [59,60].
Subsequent to the calculation of Ri, the dimensionless stability parameter ζ = z/L is computed as The constants 10 and 5 above are empirical values [55]. The latter implies a critical stable Ri of 0.2 so the theory has a limitation when the conditions are very stable. Besides, when using the bulk Richarson number approximation, the thicker the layer, the more likely to smooth out large gradients in the layer and, therefore, to misestimate the occurrence of turbulence [43].
When we derive L from Equation (12) we refer to it as L Ri , where z is a reference height ensuring the validity of the bulk Richardson model (Equation (10)) above in the layer [z bottom , z top ]. In the standard methodology ( [55]), z is taken as the sonic-anemometer observational height (11 m). However, in the case of IJmuiden, temperature and particularly the wind-speed observational heights are relatively high, 21 and 27 m, respectively. Although a rough estimation of the reference height can simply be taken as the mean height between the air observational height and sea level [33], in this case, z ≈ (z top + z bottom )/2 = 10.5 m, accurate estimation of the reference height is more involved.
In order to estimate the most accurate meaningful reference height, z , the following procedure was carried out: different estimations of L Ri (z ) were retrieved from Equation (13)  Next, the MOST model wind-speed ratio of Equation (8) was computed for all sets of the L Ri (z ) as well as the coefficient of determination, ρ 2 , between the measured and modeled wind-speed ratios. The z (equivalently, L Ri (z )) with the highest ρ 2 yielded the sought-after reference height. Thereby, z = 15.5 m was found to be the reference height ensuring the highest ρ 2 after the outlier screenings applied over L Ri (z ) and the measured wind-speed ratio (see the outlier screening and the results in Section 4.3).

Reference Measurements: Friction Velocity
The friction velocity is a form by which a shear stress may be rewritten in units of velocity such as the velocity that relates wind shear between layers of flow, e.g., for neutral (2)). Two reference methods are considered for comparison: The first solves Equation (6) for u * given the Obukhov length L Ri estimated from the bulk Richardson number using Equation (13) and mast HWS measurements at 27 m. The roughness length z 0 is computed from Equation (3). This method can be viewed as "mast" 1D (u * 1−D ).
The second uses directly the measured turbulent fluxes [43], For these computations, we used the sonic measurements at 85 m.

Data Screening
To begin with, FDWL, mast-sensors and mast-DWL data were averaged to a common uniform time resolution of 10 min, which is according to WE standards [11]. Accordingly, the wave-buoy-derived water temperature was resampled from 60 to 10 min as already mentioned in Section 3.4.
Regarding anemometers' data, the screening procedure relied on selection of the true wind direction (TWD) [40]. To do this, WD and HWS measurements from all the three anemometers at each of the three metmast measurement heights (27, 58.5 and 85 m) were checked. Figure 3 sketches a top-view of mast depicting the positions of the wind vanes and anemometers at each height. First, the middle WD (MWD) was computed as the median of the WD measurements by the three anemometers. Then, the TWD was derived following these criteria (see Figure 3a): • If the MWD lay on the red crown, the TWD was computed as the mean WD between vanes HxxB240 and HxxB120, Finally, the true wind speed was chosen as the one retrieved from the anemometer over which the TWD was within a confidence angle range (the wake-free range) of 90°± 30° (  Figure 3b) in relation its boom orientation.
Concerning the lidar instruments, HWS measurements lower than 2 m/s or higher than 100 m/s were rejected on account of the ZephIR-300 manufacturer's specifications [22]. Regarding the FDWL, it is well-known that motion-induced effects on the retrieved wind vector become prominent at a 1-s temporal resolution but are negligible at the 10 min level. For the latter case, motion effects increase fluctuations on the wind speed [21,39,61]. At the 10 min level, it was shown [62,63] that the spatial variation (SV) parameter (see below) can be used as a threshold to filter data as it represents a trade-off between key performance indicator (KPI) improvements and data availability. FDWL measurement outliers were identified by plotting the HWS differences between the metmast and FDWL measurements as a function of different FDWL internal parameters, namely, the backscatter, bearing (compass), points in fit out of the 50 radial velocities within a full scan (see below), SV, precipitation, and visibility. For each internal parameter, we considered as outliers the samples which yielded differences higher than 2 m/s in relation to metmast reference. The same analysis was conducted with the mast-DWL, which served as the reference lidar.
Backscatter refers to the intensity of the Doppler lidar echo. Bearing refers to the lidar pointing direction with respect to the Earth's magnetic North. Points in fit and SV are parameters related to the Velocity Azimuth Display (VAD) algorithm used by the lidar to retrieve the wind vector under the assumption of uniform wind. The SV is an indicator of the goodness of the VAD fitting of the radial velocities within the lidar scanning circle used to retrieve the wind vector at a given height. Therefore, the SV can be understood as a turbulence parameter describing the variation degree of the radial wind speeds within the scan [64]. High SV values are associated to measurements for which the hypothesis of uniform wind during the lidar scan is no longer valid [39]. This leads to outlier measurements, which must be removed.
Numerical analysis showed that the largest differences in HWS between the mast anemometer and FDWL measurements were attributable to bearing and SV of the FDWL. Therefore, we applied the following outlier rejection criteria: (i) HWS < 2 m/s or HWS > 70 m/s, (ii) bearing = 0 • (FDWL compass issue, see Section 4.1), and (iii) 95th percentile spatial-variation threshold.

Data Screening and Quality Assurance
Now, we discuss our results for outlier criteria (ii) and (iii) which were introduced in Section 3.6 above. Figure 4a shows the HWS difference between the cup anemometer and the FDWL measurements at 83 m as a function of the FWDL bearing. Similar results were found for the two other FDWL heights, 25 and 56 m (not shown). Figure 4b shows a similar plot but for the mast-DWL. As expected, the FWDL bearing covers the full range of motion, whereas the mast-DWL bearing remains mostly fixed at ≈ 3 • . The largest HWS differences occur for bearings equal to zero, which lacks physical consistency and might show an issue in the data acquisition system.  Figure 5 shows the histogram of SV from the FDWL at 25 m and the HWS differences as a function of the SV between the 27-m cup anemometers and the FDWL, respectively. In Figure 5a, the border line between the white and grey shades delimits the 95th percentile of the one-tailed spatial-variation distribution. This percentile corresponds to a SV threshold, SV= 0.055 and SV≈ 0.055 − 0.065 for the other measurement heights of the lidar (data not show). For the mast-DWL, the 95th percentile corresponds to the threshold SV≈ 0.07 (data not shown). The SV thresholds used here for 10 min data were found to be more conservative than those used in previous studies (SV≈ 0.1) [62,63], in which the FDWL measurements were compared against metmast values for quality assurance (QA).
The appropriateness of the SV filtering criterion is re-encountered in Figure 5b, which shows that HWS differences between the FDWL and the metmast measurements started to deviate from the ideal 0-m/s bias (dashed red line) for SVs above approximately 0.05, which justifies our choice for criterion (iii) in Section 3.6.
Regarding the anemometers, Figure 6 shows the HWS differences between the cup anemometers and the lidars as a function of direction. The reduction on the scattering of data after application of the TWD screening criterion is well noticed (black versus grey dots). Quantitatively, after outlier rejection, peak HWS differences were below ±1 m/s at all heights, which shows that the mast effects were effectively removed. The 82-day IJmuiden campaign produced a dataset with 10,833 10 min samples. After the filtering procedure described in Section 3.6, we obtained 8263 10 min samples. Figure 7 shows that both the cup anemometers and the fixed and floating lidars measured very closed HWSs with determination coefficients higher than ρ 2 = 0.996 and slopes of the regression lines between 0.985 and 1.010 at all measurement heights, as expected from previous studies [17][18][19][20][21]. These values are within the KPIs defined by the Carbon Trust's Offshore Wind Accelerator (OWA) [4].

Sensitivity to the Wind Model Parameters
In order to analyze the sensitivity of the parametric wind model presented in Section 3.2 with respect to Equation (9) model parameters, a sensitivity simulator was implemented ( Figure 8). The simulator inputs are the model parameters Obukhov length L, friction velocity u * , and roughness length z 0 , denoted as L 0 , u 0 * , and z 0 0 , respectively. These inputs were perturbed by an intensity factor ±P%(±10%). Then, the perturbed inputs were used to compute the non-neutral wind correction Ψ m z L using Equation (7). Finally, the perturbed wind profile, U(z) + ∆U(z), was simulated as a function of the perturbation intensity on the input variables through Equation (6). Five different atmospheric stability scenarios were considered: very unstable, unstable, neutral, stable, and very stable, based on Obukhov length values denoted by L 0 vu , L 0 u , L 0 n , L 0 s , and L 0 vs , respectively. The error propagated on the model wind profile was found to be insensitive to perturbations in the roughness length z 0 . Therefore, practical implementation of the algorithm embeds Equation (3) into Equation (6) as Thus, the model dimension reduces from three (L, u * , and z 0 ) to two (L and u * ), hereafter called the 2D algorithm. Likewise, slight deviations in the Charnock's constant (Equation (3)) have low impact on the retrieved wind model parameters (L and u * ). Furthermore, as expected from correction function Ψ m ( z L ) (Figure 2), the sensitivity of the model wind profile with respect to the Obukhov length was found to be proportional to 1/|L|. Further discussion is given in Section 4.4 , in line with the statistical results encountered. Figure 8. Sensitivity simulator block diagram. L 0 , u 0 * , and z 0 0 denote the nominal values for L, u * ,and z 0 , respectively. Subscripts vu, u, n, s, and vs stand for very unstable, unstable, neutral, stable, and very stable atmospheric stability conditions, respectively. The non-neutral wind correction block implements Equation (7) and the parametric wind profile model block refers to Equation (6).

Wind Shear Dependence on Dimensionless Stability
The dependence of the wind speed ratio (Equation (8)) on the dimensionless stability computed from both the 2D Obukhov length ( z L ) and from the Ri reference ( z L Ri ) enabled us to quantify the performance of the 2D algorithm. The wind speed ratio was computed between the top and bottom measurement heights of both the FDWL (83 m and 25 m) and the mast (85 m and 27 m). As described in Section 3.4, z = 15.5 m was taken as the reference height. We used 8,263 10 min clean data records (Section 3.6). Additionally, two successive levels of outlier filtering were considered: (i) retrieved Obukhov lengths in the interval (−50 < L Ri < 10) m were rejected because these values were outside the stability limits in Table 3, and (ii) wind speed ratios were histogram-filtered. To achieve the latter, the dimensionless stability range was divided into 0.01-width bins and, in each bin, the wind-speed ratios outside of the µ ± 1σ (with µ being the mean and σ being the standard deviation of the wind-speed ratio distribution in the bin) were rejected as outliers. Figure 9a,b plots FDWL and metmast wind-speed ratios as a function of dimensionless stability parameter z L , respectively. Figure 9c plots metmast wind-speed ratio as function of z L Ri . The FDWL wind-speed ratio as a function of z L (Figure 9a) yielded the best statistical indicators either if all or nonoutlier samples were considered (ρ 2 all = 0.965, RMSE all = 0.029; ρ 2 1σ = 0.988, RMSE 1σ = 0.016, respectively). This was expected because the FDWL-derived HWSs were used to retrieve the estimated Obukhov lengthL using the 2D method, which, in turn, relies on the Businger-Dyer functions. Figure 9b provides tiebreaker proof confirming the successful performance of the 2D algorithm. Thus, in Figure 9b, the mast-derived wind speed ratios are compared to the estimatedL values, and the statistical indicators found (ρ 2 1σ = 0.974 and RMSE 1σ = 0.023) are virtually coincident with those in Figure 9a. Additionally, when comparing the distribution of the points along the horizontal axes in Figure 9b,c ( z L and z L Ri , respectively), figures that share identical mast data, the spread along the z L axis (2D algorithm) is narrower than along the z L Ri (Richardson reference). Again, this confirms the superior performance of the 2D algorithm. Regarding Figure 9c, and despite the clear dependence obtained between the wind-speed ratio and dimensionless stability z L Ri , the huge number of outliers suggests that not all measured wind profiles were acceptably modeled by the MOST model (Equation (6)).
In non-MOST situations, e.g., low-level jet and storm events, the estimates of the 2Dalgorithm may be highly biased. The residual norm between the measured and fitted wind profiles has been used to quantitatively assess the 2D-algorithm accuracy. The residual norm is indicative of the quality of the fitting and is defined as the norm of the error between the parametric and the measured wind profiles. Most of the records showed residual norm values between 0 m/s and 2.21 m/s. Occurrences with a residual norm higher than the residual norm 95th percentile, corresponding to a value of 0.42 m/s, were filtered out (413 records). In the vast majority of cases, the 2D-algorithm is able to converge because most of experimental wind profiles measured by the FDWL are well predicted by MOST. The mast wind-speed ratio plotted as a function of z L Ri . Black and grey dots represent valid and outlier samples, respectively. z = 15.5 m. ρ is the coefficient of determination. Subscripts "all" and "1σ" represent all samples and samples at µ ± 1σ (see body text), respectively. Red line, wind-speed ratio reference model (Equations (8) and (7)). Blue lines, stability classification thresholds from Gryning et al. [51] (see Table 3). vu, very unstable; u, unstable; nnu, near-neutral unstable; n, neutral, nns near-neutral stable; s, stable; and vs, very stable.

Performance Statistics: Friction Velocity and Stability
The performance of the 2D algorithm with regard to friction velocity u * was assessed by comparing the velocity estimated using the 2D algorithm,û * with the two reference velocities described in Section 3.5, namely, the reference u * sonic computed by the sonic anemometer and the mast-derived u * 1−D computed by the 1D method ( Figure 10).  Figure 10a shows that the points (u * sonic ,û * ) were scattered widely. This was further corroborated by numerical analysis, which yielded a linear regression (LR)û * = 0.813 · u * sonic + 0.113, coefficient of determination ρ 2 = 0.731, and RMSE = 0.077 m/s. These comparatively modest indicators are due to the fact that the sonic anemometer was located at a height of 85 m, which may well be beyond the surface layer. Above the surface layer, the proposed model is not valid, which leads to biased estimations [43]. Worse results were attained when considering alternative values for Charnock's parameter, α c , which is a site-dependent parameter here assumed as α c = 0.012 based on previous studies [33,44,65,66] When comparingû * against mast-derived u * 1−D (Figure 10b), a much better agreement is found. This is supported by a LR as good asû * = 0.965 · u * 1−D + 0.022, which is virtually the ideal 1:1 line; ρ 2 = 0.956; and RMSE = 0.031 m/s. The performance of the 2D algorithm for the estimated Obukhov length,L, was tested against the Richardson reference, L Ri , for the whole campaign. For this task, two wellaccepted stability classification criteria in the literature were considered: (i) Van Wijk's (Table 2), and (ii) Gryning's criteria (Table 3).
Van Wijk's classification. Figure 11 compares classification results among the Obukhov lengths estimated by the 2D algorithm,L, and by the Ri estimate, L Ri . The type criteria are provided in Table 2. Figure 11a,b shows similar but not identical results. For example, the hourly evolution of the stable and neutral classes is nearly identical for both estimators. However, this is not the case for very unstable class, which was overestimated by the 2D algorithm to the detriment of the unstable class when compared to the Richardson reference. Similar results can be derived from Figure 11c,d. Thus, the 2D method yielded 33% of very unstable cases, whereas the Richardson reference yielded 25%. Gryning's classification. A similar analysis was conducted for the Gryning classification (Table 3). This criterion essentially breaks down the Van Wijk stability classes into subclasses in order to create room for the near-neutral stable and near-neutral unstable subclasses. Figure 12 compares the classification results yielded by the 2D algorithm,L, and Ri estimate, L Ri . Figure 12a,b shows similar temporal behavior for the very stable and stable classes over the course of the day. The latter class became prominent during the last hours of the day. In contrast, as occurred with the Van Wijk classification, the neutral class and all the unstable classes estimated by the 2D algorithm exhibited larger variations compared to the Richardson reference. When comparingL to L Ri in the pie charts, 5% (L chart) vs. 3% (L Ri chart) of the cases were classified as very unstable, 11% vs. 15% were classified as unstable or near-neutral unstable. In contrast, the following classes remained essentially the same: 31% (L panel) vs. 31% (L Ri panel) for the neutral case, 19% vs. 19 % for the near-neutral stable case, 28% vs. 27% for the stable case, and 4% vs 4% for the very stable class.
(c) (d) Figure 12. Overall campaign stability classification results (Gryning et al. [51] criterion) using the same format as in Figure 11.
Next, we quantified classification performance by rating the one-to-one stability class correspondence betweenL and L Ri . This one-to-one correspondence is expressed via a confusion matrix, in which each matrix row represents the instances in an actual reference stability class (L Ri ), and each column represents the instances in an estimated stability class (L). Figure 13 depict the matrix obtained when considering the Gryning classification criterion. An ideal predictive method would have all instances along the main diagonal of the confusion matrix. Thus, bluish cells represent instances in which the estimated stability class matches the Richardson reference class, whereas reddish cells represent misclassifications. In order to quantitatively assess performance of the 2D algorithm, we define the hit rate (HR) as  The miss rate is simply the complementary function 1 − HR. Figure 13 depicts the confusion matrix considering Gryning's stability classes (Table 3). We found that the estimated classes lay along a band formed by the main diagonal, the first diagonal below this, and the first diagonal above the main diagonal. Thus, an overall HR of 62.59% was obtained and the remaining 29.21% corresponded to the classes estimated by the 2D algorithm that were adjacent to the correct ones given by the Richardson reference. The summary matrix shows a very high miss rate in unstable atmospheres with the HRs ranging from 18.8% to 40.7%. The low HRs attained for the nnu, u, and vu classes occurred due to the low sensitivity of the wind shear model in Equation (6) with respect to the Obukhov length for unstable atmospheres. This is re-encountered in Figure 14a, in which the median wind profiles and related 25th to 75th percentiles are plotted parameterized by stability condition. For the nnu, u, and vu classes, the percentile bars overlap (dark green, red, and magenta traces, respectively), which makes it impossible for the 2D algorithm to discern them. In contrast, "s", "nns", and "n" classes scored HRs of 83.7%, 58.1%, and 69.3%, respectively, on account of their higher sensitivity, which, in turn, led to higher discernability among these classes. It is also worth noting the low HR (35.9%) achieved by the "vs" class. One reason accounting for that is the comparatively short span of Obukhov lengths included in the "vs" class (10 ≤ L < 50) in Table 3. This was a too short span for the NSLQ solver to estimate the Obukhov length at such fine level of accuracy, which eventually led to miss-classification between the "vs" and "s" classes. This is also evidenced in Figure 14a by overlapping percentile bars for these classes.
Analogously to Figure 13, Figure 15 represents the confusion matrix considering the Van Wijk stability classes ( Table 2). An overall HR of 66.07% and a misclassification rate of 33.93% were observed. The estimation performance varied widely, ranging from HRs of 26.2% for the u class to 88.3% for the vs class. Again, the low HRs were due to the low sensitivity of the wind shear model with respect to the Obukhov length in unstable conditions, as evidenced in Figure 14b.  (Table 2) in the same format as in Figure 13.
By comparing Figures 15 and 13 (equivalently, Figures 11 and 12) Van Wijk's classification criterion shows better performance than Gryning's. Although both have similar HRs (66.07% for Van Wijk and 62.59% for Gryning), the Van Wijk confusion matrix (Figure 15a) scores higher HRs than Gryning's in all but the u class, as indicated by the reddish cell in the summary matrix. This, however, warrants some comments. For example, the good match between the estimated and reference stability classes in the pie chart in Figure 12 for the Gryning classification hides a population inversion between the right-and left-adjacent classes to the main diagonal ( Figure 13). This also applies to Figure 15a. Experiments also showed that wind profiles often became indistinctly misclassified into adjacent classes, as shown by the overlapping error bars in Figure 14. A method to mitigate this issue involves collapsing all unstable subclasses into an aggregated class named "unstable". Specifically, this aggregation leads to vs, s, n, and u classes in the Van Wijk classification (Figure 15b), and to vs, s, nns, n, and u classes in Gryning's classification (figure not shown). As a result, a much better agreement between the estimated and reference stability classes was achieved, as indicated by the higher HRs attained in all classes. Quantitatively, the overall HR increased to 73.6% for the Van Wijk classification. Similar results were achieved for the Gryning classification (figure not shown), which scored an HR of 72.6%.
It is important to mention that although we use L Ri as the reference, it cannot be considered as the ground-truth reference for one-to-one comparison withL, because it was not derived exactly as proposed by Grachev et al. [55] due to the absence of turbulent fluxes measurements in the surface layer and the different instrumental measurement heights available. However, it is a good indicative of the method goodness in terms of atmospheric classification. In order to validate the method against a third party, the stability estimated by the 2D-algorithm was compared against the stability estimated by means of virtual potential temperature gradient (dθ v /dz) between 90 and 21 m. Virtual potential temperature was retrieved using temperature, pressure and relative humidity data from the IJmuiden mast via Equation (11). The potential temperature gradient indicator is only able to discern between two stability classes: unstable (dθ v /dz < 0) and stable (dθ v /dz > 0). Considering these categories, HRs higher than 82% were found for the comparison between L and dθ v /dz, further proving the goodness of the 2D algorithm.

Summary and Conclusions
In this study, we suggested an 2D algorithm that accurately estimated both the Obukhov length and the friction velocity, hence, correctly determined atmospheric stability from FDWL measurements at four different heights. We used the spatial variation within the lidar scans as an important filtering parameter of the FDWL ( Figure 5). The parameter is a proxy of turbulence; high spatial variation is associated to nonuniform wind during the scan. The latter can originate either from apparent turbulence (i.e., induced by sea motion) or from true wind turbulence [62]. When comparing FDWL to mast-based wind speeds, we found that the mean winds were biased less than 1.5% (LR, U FDW L = 0.985 · U mast + 0.018) with a determination coefficient, ρ 2 = 0.996, and a RMSE as low as 0.25 m/s. Similar values were found for the fixed DWL on the mast.
The algorithm assumes that the measured FDWL wind profiles follow MOST and the Charnock's relation is used to parameterize the roughness length. Over the sea, with such low roughness, the sensitivity of the roughness on the estimated parameters can be neglected. A parametric sensitivity study showed that unstable wind profiles exhibited lower sensitivity than stable wind profiles to variations of the Obukhov length ( Figure 14).
The 2D-estimated Obukhov length and friction velocity were compared with reference to the mast-derived Obukhov length using the Richardson number, to the sonicanemometer-derived friction velocity, and the 1D-model-derived friction velocity. Thus, the 2D-estimated friction velocity was largely in agreement with the sonic-anemometer-derived and the 1D-model-derived friction velocities, with coefficients of determination of 0.72 and 0.94, respectively ( Figure 10).
We also examined the performance of the algorithm by classifying atmospheric stability in a number of classes. When comparing the relative frequencies of occurrence of each class of the Van Wijk classification, differences occurred between 0% and 10% ( Figure 11) and between 0% and 4% for Gryning's classification ( Figure 12). Notwithstanding these results, analysis through confusion matrices showed HRs of 65.24% for the Van Wijk and 62.63% for the Gryning classification. Higher HRs were attained for stable regimes than for unstable ones due to the lower sensitivity of the 2D algorithm to the Obukhov length in unstable regimes. This issue was addressed by collapsing the unstable sub-classes into a unique aggregated class named "unstable". For the Van Wijk's, this reclassification improved the HR up to 72.9% and up to 72.4% for the Gryning's. The confusion matrix study also showed that although simple pie chart statistics ( Figure 12) indicate a very good one-to-one correspondence between the estimated stability class (througĥ L) and the Richardson reference class (through L Ri ), this correspondence is only apparent because cross-correspondences occur frequently between classes adjacent to the main diagonal (Figures 13 and 15). This effect became more prominent for classes falling in an Obukhov length interval in which Businger-Dyer's correction function had a flatter derivative (i.e., less gradient). Other misclassifications were attributable to anomalous profiles such as those associated with low-level jets. Moreover, the bulk Richardson number methodology used to compute the gradient Richardson number is just a two-point approximation of the derivative of the local wind and temperature gradients in the surface layer. Finally, the performance of the algorithm in terms of stability classification was compared against the virtual potential temperature gradient method acting as a proxy of stability. HRs higher than 82% were encountered, further validating the algorithm performance.
Overall, we aimed to show the potential of FDWLs for offshore wind resource assessment as a standalone instrument and the ability of an algorithm to estimate atmospheric stability from the FDWL wind speeds only. As further steps, we would like to evaluate the algorithm's performance against direct measurements of both momentum and heat fluxes taken over the same range of heights. Additionally, a comparison between the 2D-algorithm and other methods based on MOST relationships is still required.