WRF-LES Simulation of the Boundary Layer Turbulent Processes during the BLLAST Campaign

Abstract

A real case long-term nested large eddy simulation (LES) of 25-day duration is performed using the WRF-LES modelling system, with a maximum horizontal grid resolution of 111 m, in order to explore the ability of the model to reproduce the turbulence magnitudes within the first tens of metres of the boundary layer. Sonic anemometer measurements from a 60-m tower installed during the Boundary Layer Late Afternoon and Sunset Turbulence (BLLAST) field campaign are used for verification, which is focused on the turbulent magnitudes in order to assess the success and limitations in resolving turbulent flow characteristics. The mesoscale and LES simulations reproduce the wind speed and direction fairly well, but only LES is able to reproduce the energy of eddies with lifetimes shorter than a few hours. The turbulent kinetic energy in LES simulation is generally underestimated during the daytime, mainly due to a vertical velocity standard deviation that is too low. The turbulent heat flux is misrepresented in the model, probably due to the inaccuracy of the sub-grid scheme.

1. Introduction

A good representation of the atmospheric flow field over complex and non-homogeneous terrain requires numerical modelling using spatial resolution of the order of hectometres or less, broadly referred as large eddy simulation (LES), to include the land use and terrain modifications on the atmospheric boundary layer (ABL) flow field. Moreover, in order to obtain a complete representation of the physical processes involved, a different range of scales needs to be resolved, from small motion scales to mesoscale and synoptic scales. Dynamical downscaling including multi-scale simulations down to LES domains which include resolved turbulence interactions is still an ongoing research topic.

In the Weather Research and Forecasting (WRF) model, turbulent stresses and fluxes are usually modelled using 1D PBL parameterisations, which assume horizontal homogeneity, an approximation that can be done if the grid size is relatively large (e.g., few kilometres to tens of kilometres). In this approach, the equations of the flow field are resolved using the Reynolds-averaged Navier–Stokes (RANS) technique where turbulent motions are fully parametrised [1]. Instead, if the horizontal cell grid size decreases to a few km or tens or hundreds of metres, the assumption of horizontal homogeneity is not valid and the effect of horizontal gradients of turbulent stresses and fluxes should be included [2]. Turbulence is explicitly resolved using the large eddy simulation (LES) approximation, where the large eddies are resolved and the small eddies are modelled in a sub-grid scheme. The cross-scale coupling system WRF-LES [3] allows the simulation of the flow field in a wide range of scale processes, from mesoscale to microscale, being able to resolve the mean quantities and turbulent stresses and fluxes. However, it is known that in the coupling model domain configuration, there will always be a range of grid resolutions where certain processes are neither sub-grid nor resolved, the so-called “grey zone” or “terra incognita” [4]. As recently suggested by [5,6,7], comparison and intercomparison projects of real long simulation time periods are needed in order to establish the progresses in the cross-scale simulations and identify the necessities for the future modelling strategies.

During the last decades, researchers have used LES approximation to study specific problems such as air pollution dispersion in urban street canyons [8,9,10], flow distortion caused by wind turbines [11,12,13], large fires [14,15] or cloud and precipitation processes development [16,17,18,19]. Many other studies have used real case LES to simulate atmospheric flows in complex terrain areas [20,21,22] where grid sizes of few hundreds of metres are needed to resolve the terrain elevation and LES approximation is required. Other authors suggested the need to use LES to resolve small-scale turbulent processes in complex terrain such as downslope windstorms [23], low-level turbulence and rotors [24,25], precipitation dynamics [26,27,28,29] or other small scale processes. Recently, the authors of [30,31] used WRF-LES with domains nested within mesoscale domains and examined the resulting flow field. However, none of them regarded the turbulent quantities, nor explored the ability of the simulation to reproduce the fluxes, variances and turbulence intensity. In addition, a number of previous LES studies are conceived as idealised theoretical studies (e.g., [32,33,34,35]), allowing to isolate specific atmospheric phenomena, while others are designed as real case studies, which need to be compared with proper high-resolution observational data sets. Some recent case studies consider simulation periods of one or few days (e.g., [22,36,37]) which allows sensitivity tests varying the model configuration but restricts the verification process and hampers the evaluation of the general quality of the simulations. In this paper, we report a large eddy simulation using the WRF-LES system spanning 25 days during a field experimental campaign.

In order to explore the accuracy of WRF-LES simulations in reproducing turbulent quantities, we use realistic LES simulations of 111-m horizontal grid size here and we compare the results with the observational data set from the Boundary-Layer Late Afternoon and Sunset Turbulence (BLLAST) campaign [38], which took place during early summer 2011 in France, north of the central French Pyrenean foothills. Therefore, the objective of this paper is to compare a long-term (25-day) real case WRF-LES simulation with observations including mean and turbulent quantities in the first tens of metres above ground, exploring the ability of the model to reproduce the magnitude, time evolution and associated energy of the atmospheric flow field.

The structure of the rest of the paper is as follows. The area of study, observations and model simulations configuration are described in Section 2 together with the methodology used to carry out the comparison. Results and discussion are given in Section 3, evaluating the mean and turbulent quantities. Finally, in Section 4 conclusions are formulated along with final remarks.

2. Site, Data and Methods

2.1. Site and Tower Measurements

The BLLAST site is located on the Plateau de Lannemezan, north of the central French Pyrenean foothills. Locally, it is a relatively flat terrain area but it is influenced by complex terrain features, mountains and valleys in its surroundings. In particular, the site is about 11 km north of the base of the first range of the Pyrenees mountain massif with peaks of around 1500 to 2000 m a.s.l. and is about 40 km far from the highest peaks of the mountain range, exceeding 3000 m a.s.l. (Figure 1). The Pyrenees mountains have a large influence on atmospheric flows of the region both at mesoscale and synoptic scale, as highlighted in a number of numerical meteorological case studies-see for example [39,40,41,42].

A 60 m tower called Valimev Tower was placed at the surface-site 4 (ss4) in Site 1 defined in the BLLAST campaign [38] at 43.1242${}^{\circ }$ N, 0.3625${}^{\circ }$ E (Figure 1). The Valimev Tower sensors were located at 29.4 m, 45.8 m and 61.4 m (named hereafter $z_{30}$, $z_{45}$ and $z_{60}$, respectively), equipped by 3D sonic anemometers (

Table 1. Tower levels, sensors and data period.

Level Name	Height (m)	Sensors	Data Period
$z_{30}$	29.4	Campbell Csat3 3D sonic anemometer	14 June 2011–8 July 2011
$z_{45}$	45.8	Gill master pro 3D sonic anemometer, wind vane
$z_{60}$	61.4	Campbell Csat3 3D sonic anemometer	15 June 2011–8 July 2011

), which are able to measure the three wind components, u, v and w and the temperature at 10 Hz frequency. The data set for this tower starts on 14 June 2011 for the $z_{30}$ and $z_{45}$ levels, and slightly later for the $z_{60}$ level, between 10:00 and 11:00 UTC on 15 June 2011, and finishes on 8 July 2011 for the three levels. The calculated standard deviations of the three wind components and the averaged heat flux with a 5 min temporal averaging are used. We applied the corrections and filters indicated in [43]. Additionally, sonic measurements were not considered during rainfall periods or when ${\sigma }_u$, ${\sigma }_v$ or ${\sigma }_w$ were greater than 2 m s⁻¹, mostly before or after rainfall periods. In summary, the sonic anemometer measurements obtained from the Valimev Tower during the BLLAST campaign enables to evaluate the model simulations during 25 days (24 days for z${}_{60}$).

2.2. WRF-LES Simulations

A long-term simulation is run with the WRF model version 3.6.1 [44] following a configuration previously developed and tested at 85 sites all over the world for wind resource application [45]. Five nested domains centred at Site 1 (see Section 2.1) are configured. The outermost and innermost domains have a horizontal resolutions of 9 km (D1) and 111 m (D5), respectively. The relationship between domains is a ratio of 1:3 as generally used. Simulation is divided in two parts. First, a mesoscale run from 9 to 1 km with two-way nesting feedback enabled. This part is run employing so-called “mesoscale mode” (MESO), which means applying the most-commonly used Reynolds-averaged Navier-Stokes (RANS) decomposition. Second, a microscale run from 333 m to 111 m by using the mesoscale output as initial and lateral boundary conditions (LBC). These simulations are performed in “microscale mode” using the LES approach instead of a PBL scheme (see

Table 2. WRF model set-up. The number of grid points in x, y and z is marked with nx, ny, nz.

Domain	nx × ny × nz	Grid	Mode	PBL Scheme
D1	99 × 99 × 38	9 km	RANS	MYJ
D2		3 km	RANS	MYJ
D3		1 km	RANS	MYJ
D4		333 m	LES	-
D5		111 m	LES	-

). A two dimensional Smagorinsky closure is applied in MESO domains [46]. For LES domains a three-dimensional LES turbulence closure is used as a sub-grid scheme, including a prognostic equation for TKE [46,47]. One challenge in meso-micro scale coupling is the non-turbulent state of LBCs. In order to fix this issue, the perturbation potential method of LBC is applied [48,49]. The model top is set at 50 hPa, and a 5 km damping layer is included in order to prevent wave reflection. A fixed number of vertical levels (38) are used in all domains, with the first level located at around 8 m above ground level and 15 levels included in the first 1000 m of the atmosphere. Although several papers using WRF-LES in recent literature settle a greater number of vertical levels (e.g., [6,22,35,37]), our choice of 15 levels within the first km is a compromise choice between between computational cost and simulation performance for the long-term period of 25 days. Simulation tests increasing vertical levels (from 15 to 22 vertical levels in the first km of atmosphere) for a 4-day simulated period (1–4 July 2011) did not reveal much better results (see Appendix A).

The WRF physics package used here includes the rapid radiative transfer model (RRTM) for longwave radiation [50], the Dudhia parametrisation for shortwave radiation [51] and the Noah land-surface parametrisation [52]. Domains D1, D2 and D3 use the revised MM5 surface layer parametrization [53] and the Mellor–Yamada–Janjić (MYJ) model [54] for boundary-layer parametrisation, a widely used local scheme that has a prognostic equation for the TKE. Cumulus are not parameterised in these domains. Domains D4 and D5 use WRF Single-Moment 3-class scheme microphysics parametrisation [55] and the Monin-Obukhov similarity theory for surface layer parameterisation [56]. Static data of topography elevation and land use are adapted from global databases to the WRF preprocessing system (WPS). The Shuttle Radar Topography Mission (SRTM) [57] with 3 arc-seconds resolution (∼90 m) global coverage is used for the topography and the the GlobCover [58] for land use database.

The simulated period spans for 25 days, from 14 June 2011 to 8 July 2011, when the BLLAST campaign took place. The initial and boundary conditions are given by ERA5 [59] with horizontal resolution of 30 km, 38 vertical levels and temporal resolution of 1 h. For this study, mesoscale and microscale outputs at 1 km and 111 m are considered. Mesoscale and microscale outputs have a different postprocessing. In the first case, WRF output is set to 10 min. In the microscale case, outputs are saved every 0.25 s and a 10-min average and standard deviation for the main fields are computed later.

2.3. Quantities and Statistics

Mean and turbulent quantities are evaluated at three levels of the tower, $z_{30}$, $z_{45}$ and $z_{60}$. The evaluated mean quantities include the horizontal wind speed and wind direction. On the contrary, mean temperature from the sonic anemometers are not considered, as this variable can be underestimated due to the CSAT firmware, because the speed of sound is also underestimated, specially in high winds [60]. The verified turbulent quantities are the standard deviations of the three wind components (${\sigma }_u,{\sigma }_v$ and ${\sigma }_w$), gathered with the turbulent kinetic energy (TKE) calculated as $TKE=0.5({\sigma }_u^2+{\sigma }_v^2+{\sigma }_w^2)$ and the heat flux ($F_H$) expressed as $F_H=\overline{w^{\prime }{T}^{\prime }}$.

The model verification has been performed at 10 min time steps, a temporal resolution deemed appropriate for this type of study. Three verification scores have been considered, namely the mean bias (MB) (Equation (1)), root mean square error (RMSE) (Equation (2)) and regression coefficient (R) (Equation (3)), given by:

$$MB=\frac{1}{n}\sum _{i=1}^n{m}_i-o_i$$

(1)

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{(m_i-o_i)}^2}$$

(2)

$$R=\frac{{\sum }_{i=1}^n(m_i-\overline{m})(o_i-\overline{o})}{{\sum }_{i=1}^n{(m_i-\overline{m})}^2}$$

(3)

where $m_i$ and $o_i$ are, respectively, the model and the observation values at time step i and n is the total number of time steps.

Contingency tables were used to evaluate the horizontal wind speed according to the following verification scores: threat score (TS), probability of detection (POD) and false alarm ratio (FAR)-see Appendix B.

3. Results and Discussion

The results presented herein are divided in two sections. Firstly, the model verification of the mean quantities is explored, stratifying results by height levels. Secondly, the turbulent quantities are evaluated, regarding their temporal evolution, data distributions, daily cycle representations and turbulence spectrum description.

3.1. Model Performance: Mean Quantities

Figure 2 shows the distribution of horizontal wind speed and wind direction through wind roses at three different heights, given by the tower observations (OBS), the mesoscale simulation of 1-km horizontal grid (MESO) and the LES simulation of 111-m horizontal grid (LES). Horizontal wind speed below 0.5 m s${}^{-1}$ are considered as calm winds and not plotted, as the wind direction cannot be determined reliably at these speeds. The dominant wind component from measurements is NW, which is also reproduced in the simulations, increasingly as the height above ground increases. According to previous studies from BLLAST campaign (e.g., [61,62]), northern and western components are mostly associated with large scale circulation, such as frontal passages and strong synoptic flows. The observed light winds from the N and NE directions corresponding to daytime mountain breezes [63] are misrepresented in the MESO and LES simulations. Instead, the S and SE nocturnal flows are well captured by the LES simulation, although simulations tend to exaggerate the frequency of SE flows, specially LES, and to shift slightly the wind towards the east.

In Figure 3 time series of horizontal wind speed and direction observations at $z_{30}$ are compared with MESO and LES simulations. Model simulations are able to reproduce the wind speed evolution tendencies (Figure 3a), slightly overestimating the values, as indicated by the Mean Bias (MB) which is positive in all levels and simulations (

Table 3. Statistics for horizontal wind speed.

Level	z${}_{30}$		z${}_{45}$		z${}_{60}$
Statistic [units]	MESO	LES	MESO	LES	MESO	LES
MB [m s${}^{-1}$]	0.92	0.64	0.43	0.29	0.54	0.09
RMSE [m s${}^{-1}$]	1.35	1.25	1.10	1.27	1.45	1.52
R	0.52	0.45	0.61	0.46	0.66	0.58

). Simulations show relatively low values of RMSE, between 1.25 and 1.52 m s${}^{-1}$ on the three levels (see Table 3), meaning that the wind speed values are relatively well estimated. The linear regression coefficient (R) ranges between 0.45 and 0.66 in different simulations and levels, and is slightly better for MESO than for LES simulations. Moreover the linear regression improves as the height above ground increases, so that z${}_{60}$ gives the best results. Contingency tables and verification scores for different wind speed thresholds at z${}_{30}$ confirm that MESO simulation gives slightly better results than LES (Figure 4). This is illustrated by the fact that the probability of detection (POD) is greater than the false alarm ratio (FAR) for a wider wind speed range in MESO than in LES simulations. The wind speed threshold at which POD is over FAR is around 5.5 m s${}^{-1}$ in MESO and 4.25 m s${}^{-1}$ in LES simulations. These wind speed thresholds increase in height for z${}_{45}$ and z${}_{60}$ although the difference between MESO and LES is maintained (not shown). Threat scores (TS) are very similar between MESO and LES simulations, giving a comparable fraction of observed and forecasted events that were correctly predicted.

The wind direction evolution time series plot shows that MESO and LES capture the main changes in wind direction, which were shown in the general statistics by the wind roses. In particular, simulations reproduce south-eastern dominant winds during the hot period (from 25 to 27 June) defined by [64] (Figure 3b) or the wind shift during intensive observation periods (IOP) days (19–20 June, 24–25 June, 1–2 July, 2–3 July) from north (daytime) to south (night-time) components. Statistics for wind direction show negative mean bias meaning that, in average, modelled wind direction is shifted counterclockwise, greater for LES than for MESO simulations, similar at all levels (

Table 4. Statistics for wind direction.

Level	z${}_{30}$		z${}_{45}$		z${}_{60}$
Statistic [units]	MESO	LES	MESO	LES	MESO	LES
MB [deg]	−11.16	−16.24	−5.94	−9.65	−13.01	−16.60
RMSE [deg]	35.37	51.37	35.30	48.24	35.08	43.62
R	0.86	0.63	0.88	0.67	0.96	0.76

). The square root of the average of squared errors shows a greater difference in LES modelled wind direction than in MESO in all levels, with magnitudes between 35 and 51 degrees. Therefore, MESO simulation, at three levels, performs slightly better the wind direction than LES. Linear regression coefficient for wind direction, calculated using directional statistics following [65,66], shows a very good regression between MESO and OBS (increasing from 0.86 at $z_{30}$ to 0.96 at $z_{60}$). The R is slightly worse for LES.

In general, the wind statistics show similar behaviour between MESO and LES simulations (Table 3 and Table 4), slightly better for MESO, meaning that the mesoscale forcing, driven by the reanalysis forcing, is the responsible for the changes in wind speed and wind direction. The LES simulation does not improve the wind field, instead, it shows higher variability in wind speed and wind direction values. As the wind time series in LES are calculated using 10 min temporal averages from 0.25 s outputs, they include more fluctuations and variability than MESO time series, therefore, statistics can easily be worse for LES. In fact this could be somewhat expected due to the ’double penalty’ effect [67,68], caused by miss-location (in this case temporal shifts) of forecast features found when comparing higher with lower resolution forecasts with deterministic verification scores. The added value of LES is, then, the ability to resolve turbulent quantities (see next Section 3.2).

3.2. Model Performance: Turbulent Quantities

One of the most important magnitudes in the ABL is the turbulent kinetic energy (TKE), which defines the status of the flow and the energy available in the atmosphere for mixing, i.e., the turbulence intensity. The temporal evolution of the TKE at z${}_{30}$ is shown in Figure 5. While in MESO simulation the TKE is obtained from the parametrisation of the PBL scheme, in LES this magnitude is explicitly calculated (see Section 2.3). For the whole period, the MESO and LES simulations reproduce the daily cycle of the TKE and reproduce the turbulence intensity changes through time. Although the correlation is very good at all levels (

Table 5. Statistics for turbulent quantities. Note that the regression coefficient is dimensionless.

Magnitude [Units]	TKE [m${}^2$ s${}^{-2}$]						F${}_{\mathit{H}}$ [K m s${}^{-1}$]
Level	z${}_{\mathbf{30}}$		z${}_{\mathbf{45}}$		z${}_{\mathbf{60}}$		z${}_{\mathbf{30}}$	z${}_{\mathbf{45}}$	z${}_{\mathbf{60}}$
	MESO	LES	MESO	LES	MESO	LES	LES	LES	LES
MB	−0.11	−0.16	−0.09	−0−17	−0.10	−0.15	0	0	−0.01
RMSE	0.24	0.33	0.24	0.32	0.25	0.33	0.03	0.04	0.05
R	0.97	0.45	0.91	0.49	0.93	0.48	0.39	0.17	0.15

), MESO simulation, as expected, is not able to capture the observed turbulence peaks, while LES approximates better the enhanced TKE periods (18, 22, 28 June and 3 of July) but generally underestimating the turbulence intensity. These days correspond to non-IOP periods, when large scale winds from W and N dominated the region.

An examination of the data distributions (probability distribution functions, or PDF) from the whole period LES simulation and observations allows a global comparison of the behaviours and discrepancies between the two data sets. For example, LES simulation generally underestimates TKE values (Figure 6a), either during night-time or during daytime, and extreme TKE values are not well captured. This is also shown by the negative MB at all levels (Table 5), while R coefficients are maintained around 0.45–0.49. The horizontal wind standard deviation ranges are underestimated (Figure 6b), but most significantly, the vertical velocity standard deviation (Figure 6c). While the model hardly reaches values of ${\sigma }_w$ greater than 0.2 m s${}^{-1}$ (more than 90% of values are below 0.2 m s${}^{-1}$), more than the 70% of observations range between 0.2 to 0.6 m s${}^{-1}$ (Figure 6c). Very low values of the kinematic heat flux ($F_H$) are captured by LES (Figure 6d). Positive peaks are often underestimated, corresponding to central hours daytime maxima during IOP days. Linear regression coefficient is relatively low, decreasing as the height increases (Table 5). These results are only displayed for z${}_{30}$, although the other two heights (z${}_{45}$ and z${}_{60}$) show very similar PDF distributions for observations and simulations. Statistics for these heights are shown in Table 5.

A closer look to the daily cycle of the turbulent quantities confirms the general underestimation of TKE and heat flux values, calculated for 10 IOP days (Figure 7). The TKE daily cycle (Figure 7a) is generally underestimated by LES, a feature also observed by [12], but opposite to the results from [6]. These authors and other recent existing literature [22,35,37] use a greater number of vertical levels which can be one reason for these discrepancies. In addition, the sub-grid scheme may not be providing the fully correct energy contained in the eddies smaller than the grid, which can be another source of errors. While morning median and peaks are too low in the simulation, afternoon median gets closer to observations, which means that turbulence in the simulation probably needs more time to get fully developed. Night-time TKE values are maintained in the simulation, although slightly below the observed ones.

The daily cycle seen in heat flux observations is misrepresented by the LES simulation (Figure 7b). During night-time the simulation maintains the heat flux close to zero, mostly positive, whereas in observations the median heat flux is negative, with values around −0.02 K m s${}^{-1}$. During daytime the magnitude of F${}_H$ is noticeably underestimated, and extreme positive values are never reached, a similar behaviour than observed by [2]. Greater values of F${}_H$ are given by LES during IOP daytime central hours in comparison to non-IOP days. Therefore, transport due to thermal origin of turbulence is reproduced by LES to some degree. Instead, no remarkable differences are noticed in the comparison of TKE for the non-IOP days.

The general inaccuracy of the model in reproducing the heat flux can be attributed to different explanations. Firstly, we have seen there is a general underestimation of the standard deviation of vertical velocity component (Figure 6c), which can lead to an underestimation of the heat flux, and the TKE as well. In addition, a great part of the heat flux may be unresolved given the grid resolution. Thus, the sub-grid part of the heat flux becomes significant, which is not resolved by the simulation, as Figure 7b only presents the resolved heat flux. Indeed, in our 111-m horizontal grid size simulations eddies that trigger turbulence may be smaller than the effective resolution, for instance in stably stratified ABL (e.g., [36]). Then, the sub-grid scheme of the model can play an important role, as most part of the heat flux is unresolved and has to be parameterised and it still exists a lot of uncertainty to represent adequately the unresolved sub-grid turbulence. In addition, the grid geometry lead to high anisotropy, breaking the assumption of the LES simulation. On the other hand, the temperature perturbation method used to couple the mesoscale to LES domains can be introducing artificial unrealistic potential temperature fluctuations.

3.3. Velocity Field Spectrum

The power spectral density of the turbulent velocity field, or turbulence spectrum, gives information about the energy associated with each time scale. In this case we can explore the energy associated to limited time scales, between 25 days and 10 min. Figure 8 shows that the energy of eddies with lifetime below 2 or 3 h are maintained in LES simulation, similarly as in observed wind speed spectrum (OBS), but decays too fast in MESO simulation. Therefore, MESO simulation is not able to resolve sufficient energy for frequencies higher than 10${}^{-4}$ (i.e., lower than 2 h), producing a strong decay of the energy. Instead, LES is reproducing the eddies lifetime with similar energy than in observations. The lack of energy identified for frequencies of a few hours in LES can be attributed to the unresolved flow field in the “Terra incognita” scales, probably not fully resolved in the intermediate LES domain.

On the left part of the turbulence spectrum, the energy in the smallest range of frequencies (larger scales) is slightly overestimated in MESO an LES simulations, in comparison with OBS time series, a limitation probably dragged from the reanalysis. In addition, large scales, from several days to few hours are reproduced similarly in both MESO and LES simulations.

4. Conclusions

This study demonstrates the ability of the WRF-LES modelling system to simulate the mean and turbulent quantities of the flow field during the 25 days when BLLAST campaign took place. The long WRF-LES simulation and the field campaign observational data sets provide a unique opportunity to identify the benefits and shortcomings of the mesoscale (MESO) and LES simulations for such a long time period.

The large scale winds are well represented by simulations, and the enhanced turbulence associated with the strong winds as well. Mesoscale flows from north during daytime veering to south during night-time are quite well reproduced by both MESO and LES simulations. The general statistics for wind speed and direction are slightly better for MESO than for LES, although the differences are small and statistics can harm LES because of the greater fluctuations of the time series. Few differences in wind speed and wind direction are seen between levels, better for the highest level, z${}_{60}$, to some degree.

The daily cycle of the turbulent kinetic energy (TKE) is reproduced by LES but values are generally underestimated, mostly during daytime morning and central hours. Non-IOP periods with enhanced turbulence are reproduced by the simulation, but still lowering the intensity. Although the horizontal wind standard deviations are lower in LES than in OBS, the simulation mostly fails in reproducing too low vertical velocity standard deviation. In addition, LES underestimates the daytime heat flux peaks and is not able to reproduce the mostly negative heat flux at night. As shown by [36], simulations at smaller grid sizes would probably resolve a larger part of turbulence, helping to better approximate the turbulent processes. However, there is still a lot of uncertainty in the role of the sub-grid scheme, when the unresolved scales are significant. Despite the test performed with a higher number of levels in the boundary layer indicated similar results, future work should be devoted to test further model configurations with different physical parametrisation schemes and higher number of levels.

In terms of the energy associated with the wind speed resolved in the different time scales, it is shown that LES is able to reproduce sufficient energy for all the turbulence spectrum frequencies, while MESO is limited to time scales larger than a few hours. For higher frequencies, MESO produces a too strong decay of energy where eddies of lifetime smaller than few hours are not reproduced. Therefore, as already documented in previous research, this study also confirms the ability of LES simulation to represent the small scale motions of lifetime around few minutes.

This work provides an insight in the use of high resolution numerical modelling using WRF-LES. WRF-MESO still provides slightly better results for horizontal wind speed and wind direction. The major shortcomings of WRF-LES seem to be related to the grid resolution, which establishes the limit for the scale of resolved turbulence. Another key factor seems to be the underestimation of the unresolved turbulent quantities, given by the sub-grid scheme, that may be important at our LES grid size of 111-m. Further research would include testing for real simulations including sensitivity studies with different grid resolution, increasing vertical levels and exploring the role of the sub-grid scheme. In addition, as pointed out by [2], a fully 3D PBL turbulence parameterisation is probably needed to correctly resolve the flow field mean and turbulent quantities in intermediate domains in the range of the “Terra incognita” [69], to provide better inflow conditions for nested LES domains.