Using Observed and Modelled Heat Fluxes for Improved Extrapolation of Wind Distributions

Modelling the horizontal and vertical variation of wind speed is crucial for wind energy applications. A model frequently used for this purpose is part of the Wind Atlas Analysis and Application program (WAsP). Here, we modify the model in WAsP to account for local atmospheric stability parameters. Atmospheric stability effects are treated by using the impact of a temperature scale on the geostrophic drag law and the diabatic logarithmic wind profile. Using this approach, wind-direction dependent mean and standard deviation of a surface-layer temperature scale and a mean boundary-layer height scale can be obtained from either numerical weather prediction model output or observations to improve vertical extrapolations of Weibull wind speed distribution parameters. The modified atmospheric stability model is validated at six flat sites in northwestern Europe. The surface-layer temperature scale is available from sonic anemometer measurements at three of the sites. At all sites the temperature scale is also estimated from reanalysis data and from mesoscale model output. The modified model improves the aggregated estimations of power density distributions when extrapolating to nearby locations from 5.2 to 3%, when using the temperature scale derived from either observations or mesoscale/reanalysis output compared to the current model.


Introduction
Knowledge of the change of wind speed with height and horizontal position is critical for many applications. For example, before developing a wind farm, a detailed estimate of the wind climate at each prospected location of a wind turbine has to be known. The annual energy production and its uncertainty are needed to assess the economic viability of the wind farm. In this process, a reduction of 1% in wind speed uncertainty can lead to a 3-5% change in the net present value of a wind farm (Lee and Fields 2021). Dependent on the size and lifetime of the wind farm, such an uncertainty reduction may represent millions of euros.
B Rogier Floors rofl@dtu.dk In a typical onshore wind farm, both microscale flow effects, such as flow accelerations over hill tops, and larger scale effects, such as the boundary-layer height, may impact the flow. Troen and Petersen (1989) combined models for both microscale and geostrophic flow effects and provided a physically sound way of extrapolating wind statistics from a site to nearby locations. The combination resulted in a software called the Wind Atlas Analysis and Application Program (WAsP), which has been used for decades by the wind energy industry for wind resource assessments. The input in their method is a long-term measured wind speed histogram, which is 'cleaned' from the microscale effects of elevation (Troen 1990), internal boundary layers (Sempreviva et al. 1990) and obstacles (Perera 1981), resulting in a generalized wind climate (GWC) over flat, homogeneous terrain. The GWC contains the scale, A, and shape, k, parameters of a Weibull distribution at a number of standard heights and surface roughness length z 0 . The GWC can then be used to predict the flow at a nearby location when microscale corrections at the site of interest are applied.
The orographic sub-model in WAsP generally fails to model the flow in complex terrain, e.g. when terrain slopes > 30% occur in vicinity of the site (Petersen et al. 1998). Therefore, we restrict ourselves in this study to the operational envelope in which the default WAsP model is valid as defined in Petersen et al. (1998).
Atmospheric stability largely influences the behavior of the vertical wind profile in the boundary layer (Businger et al. 1971;Holtslag 1984;Kelly and Troen 2016). Traditionally, the effect of stability on the wind profile is modelled by the logarithmic wind profile with a correction term Ψ that is a function of atmospheric stability, i.e, the diabatic wind profile (Holtslag 1984). The extrapolation of A and k with height depends on the mean and standard deviation of the wind speed. Troen and Petersen (1989) suggested a model that can predict A and k at a given height.
In Troen and Petersen (1989), atmospheric stability effects are treated by determining the impact of varying surface heat flux on the geostrophic drag law and the diabatic wind profile. The effects are then removed from the observed winds and reapplied to the winds at the site of interest. This allows for extrapolation of wind climatologies from onshore to offshore sites, which might have very different stability regimes. Currently, the model relies on 'default' surface heat flux values to estimate the atmospheric stability conditions offshore and onshore, although these are site specific. Therefore, users often manually 'tweak' the default values of the stability model until they match the observed and the modelled wind profile, when there are observations at more than one vertical level. For example, Mortensen et al. (2014) found that for 10 masts in South Africa, the vertical profiles of observed mean winds were well predicted by WAsP when the mean heat flux was varied within the range −60 to 50 W m −2 . However, over flat terrain, the surface roughness length largely affects the wind profile. Therefore, tweaking might potentially compensate other errors in the model setup and thus a more physical approach is needed. Kelly and Gryning (2010) and Kelly and Troen (2016) developed a method to use observed distributions of heat fluxes to better predict the vertical wind profile, but their method did not reach operational use.
In this study, we develop a wind extrapolation method that can use atmospheric stability measures from observations or numerical model output. Since the 90 s, the advent of numerical weather prediction (NWP) models allowed general circulation modelling at increasingly high resolutions. Particularly reanalysis products, such as CFSR or ERA5 (Saha et al. 2014;Hersbach et al. 2020), nowadays provide global information about the flow at scales of ≈ 20-50 km. Alternatively, mesoscale models with a grid spacing of 1-10 km, which use the outputs from a global model as boundary conditions, can provide even higher resolution representations of flow phenomena. The Weather Research and Forecasting (WRF) model is one of the most frequently used mesoscale models because it is open source and has a large user community. Hahmann et al. (2021) showed that using the WRF model output compared to that of ERA5 led to much smaller errors in the predicted wind distributions compared to measurements by analyzing wind speed observations from 19 masts in South Africa. One may therefore expect that such high resolution models are important to estimate local stability parameters.
Despite that both mesoscale and reanalysis data might be too coarse to capture microscale flow effects, we hypothesize that they can still provide better local estimates of some of the flow properties than using constants. Particularly, Troen and Petersen (1989) assumed atmospheric stability and boundary-layer height to be constant due to lack of measurements and simulation output at the time. The effect of baroclinicity was examined in Floors et al. (2018) by estimating the geostrophic wind shear from reanalysis data, which improved the accuracy of wind extrapolation with height when compared to measurements within the entire boundary layer.
In Sect. 2 we present a modified stability model that predicts the variation with height of the mean and standard deviation of the wind speed using classic atmospheric stability scaling parameters. Section 3 presents the sites and measurements that are used to validate the model. In Sect. 4 we present the details of the inputs and the setup of the new and the original WAsP stability model. In Sect. 5 we then investigate the behaviour of the temperature scale, which is a new model input, and use cross-predictions at six different sites to validate the model. In Sect. 6 we conclude the study and provide ideas for future research.

Theory
The wind profile in the planetary boundary layer (PBL) is largely influenced by the surface sensible heat flux H (W m −2 ), i.e. the vertical flux of heat that is transported from or to the surface by turbulence. Using Reynolds decomposition, representing turbulent fluctuations using primes and denoting a time-average (typically over 10-30 min) with an overbar, the vertical kinematic heat flux (K m s −1 ) is given by, where w is the vertical wind component, θ the potential temperature, c p the specific heat of dry air and ρ the air density. It is common to use the virtual kinematic heat flux w θ v instead of w θ to account for the impact of humidity, because humidity impacts the buoyancy of air (Sjöblom and Smedman 2002). Then, Eq. (1) becomes: where E is the moisture flux (kg m −2 s −1 ). Troen and Petersen (1989) accounted for the effect of atmospheric stability on the geostrophic drag law by considering its impact on the friction velocity u * under barotropic atmospheric conditions, by using a first-order expansion of the geostrophic drag law in H (see Appendix), 78 R. Floors et al. where G is the geostrophic wind speed, f the Coriolis parameter, g/T the buoyancy parameter and: where κ is the von Kármán constant, A 0 and B 0 denote the geostrophic drag law constants in neutral conditions (Blackadar and Tennekes 1968) and A μ and B μ denote the same 'constants' but varying with respect to the Monin-Kazinski stability parameter (e.g. Zilitinkevich and Esau 2005), where L is the Obukhov length. Equation (3) is then used to calculate the deviations on u * from its value under neutral conditions caused by both a mean heat flux H off and a fluctuating part H rms (see Appendix for the derivation). The subscript off refers to an offset value caused by the mean, while rms denotes a root-mean-square value caused by the fluctuating part. In this way, we avoid computation of the wind profile for each instance in time, e.g., each 10-min period. Instead, we can obtain a correct climatological behaviour of the wind profile, which has an enormous computational advantage.
Because H is highly dependent on wind speed (Mahrt 2017), here we apply a similar first-order expansion of the gestrophic draw law to that in Troen and Petersen (1989) but in terms of a normalized heat flux T * (see Appendix), so Eq. (3) becomes, Note that some authors define T * with the same sign as L and opposite of that of H , but here we omit the minus to allow easy comparison with the original formulation that is expressed in terms of H . In stable conditions, T * is also known as the surface-layer temperature scale θ * (Monin and Obukhov 1959): whereas in unstable conditions we can use a velocity scale w s for normalization, Various forms for w s have been proposed, as it should account for combined mechanically and convectively driven boundary layers and obey the free-convection limit, in which the scaling velocity w s = w * = ( g T w θ v h) 1/3 (McBean 1976). In first-order closure schemes the velocity scale at the top of the surface layer is often specified as: where c φ is a constant in the expression for the non-dimensional wind shear. By normalizing with this velocity scale, boundary layers that are characterized by a combination of wind shear and buoyancy lie within the range of fully mechanically to buoyancy driven PBLs. The transition between stable, Eq. (7), and unstable, Eq. (8), is continuous. In ERA5 c φ = 15 was adopted (ECMWF 2016), whereas Hong et al. (2006) adopted c φ = 8. The reason for adopting Eq. (9) is further discussed in Sect. 5.2. As in Troen and Petersen (1989), an expression for the reversal height z m is used, which corresponds to the height where first-order effects of surface heat flux variations vanish and therefore the vertical profile of the wind speed variance has a minimum, i.e.: where Ψ denotes the integrated dimensionless wind shear (Holtslag 1984). As shown in the Appendix, this equation results in: where a = 4.8 (Cheng and Brutsaert 2005). Because the Ψ -function is only valid in the surface layer, we discuss the influence of the limited boundary layer height later in this section. According to Eq. (11), z m is a weak function of z 0 . This was confirmed by the analysis of observations of vertical profiles of the Weibull shape parameter at offshore, coastal and onshore sites . We note that z 0 is used in connection with the geostrophic drag law and in WAsP it is determined by using an exponential weighting of the logarithm of the roughness lengths of all land-cover segments within each wind direction sector (Troen and Petersen 1989). Consequently, z 0 is a variable determined over ≈ 10 kilometers and is not the roughness length at the point of interest.
According to Troen and Petersen (1989), the deviation of the wind speed from the neutral wind speed U 0 at z m can be calculated as, where Ψ is computed as a function of z m /L off and z m /L rms instead of z/L, where L off is the Obukhov length corresponding to T * off and L rms the Obukhov length corresponding to F rms T * rms , where F rms is a form factor. Troen and Petersen (1989) suggested F rms = 0.6. Δu * rms /u * 0 denotes the results of Eq. (6) when T * is replaced with T * rms multiplied with F rms . A range of Ψ functions have been reviewed by Högström (1988). Here we still adopt the formulation from Jensen et al. (1984) for unstable conditions, but that of Cheng and Brutsaert (2005) for stable conditions, where b = 2.5. Holtslag (1984) demonstrated that a Ψ -function that is linear in z/L, as proposed in Troen and Petersen (1989), does not yield good results for wind profile modelling in very stable conditions. Equation (14) also has the advantage that it uses a from Eq. (11) so we do not introduce new parameters. We assume that Eqs. (13) and (14) are valid at z = z m . Equations (6) and (12) are used to obtain a vertical profile of deviations from a neutral state of both the mean and standard deviation. The deviation of mean wind speed U (z) from its neutral value U 0 (z) is specified as, whereas the deviation of standard deviation σ U (z) from its neutral value σ U 0 (z) is according to Troen and Petersen (1989) specified as, where Δu * off and Δu * rms denote the results of Eq. (6) when evaluated using an offset (mean) or root-mean-square value of T * , respectively, and F(z) is a profile function that will be discussed below. In Eq. (15), the first three terms in the large brackets are the same as in Troen and Petersen (1989) and the last term is adapted from the wind profile model from Ghannam and Bou-Zeid (2021). The latter model reduces to the logarithmic wind profile in the limit G f z 0 ≡ Ro → ∞, but accounts for the effect of the boundary-layer height and baroclinicity. Here we are studying the effects of atmospheric stability only and therefore the impact of baroclinicity is ignored. We note that the influence of baroclinicity on the wind profile can be pronounced and the model can correct for this using A μ and B μ values (Eq. (4)) that vary with respect to the mean geostrophic wind shear vector Ghannam and Bou-Zeid 2021). Troen and Petersen (1989) derived a profile function from the first order expansion of U (z) given by: where λ limits the profile function for increasing z/h: where p = 1.5. This term is required because Ψ in Eq. (12) is too large at heights above 100 m under stable conditions (see for example Fig. 4 in Holtslag (1984)). Troen and Petersen (1989) estimated h using the relation: where h * = 0.165, in agreement with the review from Hess (2004). Zilitinkevich and Esau (2005) showed that the lapse rate just above the boundary layer and baroclinicity also impact h. Because Eq. (15) only specifies the deviation from the neutral state, we obtain the wind speed at a certain height z 2 from another height z 1 as: where the ratios U (z 2 )/U 0 (z 2 ) and U (z 1 )/U 0 (z 1 ) are computed from Eq. (15). For the standard deviation we have: where the ratios σ U (z 2 )/σ U 0 (z 2 ) and σ U (z 1 )/σ U 0 (z 1 ) are computed from Eq. (16). Note that in Eqs. (20) and (21), we have only explicitly specified the z dependence, but one can also apply a varying surface roughness. Equations (20) and (21) are thus used to correct a wind climate input given by the Weibull parameters A 1 and k 1 to a wind climate at another height. The mean and standard deviation of a Weibull wind speed distribution are given by: respectively, where Γ is the gamma function. At the new height z 2 , the new k 2 is found by using a simple bisection algorithm using the corrected first and second moment of the input Weibull distribution. We can now readily obtain the power per unit area swept by wind turbine blades, P (also referred to as the power density) from the corrected A and k parameters by using the mean cube of the Weibull distribution, We further illustrate Eqs. (20) and (21) in Sect. 5.1. Using the power curve of a wind turbine, which is typically supplied by the manufacturer and which is assumed to be a function of the hub-height wind speed only, the wind speed U at hub height (z 2 in Eq. (20)) can be used to calculate the annual energy production.

Measurements
Measurements at six sites in northwestern Europe are used in this study (see Table 1), because tall masts with sufficient levels of measurements and at least three years of good quality data are scarce. Two of the sites are offshore, two are coastal, and two are inland (see Figs. 1 and 2). The height of the measurements at offshore and onshore sites refers to heights above mean sea level (amsl) and above ground level (agl), respectively. The measurements at IJmuiden (offshore) were obtained with a continuous-wave wind lidar (Zephir 300) from 90 to 315 m, which was installed at the platform (Kalverla et al. 2017). This dataset was supplemented with measurements from cup anemometers and wind vanes at 27 and 58 m from the platform's meteorological mast. The mast is located ≈ 85 km from the coast (Fig. 1). The data are available from 1 January 2012 to 1 January 2016. During this period the data recovery rate was 83.9%. We only use measurements where both the lidar, cup anemometers and wind vanes simultaneously report valid values. For all the other sites reported from now on, we also use only periods where all sensors report valid measurements simultaneously.
The FINO3 platform is located ≈ 80 km offshore west of Jutland, Denmark. The wind speed is measured at 50, 70 and 90 m by cup anemometers mounted on 3 booms in different directions to avoid flow distortion. The wind direction was measured by a sonic anemometer at 60 and 100 m. This dataset is described in more detail in Peña et al. (2015). After 2014, several wind farms were constructed in the vicinity of the mast, leading to large wake effects. Therefore we only use measurements from 1 November 2011 to 1 November 2014. No data from the calendar year 2012 were used, because it contained large gaps in the time series. After removing this year, the recovery rate was 83.3%.
The wind speed measurements at Cabauw are described in Van Ulden and Wieringa (1996) and are corrected for flow distortion using measurements from all three booms dependent on the wind direction (Bosveld 2018). Here we use measurements from 1 January 2015 to 1 January 2018 (recovery rate 97.2%). In addition, we use measurements from the sonic anemometer at 10 m, which provides high-frequency time series at a sampling rate of 10 Hz (Bosveld 2018). The area around the mast is characterised by meadows and scattered higher vegetation and villages (see Fig. 2). The measurements at Høvsøre are described in Peña et al. (2016) and have been corrected for flow distortion using the method in Lindelöw-Marsden et al. (2010). Because there are episodes where the cup anemometers were frozen due to icing, measurements where the standard deviation within a 10-min period in wind speed and wind direction were smaller than 1 × 10 −5 m s −1 and 0.15 • , respectively, are filtered. A Metek USA-1 sonic anemometer at 10 m is used to determine the surface-layer fluxes at a sampling rate of 20 Hz. To separate the mean wind field from the turbulence fluctuations, a double rotation technique is used (Wilczak et al. 2001). This method rotates the coordinate system such that the v and w components are zero, while the u component is equal to the value of the mean wind speed over the flux averaging interval. In addition, a cross-wind correction according to Liu et al. (2001) is applied. The crosswind-correction accounts for contamination of the sonic temperature by wind components perpendicular to probes. Data from 1 January 2015 to 1 January 2018 are used here (recovery rate 85.0%). The land cover around the site is dominated by open grassland with scattered villages, while to the east there is a larger forested area. The North Sea coastline is located 2 km to the west of the mast, resulting in a distinct internal boundary layer at this site (Floors et al. 2011).
The Østerild site has two masts of 250 m height: the northernmost mast is located ≈ 3.5 km south from the North Sea coastline, whereas the southernmost mast is located ≈ 4 km south of the northern mast (Peña 2019). To obtain a wind speed and direction time series, the same quality control procedure was performed as at the Høvsøre site. The terrain around the masts is rather heterogeneous with nearby forests, meadows and villages. A Metek USA-1 sonic anemometer installed at 37 m height was used to compute the surface-layer fluxes. We use the data from 1 January 2017 to 1 January 2019. Because we require all sensors on both masts available simultaneously, the recovery rate was 80%.
The Risø measurements are from a mast that is situated at a peninsula. Most of the peninsula is gently rolling land with the 117-m meteorological tower erected on a 6-m hill. The mast is located on a meadow with rows of trees and buildings north of the mast and to the south a shallow bay. The measurements from 1 January 2012 to 27 December 2016 are used, with a recovery rate of 97.9%.

Methods
In the current study, we will compare four different profile model configurations implemented within the WAsP model chain. All version of WAsP require the following inputs: -A histogram of observed wind speed and direction at a given location and height -A roughness length vector map -An elevation vector map To vertically extrapolate wind statistics, the current profile model in WAsP (labelled as WAsP 12.7 default throughout the rest of this paper) uses default values of H off of −40 and −8 W m 2 over land and sea, respectively. Thus, the current version assumes slightly stable conditions on average. For H rms , 100 and 30 W m 2 are adopted over land and sea, respectively. Sea and land are treated separately because atmospheric stability conditions are usually largely different due to the heat capacities of water and land surfaces. In coastal areas both offshore and onshore stability conditions might affect the site conditions, so the deviations from neutral conditions, U /U 0 in Eq. (15) and σ U /σ U 0 in Eq. (16), are found by linear interpolation based on the fraction of land and sea for a certain wind direction sector. Wind speed histograms are generated for each height and site (e.g., mast or lidar) from the 10-min time series of mean wind speed and wind direction, by discretization using a bin width of 1 m s −1 and 30 • , respectively.
To create a roughness length map, 40 × 40 km areas surrounding the sites were extracted from the CORINE database (Copernicus Land Monitoring Service 2022). This dataset stores polygons of land cover throughout Europe, classified in 44 different landcover types. To each of these types, we assign a characteristic z 0 based on the revised values recommended in Table A4 in Floors et al. (2021) (Fig. 2). Because both masts at Østerild are located near the forest edge, a displacement height of 10 m is added to landcover classes that contain forests. The landcover polygons are then converted to roughness lines with a left and right-hand side roughness length, which is the input required for the roughness analysis in WAsP (Floors 2023). This conversion was done using the QGIS software (QGIS Development Team 2022).
The landcover to roughness length approach is used extensively in wind engineering (Indasi et al. 2016). For the offshore sites, we assume homogeneous terrain in all directions and z 0 = 0.0002 m as recommended by Wieringa (1993), although z 0 depends on wind speed (Charnock 1955).
The elevation vector map for all sites in Denmark was obtained from the national digital terrain elevation (DTM) model (Styrelsen for Dataforsyning og Effektivisering 2022), which is obtained from lidar scans. The vector lines had a equidistant spacing of 2.5 m. At Cabauw, the raster DTM model from the Netherlands was used (Nationaal Georegister 2022). It has a resolution of 5 m and was converted to vector lines using the QGIS software. All the inputs described so far are the same for all the WAsP versions within this paper.
The new version of WAsP that uses Eqs. (6)-(14) to model the wind profile is labelled WAsP 12.8 hereafter. In addition to the inputs listed before, this new version requires sectorwise (i.e. binned by wind direction) mean h * (Eq. 19), T * off and T * rms values.
Because measurements of heat fluxes with a sonic anemometer are rarely available and we are aiming to develop a method that can be easily applied anywhere in the world, the ERA5 data were adopted as a practical choice to estimate T * in absence of measurements. The ERA5 data are generated using version CY41R2 of the Integrated Forecast System (IFS) (Hersbach et al. 2020). The model has a grid spacing of ≈ 31 km. Regular latitude/longitude grid data are downloaded at 0.25 • (≈ 30 km) resolution (Copernicus Climate Change Service Climate Data Store (CDS) 2017). We find the ERA5 grid point over land and over sea nearest to the sites. From these grid points, time series of H , E, u * , h, θ and ρ were extracted to compute T * and h * (Eqs. 7-8 and Eq. 19, respectively). Because high wind speeds are much more important in the calculation of P (see third power in Eq. 24), we compute T * off , h * off and T * rms conditional on the highest 50% of wind speeds at 100 m from the ERA5 data. This percentage was chosen because it closely corresponds to the frequency of occurrence of wind speeds higher than the observed average speed, which is used for Weibull fitting in WAsP (Troen and Petersen 1989) and most important for calculating the annual energy production of a wind turbine. The ERA5 wind direction at 100 m is used to create the wind-directions bins. The version of WAsP that uses stability inputs from ERA5 is labelled WAsP 12.8 ERA5. We use the data from 1 January 2011 to 1 January 2021 and we assume this period is representative of the climatological average (Floors 2022).
In addition, output from the WRF model version 3.8.1 was used (Skamarock et al. 2008). The numerical setup is extensively described in Hahmann et al. (2020), but it is briefly summarized here. It uses three nested domains with a resolution of 27, 9 and 3 km, respectively. The simulation were reinitialised from ERA5 boundary conditions every 7 days and ran for the period 1989 to 2018. Each seven day period was started 24 hrs early and this spin-up period was removed from the simulations. The WRF setup included the MYNN PBL scheme Niino 2006, 2009) and the Janjic (2019) surface-layer scheme. We use the output from 1 January 2010 to 1 January 2019. The methodology to derive T * off , h * off and T * rms is the same as that described above for the ERA5 data. This version with stability inputs from the WRF model output is labelled WAsP 12.8 WRF.
Finally, to obtain a measured T * , a sonic anemometer was available at Høvsøre, Østerild and Cabauw (see Sect. 3). At those sites we can estimate u * = u w 2 + v w 2 1/4 and w θ v in Eq. (7) from the measurements. Sjöblom and Smedman (2002) showed that w θ v can be approximated well from w θ measured by a sonic anemometer. Furthermore, h is required in Eq. (8) but it is not measured directly at any of the sites. Therefore we estimated it by linearly interpolating the PBL height from the half-hourly WRF time series to the 10-minute time series from the sonic anemometer. In addition, h and u * from the WRF output were used to estimate h * . This version with blended inputs from observations and WRF is labelled WAsP 12.8 Obs.. It represents our best attempt to vertically extrapolate Weibull parameters given the current measurements. All versions of WAsP described above can predict the distribution of wind speeds at a certain height using the wind distribution at another height or site, i.e. a cross-prediction. This sector-wise predicted wind distribution is then compared to the observations at that height. A Weibull fit is done to the observed histogram of wind speed to obtain the sectorwise Weibull parameters. Because WAsP is usually used to extrapolate winds to hub height, which is generally equal or higher than the height of the input observations, here we are going to perform upward cross-predictions only.
We can find P from the Weibull parameters A and k for N d sectors as: where φ is the sector frequency and we assume ρ = 1.225 kg m −3 because most of the sites are located near sea-level and we are interested in the effects of atmospheric stability only.
In this work, the metric we use to evaluate the results is based on the relative error of power density, where P is obtained from the modelled (mod) or observed (obs) Weibull-derived values from Eq. (25). The statistic later used in Sect. 5.6 is the mean absolute relative error, |ε P |, because it contains information about both the mean error (bias) and the scatter around the mean (width of the distribution).

Profiles of Weibull Scale and Shape Parameters
To illustrate the effect of T * on Eqs. (20) and (21), we assume that at 40 m the wind speed is described by a Weibull distribution with A = 7.4 m s −1 and k = 2.1 and set z 0 = 0.05 m, which are used as input to the model. The vertical profiles of the Weibull scale parameter (Fig. 3a) and shape parameter (Fig. 3b) are shown for a range of T * off values from −0.15 (stable) to 0.05 K (unstable), whereas T * rms = 0 K for all profiles. At 40 m all profiles collapse, because this is the input height. For T * = −0.15 K, the Weibull A strongly increases with height, whereas for unstable atmospheric conditions (T * = 0.05 K) the increase with height of A is much smaller. This is directly related to the much smaller variation of the stability correction as a function of z/L in Eq. (13) used in unstable conditions compared to that in Eq. (14) used in stable conditions. The dotted lines in Fig. 3 show Eq. (20) with T * off = 0 K and f = 0, i.e., the atmosphere is 'truly neutral' according to the terminology used in Zilitinkevich and Esau (2005), and thus a logarithmic wind profile is observed at all heights. When the last term in Eq. (15) is non-zero and T * off = 0 K, A increases slightly more with height compared to the case in which the term equals zero, i.e. a 'conventionally neutral' atmosphere (Zilitinkevich and Esau 2005). This correction term has also been studied by Gryning et al. (2007), Kelly and Troen (2016), Ghannam and Bou-Zeid (2021) and accounts for the effect of limited boundary layer height and the turning of the wind with height. It is required because traditional Monin-Obukhov similarity theory is not valid above the surface layer.
In Fig. 3b when T * off = 0 K, k becomes constant with height when using Eq. (16) with f = 0, while for stable and unstable conditions k increases and decreases with height, respectively.
The profiles in Fig. 3a, b behave the same as the logarithmic wind profile with a stability correction term, e.g. those in Holtslag (1984), but with Eq. (16) we can also model the behaviour of the k parameter. Gryning et al. (2016) found that k increases up to the reversal height and then typically decreases. This behaviour mostly depends on the T * rms term in Eq. (16) and therefore we also show a typical range of T * rms values from 0 to 0.15 K in Fig. 3c, d.
In Fig. 3c the A parameter increases more with height the higher the T * rms values. This is because of the shape of the Ψ function, whose dependence on z/L is higher under stable than unstable conditions; thus the A profile is displaced to the stable side despite of T * off = 0 K (see Eq. (12)).

Fig. 3
Weibull A and k parameters as a function of height for a range of T * off (panels a and b) and T * rms (panels c and d). In panels a and b T * rms = 0 K, whereas in panels c and d T * off = 0 K. The solid lines use Eqs. (20) and (21) to get A and k, and the dashed line represents the same equations but with the last term in Eq. (15) equal to zero. The input height for the Weibull parameters is 40 m In Fig. 3d k increases up to the reversal height and then decreases upwards due to the impact of |F(z)| in Eq. (16). Higher values of T * rms lead to higher vertical gradients of k. When T * rms = 0 K, k is constant with height.

Distributions of T *
Here, we aim to investigate whether measured or modelled T * show similar distributions and can therefore be used as inputs in Eqs. (15) and (16). To compute the terms Δu * off /u * 0 in Eq. (15) and Δu * rms /u * 0 in Eq. (12) and (16), we first need to obtain T * off and T * rms and use those in Eq. (6). If T * does not follow a symmetrical distribution, we will not be able to calculate a representative mean and rms value from a time series of T * because a long tail of the distribution can strongly impact the calculated mean and rms. By using T * from Eq. (9) in the calculation of L for unstable conditions, the method is close to exact near neutral. However, it deviates for larger heat fluxes, where it has the desirable (albeit somewhat arbitrary) effect of reducing the influence of the very unstable tail of the distribution. For example, the distribution of H or θ * has highly asymmetrical distributions and could therefore not be directly used.  Figure 4a-c shows the histograms of observed time series of T * for Cabauw, Høvsøre and Østerild, where surface-flux measurements were available. T * obtained from sonic anemometers are compared with those obtained from ERA5 and WRF output to investigate how well they represent atmospheric stability conditions. A Gaussian distribution with the mean and standard deviation of T * is also shown.
The observed histogram of T * resembles a Gaussian distribution for these three sites, although the number of unstable conditions skews the histogram slightly. T * off is smaller (∼ −0.05 K) at Cabauw than at Østerild and Høvsøre (∼ −0.02 K). The standard deviation is slightly higher at Cabauw (∼ 0.06 K) than at Østerild and Høvsøre (∼ 0.05 and ∼ 0.04 K, respectively).
The Cabauw mast is more inland (see Fig. 2a), and is characterized by a stronger diurnal cycle and relatively more stable conditions. Because Høvsøre and Østerild are located less than 7 km from the shoreline, they present a more maritime climate with less frequent very stable and unstable conditions. The shape of the bin-wise distributions of T * from observations and the WRF model is rather similar. The Gaussian distributions based on the WRF model output and the observations show that T * off and T * rms for Cabauw and Høvsøre are rather similar, but for Østerild T * off from WRF is more negative than the observed one. Also the spread of T * from the WRF model is higher than observed values. The observed histogram and that obtained from ERA5 output are very similar at Cabauw and Østerild, but at Høvsøre ERA5 has a lower T * off .
The WRF model has higher resolution than ERA5 and should estimate surface heat fluxes more accurately. However, ERA5 is globally available for a long period. When using ERA5, a long-term reference climate for the life time of a wind turbine (∼ 25 years) can therefore be estimated without requiring simulations.
Some of the deviations between the distributions in Fig. 4 may be because the WRF or ERA5 output and the observations are not concurrent (see Sect. 4). However, year-to-year variation of T * off values, which is the stability parameter to which the wind profile is most sensitive, was generally less than 10%. Because all observational datasets have at least 3 years of data, the deviations from the climatological mean are expected to be small.

The Impact of Wind Direction on T *
Because the wind resource at a site is highly dependent on the wind direction and the upstream conditions, it is common to model the wind profile in a sector-wise fashion. In this section, we investigate if and how T * depends on the wind direction sector for the onshore site Cabauw and the offshore site FINO3. Figure 5a shows the mean and standard deviation of T * at Cabauw as function of wind direction binned within 30 • , using the wind direction from the observations or model output that was nearest to 100 m. The northerly wind direction is centered at 0 • . T * is negative for all sectors. For wind directions between 60 and 150 • , T * is most negative. This is in agreement with previous studies were those sectors were characterized by an abundance of stable conditions and low-level jets (Baas et al. 2009).
In Fig. 5a, the standard deviation of T * based on the values computed using the observations is generally higher than those from the WRF model and the ERA5 output. This might indicate that these models lack the resolution and variability to fully represent the variability of T * , which can be avoided using high resolution and/or WRF-LES (Floors et al. 2018;Peña et al. 2021). Figure 5b shows nearly the same results as those in Fig. 5a but for the offshore FINO3 site. Unfortunately, no observations were available for any offshore site. Based on the ERA5 output, the mean T * is mostly around 0 K, which is in agreement with the default heat fluxes in Troen and Petersen (1989), who assumed that the atmospheric stability over the sea is closer to neutral than over land. The mean atmospheric stability at FINO3 is mostly unstable, i.e. in Fig. 5b T * off > 0 K for most sectors, but for the most frequently observed southwesterly winds (240 • ), T * off < 0 K. This emphasizes that it is difficult to interpret the all-sector mean profiles without knowing the wind frequency distributions. In most sectors, the standard deviation of T * is much smaller than over land, which is also in agreement with the current default rms values in Troen and Petersen (1989); 100 and 30 W m −2 over land and sea, respectively. Atmospheric stability is controlled by different processes. Over land, the diurnal cycle generally determines the atmospheric stability. During the day, heating of the surface results in predominant unstable stability conditions, whereas during night, the surface cools down and stable conditions are common. Due to the large heat capacity of water, diurnal modulations of stability are much smaller. Atmospheric stability conditions 90 R. Floors et al. offshore are mostly determined by the temperature difference between the sea surface and the air flowing over it (Sathe et al. 2011).
For the wind directions 180-360 • , results from using the outputs of the WRF model and ERA5 show rather similar mean and spread of T * values. Under easterly winds, results from the WRF model output clearly show more unstable conditions than those from ERA5. This might be also related to resolution: for the ERA5 data there are only few grid points to resolve the land-to-sea transition, whereas this process can be better resolved by the WRF model. Both for Cabauw and FINO3, southerly winds tend to show negative T * values (stable atmospheric conditions). Generally on the Northern hemisphere, southerly flow tends to have more stable conditions, because warmer air is advected from lower latitudes over a cooler sea or land surface. This can be seen by contrasting winds from the northerly and southerly wind direction sector in the 10-yearly mean T * off dataset that is supplied along with this paper (Floors 2022).

Illustration of Vertical Profiles
To illustrate that the proposed profile model Eqs. (15) and (16) is also observed in the wind measurements, we compare the model for a single sector from the Cabauw site. This site is the most homogeneous site with heat flux observations analyzed here and therefore microscale effects due to internal boundary layers or elevation changes are small. Winds are most frequently coming from the 210 • wind direction sector at this site, so we analyze this sector with sufficient long-term observations. The land cover in this direction is dominated by grassland and z 0 ≈ 0.05 m. Fitting a Weibull distribution to the wind speed distribution from observations at 40 m gives A = 7.4 m s −1 and k = 2.1, which are used as input to the model. We use the input from the observed T * off and T * rms values at Cabauw. For the 210 • sector, T * off = −0.05 K and T * rms = 0.06 K (see Fig. 5). The comparison between the modelled and observed A and k is shown in Fig. 6.
In Fig. 6a, the WAsP 12.8 Obs. modelled profile shows good agreement with the A values derived from the observations at all heights. In Fig. 6b the modelled k profile shows good agreement below 40 m, but fails to capture the observed increase of k at greater heights. This is mostly because the reversal height z m seems too low compared to observations. z m is sensitive to the constants d A/dμ and d B/dμ (through c in Eq. (11)), which following Troen and Petersen (1989) are equal to 0.2. We investigated the effect of using different values for d A/dμ and d B/dμ using the data from all masts (not shown), but improvements were found at only some of the masts.
In addition, we also investigated the impact of the boundary layer scale h * : it was found that this parameter had a minor impact on the profile of A, but for very low T * a low h * is required to limit the high Ψ -term in equation (16) via Eq. (18). Increasing h * by 10% only led to a 1% change in wind speed at 250 m and even less on other heights.

Validation of Vertical Profiles
For the validation of the profiles at all sites we focus on the power density, Eq. (25), because of its importance for wind energy applications. The procedure illustrated in Fig. 6 is repeated for each observed wind climate, so we can obtain profiles of the power density as a function of height for all sectors combined (Fig. 7). As stated in Sect. 4, we only use upward extrapolations, so there is usually a good agreement at the first level, because there is only a single sample available (i.e. from the lowest vertical level, see Table 1). Note that for sites with many vertical levels with observations, there are more possible cross-predictions. Our aim in this section is to analyze the performance for individual sites, whereas in Sect. 5.6 we look at the aggregated performance of the four different models.
At Cabauw (Fig. 7), WAsP 12.7 default does not agree very well with the mean P derived from the wind measurements at 80 m. WAsP 12.8, which uses T * derived from the observations, agrees best with the observed P. At 200 m, WAsP 12.8 ERA5 shows an overprediction. At Høvsøre (Fig. 7b), all versions of WAsP 12.8 perform rather similar, whereas that using the default atmospheric stability model largely overestimates the power density above 100 m.
At IJmuiden (Fig. 7c), the observed P increases less rapidly with height than at the onshore sites (P ≈ 900 W m −2 at the lowest level), because the roughness length offshore is very low. The deviations at the highest level are therefore large compared to the onshore sites, but in relative terms the errors are similar (see Table 2). Both WAsP 12.7 default and all versions of WAsP 12.8 overpredict P above 140 m. The predictions from WAsP 12.8 WRF agree best with the observations at IJmuiden, but the power density is anyway overestimated at the highest heights.
At the Østerild masts, we can also use cross-predictions to predict winds and power density from, e.g. the southerly mast to the northerly mast. Therefore there are also model results at the lowest observational level. Figure 7d-e shows that WAsP 12.7 default overpredicts P at the highest vertical levels at Østerild. The strong power density 'shear' indicates that the default value H off = 40 W m −2 is too high (too stable) for the site, similarly to the findings at Høvsøre and IJmuiden. T * values computed from the WRF model output at Østerild were also too stable (Fig. 4), which explains the over predicted 'shear' from the WAsP 12.8 WRF model at the southerly mast.
At Risø and FINO3, there are not many heights with available wind observations, so it is difficult to conclude which WAsP version performs best. At Risø, WAsP 12.7 default shows strong changes of power density with height that do not appear in the observations, nor in any WAsP 12.8 versions. At FINO3, all WAsP versions overpredict P at 90 m by about 20 W m −2 (≈ 2%).

Error Metrics at All Sites
Finally we aim to provide error metrics for the four different versions of WAsP. To have a unbiased impression of the error metrics, it is important to use the same number of crosspredictions at each mast; the Østerild site has two masts with seven different heights, resulting in 98 cross-predictions, which are ∼60% of the total number of possible cross predictions. Therefore, the error metrics would be dominated by the results from Østerild only. Instead, we therefore use bootstrap sampling and randomly select three cross-predictions from each site, which is the minimum number that is available at each site. The sampling is repeated a thousand times and subsequently averaged, thereby obtaining an estimate of |ε P | for each model version which is not biased towards a specific site.
The results of the bootstrap sampled aggregate |ε P | as well as the mean |ε P | for each site from all cross predictions are shown in Table 2. Using the atmospheric stability derived from either ERA5 output, WRF model output or the observations results in lower |ε P | at most sites when compared to the results using WAsP 12.7 default, as well as when all sites are combined.  Looking at the column with all samples, there is no clear improvement of using WAsP 12.8 with WRF output instead of ERA5 output, despite the higher spatial resolution of the WRF simulations used here. |ε P | from all sites is 3.0% and 2.9% for the ERA5 and WRF versions, respectively. More sites in areas with strong mesoscale effects, such as at lower latitudes or near mountain ranges, are needed to assess whether using mesoscale outputs adds value for this type of usage when compared to using reanalysis data only.
At Cabauw, where observations of T * were available, |ε P | of WAsP 12.8 Obs. are slightly lower than those using the ERA5 or WRF-driven versions. Due to the absence of offshore sonic anemometers, offshore stability conditions are estimated with model output, which might increase the error in the predictions. This may also affect the coastal sites (i.e. Østerild and Høvsøre), where both upstream and downstream atmospheric stability conditions are needed to describe the speed-up factors from Eqs. (15) and (16). We also evaluate the mean absolute relative errors in wind speed: as expected these were lower than |ε P |, but generally showed similar improvements as those related to power density. WAsP 12.7 default shows a relative error of 2.2%, while WAsP 12.8 versions have relative errors within the range 1.3-−1.5%. As already mentioned, such a reduced uncertainty represents significant monetary value when developing wind farms.
For future studies a more representative stability estimate can be obtained by combining several grid cells from the NWP output, for example by assigning weights according to similarity in roughness or elevation (Bechmann et al. 2020). In this study, the validation of the model is performed using measurements covering a small geographical area (see Fig. 1). Further testing of the model was therefore performed on 109 masts deployed over 4 continents and in many different climate regimes. The relative errors in power density from this larger validation decreased from 6.8 to 5.3% when using WAsP 12.7 and WAsP 12.8 ERA5, respectively. The report with all these statistics is available for download (Floors 2022).
The temperature scale in the modified model was obtained from observations at three of the sites and all parameters were obtained for all six sites from the ERA5 reanalysis and WRF model output. The six sites were located onshore, offshore and near the coast in flat terrain. T * off , T * rms and h * were derived from the nearest grid point both onshore and offshore, to be able to model both conditions and to interpolate for coastal conditions.
The mean and rms values of T * from observations of surface fluxes were well approximated by those derived from the ERA5 reanalysis and WRF model output at three sites were eddycovariance measurements were available. The mean and rms of T * were found to be dependent on wind direction, with dominating stable atmospheric conditions under southerly winds. Offshore, the mean of T * was closer to zero and the rms was much smaller than the values over land, supporting the choice of the default values of the atmospheric stability parameters in the current version of WAsP.
After gathering required wind climate histograms and vector map inputs of elevation and roughness length, cross-predictions from one position to another were used to evaluate the performance of the current version of WAsP 12.7 with the default atmospheric stability model and new versions WAsP 12.8 that use the new stability model with the additional atmospheric stability parameters derived from either the ERA5 reanalysis, WRF model output or observations, respectively. On average over all sites, the new model improved crosspredictions, decreasing the relative errors in power density from 5.2 to 2.9−3.0%.
Using observations from sonic anemometers further decreased relative absolute power density errors at one of the three sites where fluxes are measured. However, the difference in errors was small compared to using ERA5 or WRF-derived stability parameters. For practical use in wind resource applications, using globally available modelled data is therefore attractive and used for the forthcoming version of WAsP.
For future studies, we will investigate the applicability of the current approach in numerical wind atlases, where instead of using an observed wind histogram we use a modelled one as input to the WAsP model chain. Furthermore, the distributions of u * and Ψ in Eq. (11) are currently only characterized by a mean and rms, so some detail of the distribution may be lost. An approach where they are estimated directly from the observed or modelled values may improve this.
Also the spatial heterogeneity of atmospheric stability needs to be investigated, for example by using a simple energy balance model that can resolve microscale changes in atmospheric stability. We will also investigate the role of atmospheric stability within the free atmosphere and its impact on the scaling height u * / f and on the geostrophic wind shear (Zilitinkevich and Esau 2005). To validate these new models, extensive measurements are required both near the surface and higher up in and above the boundary layer, particularly in areas closer to the equator where terms using the ratio u * / f become problematic. Finally, the method also has to be tested in more complex terrain, where there might be spatial variability in the stability parameters.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. deviation at this height. Similar to Troen and Petersen (1989), we assume that the net result can be modelled by a positive H (T * ) proportional to the rms deviation from the average offset. The unknown constant involved is the form factor (0.6), which depends on the H (T * ) distribution.