Evaluation of idealized large-eddy simulations performed with the Weather Research and Forecasting model using turbulence measurements from a 250 m meteorological mast

We investigate the ability of the Weather Research and Forecasting model to perform large-eddy simulation of canonical flows. This is achieved through comparison of the simulation outputs with measurements from sonic anemometers on a 250 m meteorological mast located at Østerild, in northern Denmark. Østerild is on a flat and rough area, and for the predominant wind directions, the atmospheric flow can be considered to be close to homogeneous. The idealized simulated flows aim at representing atmospheric boundary layer turbulence under unstable, neutral, and stable stability conditions at the surface, which are statistically significant conditions observed at Østerild. We found that the resolved fields from the simulations appear to have the characteristics of the three stability regimes. Vertical profiles of observed mean wind speeds and direction are well reproduced by the simulations, with the largest differences under near-neutral conditions, where the effect of the subgrid-scale model is evident on the vertical wind shear close to the surface. Vertical profiles of observed eddy fluxes are also well reproduced by the simulations, with the largest differences for the three velocity component variances under stable stability conditions, although nearly always within the observed variability. With regards to turbulent kinetic energy, we find good agreement between observations and simulations at all vertical levels. Simulated and observed velocity spectra match very well and show very similar behavior with height and with atmospheric stability within the low-frequency interval; at the effective resolution, the simulated spectra show the typical drop-off of finite differences. Our findings demonstrate that these idealized simulations reproduce the characteristics of atmospheric stability regimes often observed at a high turbulent and flat site within a direction sector, where the air flows over nearly homogeneous land.


Introduction
For many applications and, in particular, for wind energy, we would like to characterize the long-term site conditions, i.e., firstand second-order statistics of the three-dimensional velocity vector, at a number of locations and vertical levels within a given area, so that we take into account all relevant motion scales of the atmosphere. For such a purpose, a multiscale modeling 20 approach is needed, in which one starts by downscaling the large scales of atmospheric motions, from, e.g., reanalysis and global forecasts, to the regional or the mesoscales using a numerical weather prediction (NWP) model, and continuing down We focus our analysis on the accuracy of resolved and modeled atmospheric flow parameters by WRF-LESs through comparison with observations. We include in the comparison vertical profiles of mean wind speed and direction, and velocity variances and covariances, as well as turbulence spectra at various vertical levels. 60

Selection of flow cases
The measurements at Østerild are described in detail in Peña (2019). For this study, we only use statistics based on 10-min periods of measurements performed with Metek USA-1 sonic anemometers deployed at 7, 37, 103, 175 and 241 m on the meteorological tower. Figure 1-left shows the distribution of atmospheric stability conditions close to the surface (using the sonic anemometer measurements at 37 m) within the close-to-homegenous sector at Østerild. We do not use the measurements 65 at 7 m for this purpose as these are strongly influenced by the local topographical inhomogeneities (forest trees) near the mast (Peña, 2019).
We are interested in modeling three types of ABLs: near-neutral (|z/L| ≤ 0.05) but referred to as neutral for simplicity hereafter, unstable (−0.5 ≤ z/L ≤-0.2), and stable (0.2 ≤ z/L ≤0.5), where z is the height and L the Obukhov length (Obukhov, 1946). As illustrated, the surface stability conditions are predominately neutral, however, unstable and stable conditions are 70 frequently observed. The observed ensemble-average surface heat flux was very close to zero (−6.2 × 10 −5 K m s −1 ), 0.0948 K m s −1 , and −0.0276 K m s −1 , under neutral, unstable and stable atmospheric conditions, respectively. This results in z 0 = 0.2492 m. For the unstable and stable conditions, the neutral value is too high and by using z 0 = 0.0992 m together with the atmospheric stability correction based on Monin-Obukhov similarity theory (MOST) (Monin and Obukhov, 1954), the prediction of the normalized vertical wind profile is very good for the first 100 m for the three main ABL regimes. 80 We therefore choose these z 0 -values for the LESs below. MOST predictions are given as and here we choose ψ m = 0, −4.7z/L, and −(3/2) ln([ Gryning et al. (2007) for neutral, stable, and unstable conditions, respectively. latitude (57.0489 • ). The time step used for the simulations was 0.1 s. All simulations were performed using the SGS model of Deardorff (1980) with the prognostic equation for the subgrid turbulent kinetic energy.

Simulations
Simulations were initialized assuming a dry atmosphere. For the neutral ABL simulation, the initial temperature was kept constant (289.5 K) up to 700 m and then an inversion of 10 K km −1 was imposed. Such an inversion strength is a common 4 https://doi.org/10.5194/wes-2020-119 Preprint. Discussion started: 8 January 2021 c Author(s) 2021. CC BY 4.0 License.  (Pedersen et al., 2014;Mirocha et al., 2018). The height of the inversion was chosen to 105 be lower (we expected the neutral ABL to slowly grow with time) than the value for the ABL-height estimation using the parametrization of Rossby and Montgomery (1935), where C = 0.1-0.5. For a nearby site with similar climatology to Østerild, Peña et al. (2010a) found similar ABL heights by comparing the estimations from Eq. (2) using C = 0.15 with those from observations of aerosol backscatter profiles under 110 near-neutral conditions. Using the ensemble-average u * -value from the near-neutral observations (0.69 m s −1 ) and Eq. (2) with C = 0.15, z i = 845 m. The initial u x -and u y -velocity components, aligned with the x-and y-axis, respectively, were kept constant throughout the ABL, with values of 14 and 0 m s −1 , respectively. For the unstable ABL simulation, the initial temperature was kept constant (289.5 K) up to 700 m (since we expected the unstable ABL to grow faster with time than the neutral ABL) and then an inversion of 4 K km −1 was imposed. Such an inversion strength is a common choice for unstable 115 ABL modelling (Pedersen et al., 2013;Mirocha et al., 2018). The initial u x -and u y -velocity components were kept constant throughout the ABL, with values of 8 and 0 m s −1 , respectively. For the stable ABL simulation, the initial temperature was kept constant (289.5 K) up to 100 m and then an inversion of 10 K km −1 was imposed. Such an inversion strength and level are common choices for stable ABL modelling (Kosović and Curry, 2000;Muñoz-Esparza and Kosović, 2018). Using the ensemble-average u * -value from the stable observations (0.36 m s −1 ) and Eq.
(2) with C = 0.12 as suggested for stable 120 conditions in Peña et al. (2010a), z i = 353 m. The initial u x -and u y -velocity components were kept constant throughout the ABL, with values of 14 and 0 m s −1 , respectively. The initial u x -velocity for the three simulations was chosen so that it was close but slightly higher than the observed ensemble average of U at each of the stability regimes at the highest vertical level.
MOST was applied at the surface through the in-built WRF surface-layer scheme (option 1 in WRF's namelist), although a modification of the open-release scheme was performed to maintain simulations dry. The neutral and unstable ABLs were 125 simulated during 20 and 6 h by imposing a constant surface heat flux w Θ of 0 and 0.0948 K m s −1 , respectively, mimicking the ensemble-average observed values. For the stable ABL, imposing a surface heat flux is problematic because it does not guarantee a stable solution for the computed u * -values (Basu et al., 2008). Therefore, imposing a cooling rate boundary condition at the surface is a choice (Kosović and Curry, 2000). We apply a rate of −0.25 K h −1 , run for 24 h, and check the ability of the model to reach the observed heat flux and friction velocity at the surface.

130
The simulations were performed using periodic boundary conditions in both horizontal directions. We output selected variables for a vertical column in the middle of the domain every 1 s and instantaneous values of those variables every 1 h for the positions in the whole domain.

Results
For the three ABL regimes under study, we present an analysis of the simulated transient outputs (Sect. 3.1), an overall picture

Statistics on transient simulation outputs
We need to extract WRF-LES outputs to perform the comparison with the observed statistics at Østerild. The choice of the time to extract LES statistics depends on the type of boundary layer. The analysis is mostly made by performing moving averages 145 over 600-s windows based on the 1-Hz outputs of given variables. We estimate the height of the ABL z i as that in which the maximum of the vertical gradient of potential temperature occurs. u * is an output of the WRF model, which is computed within WRF's surface-layer scheme, the subscript). It is seen that turbulence was triggered slightly before 2 h. u * behaves similarly to e sgs,1 ; the former directly depends on the resolved velocity closest to the ground. Although with nearly steady moving averages, the estimated z i slightly increases during the simulation after ≈5 h. slightly increases, and so does z i . It is also observed that the wind speed becomes supergeostrophic during the simulation with a maximum between 10 and 11 h. We therefore chose this hour for computing the neutral ABL WRF-LES statistics.

Unstable conditions
Figure 4 illustrates the same information as Fig. 3 but for the unstable ABL simulation with the correspondent imposed positive heat flux at the surface. In Fig. 4-left, it is seen that turbulence was triggered much earlier compared to the neutral 160 ABL simulation. u * is generally lower than the values of the neutral ABL simulation mostly because of the lower z 0 imposed. z i increases faster with time compared to the neutral ABL simulation and is above 1000 m for the latest 3 h. e sgs,1 is lower than that of the neutral ABL simulation because the geostrophic forcing is significantly lower than in the neutral case.
In Fig. 4-right, it is clear the effect of the positive surface heat flux on the temperature profile, which explains the behavior of z i . The wind speed is, as expected, much more constant with height above ≈100 m compared to that of the neutral ABL. To 165 compute the unstable ABL WRF-LES statistics, we select the output within 3-4 h because z i is very close to 1000 m, and so it is slightly higher than that of the neutral ABL simulation, and the mean wind speed is very constant (slightly below 7 m s −1 ) and subgeostrophic up to about 1100 m, as expected due to the stable layer at the top, and increases about 1 m s −1 within the next 150 m reaching the geostrophic value. show the closeness of the latter to the observed surface heat flux at Østerild. In Fig-5-left, it is seen that turbulence was triggered slightly before 2 h similar to the neutral ABL simulation. After ≈8 h, the moving average of u * becomes relatively steady and slightly higher than the observed value at 37 m. After turbulence is triggered, z i slightly increases until its moving average 175 reaches a value ≈300 m at 8 h and remains nearly steady afterwards. U 1 , u * , and e sgs,1 reach nearly stationary state after approximately 8 h. About the same time, the ABL reaches stationary state as evident by the values of z i . The moving average of w Θ decreases nearly linearly until it reaches the observed value at 37 m a little before 8 h and remains rather steady and generally slightly higher than the observed value thereafter. Note that we verified in the time series output that the surface potential temperature decreased at the imposed cooling rate of −0.25 K h −1 (not shown).

180
In Fig. 5-right, we clearly see the effect of the cooling rate on the temperature profile. The initial imposed inversion quickly 'disappears', a strong inversion develops with time within the range ≈275-400 m, and the initial inversion strength of 10 K m −1 is recovered thereafter. As for the neutral ABL simulation, the wind speed becomes supergeostrophic with the nose at the height where the strong inversion starts. We select the range 15-16 h to compute the stable ABL statistics.

Instantaneous resolved fields
185 Figure 6 shows instantaneous cross-sections of U along the x-z plane at the y-direction midpoint and along the x-y plane at z ≈ 100 m for the three types of ABL. Similar to Mirocha et al. (2018), we find elongated low-speed structures along the   wind shear due to the tendency of the specific SGS model to overpredict the dimensionless shear . For the three cases, the mean of the simulations is always within the observed variability, which is larger than that of the simulations.

Wind speed and direction profiles
RMSEs of the normalized simulated wind speeds also reflect the qualitative findings:   RMSEs are lowest under unstable conditions (0.52 • ) and under stable the RMSE is rather low (1.47 • ). As for the wind speed, the mean of the simulations for the three cases is always within the observed variability. The simulated variability is clearly highest under unstable conditions.

Turbulent fluxes and kinetic energy 220
For the comparison with the measurements, we need to estimate the total variances and covariances from the WRF-LESs; thus we need to account for the resolved and the subgrid stresses. The total stress is given as where τ res ij is the resolved stress, i.e., −ρu i u j , and τ sgs ij the subgrid stress. The latter can be computed as,

225
where τ dev ij is the deviatoric part of the subgrid stress, which is an output of the SGS scheme in the WRF model. Note that the resolved turbulent kinetic energy e res is u i u i /2 (with implicit summation) and so the total turbulent kinetic energy is the sum of both SGS and resolved terms. Since the observed statistics are computed so that u is aligned with the wind direction at each vertical level for each 10-min period, we need to rotate both the simulated velocities and stresses to align them with the simulated direction for each simulated vertical level. conditions. Observations and simulations of the velocity variances are in relatively good agreement, particularly at ≥100 m, and the largest apparent differences are found for the uw-covariance, where the resolved value is higher than that of the observations at ≥100 m and is lower than the observed value at ≤100 m. For both the uw-and vw-covariances, the SGS term seems to overcompensate the flux when compared to the first observed level. For the u-variance, the resolved term is already 235 close to the observed value and accounts for most of the total variance, whereas close to the surface the SGS term strongly aids both v-and w-variances in matching the observed fluxes.
For unstable conditions (see Fig. 10), the apparent bias between observations and simulations is in line with that for neutral conditions for both variances and covariances. For the variances, the apparent bias is slightly higher at ≥100 m and lower at lower levels, whereas it is generally higher for the uw-covariance and lower for the vw-covariance when compared to the 240 results under neutral conditions. Note that the observed variability of fluxes is also higher under unstable compared to neutral conditions.
For stable conditions (see Fig. 11), the apparent bias between observations and simulations is generally the highest among all stability conditions for the three velocity variances, and it is low at the two observed levels closest to the ground. Note that the observed normalized u-and v-variances do not decrease much with height compared to their behavior under neutral  only. For unstable and stable conditions, this also occurs for specific vertical levels.  Figure 12-left shows the ratio of the SGS to the total turbulent kinetic energy as function of height for the three types of ABLs. As expected, the percentage of the SGS term in the total is highest for stable conditions and lowest for unstable conditions. Further, more than 10% of the energy comes from the SGS term below 120 m under stable conditions, which is about 40% of the ABL height. For neutral and unstable conditions, the SGS term contributes more than 10% of the total energy 260 within 10% and 3% of the ABL height, respectively. when compared to neutral (0.412 m 2 s −2 ), and unstable (0.415 m 2 s −2 ) conditions. As shown, the observed power spectra is very well behaved at all vertical levels with an inertial subrange slope close to −5/3 following Kolmogorov (1941). From the lowest frequencies up to ≈0.1 Hz, both simulated and observed spectra are 275 very close, which explains the good agreement between simulated and observed velocity variances in Sect. 3.4. From ≈0.1 Hz, a drop-off of the velocity spectra appears, which is typical of finite difference and discretization schemes (Skamarock, 2004).

Turbulence spectra
The effective resolution is ≈ 7∆x = 105 m, which corresponds to a frequency of ≈0.1 Hz within the wind speed range 9-12 m s −1 . It can also be observed that at 37 m, the drop-off occurs at a frequency lower than 0.1 Hz and this drop-off frequency slightly increases with height, as expected.

280
Under unstable stability conditions (Fig. 14), the observed velocity spectra also approaches the inertial subrange slope of −5/3. Compared to the results for neutral conditions, it is also clearer both when looking at the simulated and the observed spectra that turbulence becomes more isotropic the higher the vertical level and the more unstable the atmosphere is, as the three velocity spectra become close to each other.
Under stable stability conditions (Fig. 15), we also find the −5/3 slope on the observed velocity spectra. When looking 285 both at the simulated and observed spectra, we can see a clear distinction between the velocity spectra compared to unstable conditions as turbulence is more anisotropic; by fitting the spectral three-dimensional turbulence model of Mann (1994) to the observed velocity spectra and uw-cospectra from the sonic anemometer measurements at Østerild, Peña (2019) found a distinct lower turbulence anisotropy the more unstable the surface conditions were, whereas neutral conditions appeared to be slightly more anisotropic than stable conditions. Within the low frequency range and under stable conditions, we can notice more 290 flattened spectra compared to that under unstable and neutral conditions. The premultiplied spectra, i.e., f S i (f ) (not shown), peaks at higher frequencies under stable compared to unstable conditions, which translates into turbulence length scales that are larger under unstable compared to stable conditions (Peña et al., 2010b).

Conclusions and discussion
To assess the ability of high-fidelity simulations, such as LES, to reproduce the behavior of winds and turbulence within the first 295 hundreds of meters of the ABL, which is useful, e.g., for the siting of wind turbines, we need to try, first, to isolate the effects of physics parametrizations and forcing, and, second, to analyze high-quality measurements of both wind and turbulence at several heights. The reason for the former is that such parametrizations and forcing conditions influence the behavior of turbulence and so it is difficult to differentiate their effects, which are accounted for, e.g., in real-time simulations using mesoscale models, from those inherent to the abilities of LESs. Here, by using wind and turbulence statistics, and velocity spectra computed from 300 sonic-anemometer measurements on a 250-m mast over the predominant wind direction at Østerild, Denmark, we demonstrate that idealized WRF-LESs reproduce well the observed wind and turbulence characteristics, which resemble canonical flow of typical ABL regimes (unstable, neutral, and stable).
Comparison with observations reveals that under the three ABL regimes, the vertical profiles of normalized wind and di- Simulated and observed velocity spectra match very well within frequencies lower than that correspondent to the effective resolution, which explains the good agreement between simulated and observed velocity variances. Such a good match is 325 found both under the three ABL regimes and the vertical levels examined. As expected due to the nature of the WRF model, the velocity spectra shows a drop-off close to the effective resolution and so it is only the observed spectra the one that approaches the −5/3 slope within the inertial subrange. Both simulated and observed velocity spectra show that turbulence is more anisotropic the more stable the ABL and the closer to the surface; the more sheared the flow, the more anisotropic the turbulence.

330
Regarding the assumptions made for the simulations we carried out, it is appropriate to reiterate that these are idealized simulations. As such, the initial conditions may not represent observed cases. Observations of the ABL height, and observations of vertical profiles of both potential temperature and water vapor mixing ratio within the extent of the ABL are not available at Østerild. The three cases considered here are all characterized by relatively weak surface heat fluxes and strong shear, i.e., they are shear driven. Therefore, assuming a dry atmosphere, i.e., that the moisture effects on the structure of the atmospheric surface 335 layer are small, is a good approximation in the three cases. Since these are idealized simulations, the initial potential profile is well mixed up to 700 m for the neutral and convective boundary layers, and 100 m for the stably stratified boundary layer.
Capping inversions develop naturally due to surface heating or cooling, while the overlaying inversion in the free troposphere is specified. The overlying inversion (10 K km −1 ), which we used for the neutral and stable ABLs, corresponds to the standard adiabatic lapse rate for dry atmosphere. These choices are commonly used in idealized simulations as well as a value lower 340 21 https://doi.org/10.5194/wes-2020-119 Preprint. Discussion started: 8 January 2021 c Author(s) 2021. CC BY 4.0 License. than 5 K km −1 for the unstable ABL. In all three cases, the ABL height evolves during the simulation based on the combined effect of shear and buoyancy forcing.
Given that the comparisons between the outputs of the idealized simulations and the observed statistics are rather good, in future studies we want to explore the ability of a WRF-LES-based multiscale modeling system to simulate in real-time the ABL at Østerild and other sites in which high-quality measurements are also available, in more typical unsteady operating Additional validation of these approaches using a similar framework to that applied herein will further establish WRF's value to wind energy applications.
Another important element of multiscale atmospheric simulation involves downscaling through the "gray zone" or "terra 355 incognita" scales (e.g., Wyngaard, 2004). New approaches, based on scale-awareness (e.g., Honnert et al., 2011), explicit three-dimensional turbulence transport , and explicit filtering and reconstruction (Simon et al., 2019) represent promising pathways, but also require further evaluation in wind energy applications.
The WRF model's applicability over steep terrain, a known issue when downscaling due to topographic features being better resolved, is likewise being extended, using both higher-order numerical methods (Arthur et al., 2020a) as well as immersed 360 boundary methods (Lundquist et al., 2012;Arthur et al., 2018). These methods have likewise not been adequately evaluated in relation to wind energy relevant flow information.
The framework of the present study can be used to assess the utility of the WRF model in these above described settings to improve wind energy simulations in a broader range of real-world environments and operating conditions, for which smaller-scale flow information, including turbulence, in relation to environmental and meteorological variability, is invaluable 365 to supporting the continued expansion of the wind energy industry.
Code and data availability. The numerical outputs were generated with the open-source WRF model (https://github.com/wrf-model/WRF).
Both the observational and simulated data intercompared in the manuscript, as well as the input files for the WRF model simulations are available at: https://figshare.com/s/c02c0954051d67f17992 (Peña, 2020).