# Full-Scale Experimental Validation of Large-Eddy Simulation of Wind Flows over Complex Terrain: The Bolund Hill

Authors: Ashvinkumar Chaudhari; Antti Hellsten; Jari Hämäläinen
Journal: Advances in Meteorology
 (2016)

---

## Abstract

Numerical simulation of local atmospheric flows around a complex topography is important for many wind energy applications, as it can help in locating and optimizing wind farms. The aim of this paper is to simulate the atmospheric flows over the challenging and complex site of Bolund Hill (in Denmark) using Large-Eddy Simulations (LES) and to validate the simulation methodology against full-scale measurements. Wind prediction over a potential inland wind park is the main application of the validated LES methodology. For the first time, LES-based results from more than one wind direction of the Bolund case are reported and analyzed in detail herein. The LES results from each direction are compared with the Bolund field data in a quantitative way, for example, by using simulation error. According to the results comparison, it is found that the LES model provides by far the most accurate prediction for turbulent kinetic energy. The simulation error of this LES model for predicting turbulent kinetic energy is 19% smaller than that of all other models, whether experimental or numerical, employed previously over Bolund Hill.


---

## Body

## 1. Introduction

The production of wind energy is of ever-increasing interest due to risks of global warming and nuclear energy production plant accidents. Nowadays, many wind farms are being constructed in areas with complex topographies containing hills, ridges, forests, and mountains, with the intention of reaching higher wind speed. Despite the advantages, establishing a turbine over a hill presents specific difficulties, such as a considerable surge in structural loads on wind turbine blades due to inherent flow variability, large-scale unsteadiness, and flow separation, especially in complex terrains. In such terrains, the evaluation of wind resources becomes more and more challenging, and the prediction of wind-farm profitability is directly affected. Therefore, an accurate modelling of atmospheric flows over a complicated terrain has been the subject of considerable research over the past few decades. A field measurement campaign is the most accurate way to carry out wind measurements since it accounts for all environmental parameters and weather fluctuations. However, field measurements are limited to only point locations that may not provide information over the entire site. This is particularly true in the case of complex terrains. Therefore, other alternatives should be taken into account for the scale of a large wind farm. To fulfill this requirement, several numerical and experimental models are being developed to predict the wind on a larger scale. In this way, models can provide valuable information for locating and optimizing wind farms in the absence of field data. Also, models can be employed to simulate the presence of wind turbines and their interaction with the flow. However, all these modelling approaches need to be validated against full-scale measurements.Among the numerical models, the Reynolds-averaged Navier-Stokes (RANS) approach with two-equation turbulence models has been widely used for wind prediction over complex terrains [1–9]. However, it must be mentioned that the RANS models provide solutions on the mean wind speed and the mean level of turbulent kinetic energy due to the Reynolds-averaged approach [10]. Additionally, the models are dependent on a number of constants which are typically calibrated for simple, canonical flows, such as turbulent channel flow [11]. Therefore, the RANS approach still has some fundamental limitations when used to predict flows with complex phenomena. For example, as pointed out by Sumner et al. [12], the RANS formulation is inherently incapable of capturing unsteady effects, such as intermittent flow separation and eddy generation and transport, and generally gives difficulties while modelling turbulence in areas of strong separation [12]. On the other hand, the Large-Eddy Simulation (LES) approach resolves the large (and thus the most important) turbulent structures, and it models only the small-scale motions which are more universal than the large eddies. As a result, LES is able to capture the important unsteadiness in flows. In contrast to RANS, LES typically requires a much higher spatial resolution, resulting in substantially longer computational time [11]. However, due to rapidly increasing computer technology, the use of LES towards atmospheric flow simulation is becoming increasingly common nowadays. For wind energy applications, LES is a promising tool, as it can adequately resolve the most important eddies of the flow which are relevant to wake evolution and power production [13].The aim of this work is to simulate atmospheric flows over a natural complex terrain and to validate the LES methodology against real measurements. Wind prediction over a potential inland wind park is the main application of the validated LES methodology. However, unlike idealized terrains [14–17], the application of LES over realistic complex topography is not so common due to several numerical challenges, such as the realistic upstream boundary condition, heterogeneous nature of surface roughness, resolving near-wall turbulence structures, and highly turbulent wind flow. Particularly, reproducing the inflow boundary condition is one of the critical challenges, as it is generally not known before a numerical simulation. Further, the inflow condition in LES must contain turbulent fluctuations as a function of space and time, with a realistic energy distribution over the spatial directions and the simulated wave-number range. In this study, the so-called recycling technique is employed for generating the instantaneous inflow boundary condition. Recently, the recycling technique was systematically studied and validated against the field measurements [18].Several field experiments over natural terrains have been performed to provide full-scale measurements for validating flow modelling approaches. The field experiments have been carried out over Black Mountain [19], Ailsa Craig Island [20], Blashaval Hill [21], Askervein Hill [22], Kettles Hill [23], Hjardemal Escarpment [24], and so forth. Up to now, the Askervein Hill project [22] has been the most well-known field campaign for validating numerical codes and evaluating turbulence models. In the past, several studies [3, 10, 25–28] reported successful validations of their numerical models chiefly over Askervein Hill. However, it must be mentioned that the steepness of Askervein Hill is below 20∘, and thus it may be characterized as a smooth terrain which presents almost two-dimensional flow features [6]. Also, due to the gentle slope (30∘), the Hjardemal site cannot truly be considered challenging, as the complex flow phenomena, such as flow separation, do not occur.Recently, the Bolund field campaign [29–31] over Bolund Hill in Denmark has provided a new experimental dataset of the mean flow and the turbulence properties. In this work, to validate the LES methodology over a realistic complex terrain, we chose the challenging site of the Bolund. The hill is relatively smaller (12 m high) in scale than Askervein Hill, but it offers a challenging topography with steep slopes and a cliff (90∘). The selection of Bolund Hill is justified by the three-dimensional nature of the hill and by the flow separations induced by the cliff that make it more challenging to model than gentler terrains, such as Askervein Hill or Hjardemal Hill.The first ever numerical and/or experimental (at laboratory scale) validation of the Bolund flows was initiated by Bechmann et al. [30] as a blind comparison of different microscale atmospheric flow models. They employed 57 models ranging from numerical to physical, including LES, RANS, and linearized models, as well as experimental models such as wind tunnel and water channel. From the blind comparison [30], they concluded that a RANS model with two-equation closure appeared to be the best model, due to the smallest simulation error for predicting the velocity speedup and turbulent kinetic energy (TKE), even compared to LES and experimental models. Such interesting feedback from the blind test is obviously encouraging for further investigation of the flow structures over the complicated Bolund site, especially by using LES model. After the blind test, there have been a number of attempts to simulate the Bolund flows by using numerical [6, 32–35] as well as experimental (wind-tunnel) [34, 36–38] models. However, it should be emphasized that all of the previous attempts, except the RANS simulations conducted by Prospathopoulos et al. [6], are limited to only one-wind-direction results. To our knowledge, the present work is the first LES-based study presenting and analyzing the results from more than one inflow wind direction of the Bolund case. In addition, the previous studies have focused mostly on the mean flow analysis but have not focused on the turbulence statistics in enough detail. In fact, the information on turbulence prediction is one of the important parts for a complete model validation.In the present work, LES are carried out to investigate the atmospheric flow over the Bolund Hill, and the LES model is validated against the full scale measurement [29, 30]. The LES-based results from two different wind directions, 270∘ and 239∘, are reported and analyzed in detail for both the mean and the turbulence quantities. Following Bechmann et al. [30], the results are evaluated in quantitative manner, for example, by using simulation error. The results from the two different wind directions are presented together and discussed in order to learn the differences and/or similarities that emerge owing to a change in wind direction. In the final section, using the simulation error and the model-validation criteria, this paper evaluates and compares the performance of the LES model with other models previously employed over the same hill. Through the quantitative comparison, it is found that the present LES yields the most accurate prediction for the TKE, due to the smallest simulation error. The simulation error for the TKE increase is 19% smaller than in any other (experimental or numerical) models employed so far over the Bolund Hill.


## 2. The Bolund Experiment

Bolund Hill is a 12 m high, roughly 130 m long, and 75 m wide, and it is located on the north of DTU (Technical University of Denmark) in the city of Roskilde (Denmark). It is surrounded by sea, with the exception of the east side, which is connected to the mainland by a narrow isthmus. The field experiment was performed over a period of three months in 2007-2008, with the intention of providing a new full-scale dataset for validating atmospheric flow models in wind energy. Figure1 shows the Bolund orography data with the locations of measuring masts (M0 to M9) installed on the hill during the field campaign. In this field campaign, 23 sonic and 12 cup anemometers were simultaneously employed to measure the wind velocity and variances from the 10 different masts. As it can be seen in Figure 1, the measuring masts (referred to hereafter as M1 to M8) were placed on two lines, Line-A and Line-B. In addition, two more masts, M0 and M9, were utilized to provide the inflow wind conditions for westerly and easterly winds, respectively. At each mast, the data were recorded for four different wind directions: 270∘, 255∘, 239∘, and 90∘. According to Berg et al. [29], the mean wind speed during the measurements was around 10 m/s, which leads to a Reynolds number Re of 107, and, thus, the flow can be considered to be independent of the Re. Also, the field data were collected only for the neutral stability condition. In addition to offering the latest full-scale database, the Bolund field experiment also provides an interesting analysis from the blind comparison performed by Bechmann et al. [30], which allows one to analyze the validation in great detail. More detailed information about the field experiment can be found in [29–31].Figure 1
Bolund contour plot colored to illustrate the height. The black dots show the actual positions of masts (M0 to M9). The wind direction270∘ is from west to east.


## 3. LES Methodology

### 3.1. LES Model and the Code

In LES, the turbulent eddies larger than the grid size, which are the most energetic ones and responsible for the majority of turbulent transport, are resolved directly, whereas the eddies smaller than the grid size are more isotropic and are modelled using a subgrid scale (SGS) model. The governing equations for LES can be obtained by filtering the original continuity and Navier-Stokes equations. The filtered governing equations for incompressible flow can be written as(1)∂u¯i∂xi=0,∂∂tu¯i+∂∂xju¯iu¯j=ν∂∂xj∂u¯i∂xj-1ρ∂p¯∂xi-∂τij∂xj,where u¯i and p¯ are the grid-filtered values of the velocity and the pressure, respectively, ν is the kinematic viscosity, ρ is the fluid density, and τij is the subgrid scale stress defined as follows:(2)τij=uiuj--u-iu-j.Since the smaller eddies are not resolved in the simulations, their effects are modelled using the one-equation eddy viscosity SGS model proposed by Yoshizawa [39]:(3)∂kSGS∂t+∂∂xju¯jkSGS=∂∂xjν+νSGS∂kSGS∂xj+PkSGS-CϵkSGS3/2Δ,where PkSGS=2νSGSS¯ijS¯ij, S¯ij is the strain rate tensor, and νSGS=CkkSGS1/2Δ, Δ=V1/3, where V is the volume of the computational cell. Ck=0.094 and Cϵ=1.048.In this work, the unstructured finite-volume method based code, OpenFOAM [40], is employed to solve the governing equations. The LES calculations are carried out by using our own in-house OpenFOAM-based LES solver [18, 33], recently implemented into OpenFOAM. The solver uses the classical fourth-order Runge–Kutta (RK4) scheme for time integration and a projection method where the velocity field is corrected using the pressure gradient in-between the RK4 substeps. The pressure is then obtained by solving the Poisson equation. More information on the LES solver, including its validation, can be found in [18]. The one-equation model (3) is time integrated numerically using the second-order backward implicit method, and it is solved after the fourth (last) RK4 substep. All the governing equations are discretized in space using the central differencing scheme. Following the recommendations of Berg et al. [29] and Bechmann et al. [30], the present simulations do not include the effects of the Coriolis force, and the atmosphere was considered to be neutrally stratified.


### 3.2. Numerical Domain and Grid Resolution

Figure2 shows the computational setup of the simulations, including an explanation of the recycling approach, and Figure 3 shows the numerical model of the hill with a surface grid. As discussed previously, the simulations use the recycling technique to generate the time-dependent inlet boundary conditions at the inflow plane (see Figure 2). For this reason, the length of the computational domain must be extended on the upstream side in order to ensure the fully developed flow before the reference mast M0. In the stream-wise (x) direction, the domain ranges from x=-810 m to 400 m. The location of the outlet boundary was followed from the previous RANS study [6]. The vertical height (z) of the domain is set to 120 m, and the boundary-layer depth δ is also 120 m, which is approximately 10 times higher than the height of Bolund Hill H, that is, δ≈10H. The computational domain is also extended from −205 m to 185 m in the span-wise (y) direction to avoid the periodic effects on the hill part, as the periodic boundary conditions are applied in that direction. The center of the hill is located at x=y=0. In this way, the LES computational domain is of the size 1210×390×120m3 in x, y, and z directions, respectively.Figure 2
Picture of the computational domain with an explanation of the recycling technique used in the simulation of the270∘ wind direction. A similar setting is used for the 239∘ wind direction, but with the x-axis aligned with that direction.Figure 3
Numerical model of Bolund Hill with a surface grid.The computational domains are then discretized using a block-structured hexahedron mesh with finer (order of cm) grid resolution on the hill surface. To select the grid resolutions in the present case, some guidelines were followed from the previous work [6, 35] dealing with the same hill. The computational domains were meshed using the commercial grid generator called Pointwise [41]. In both LES cases, the grid created was nonisotropic with the finest grid spacing of Δx=0.31m and Δy=0.40m being kept on the vertical escarpment (see Figure 3). On the rest of the hill part, the stream-wise resolution Δx varies from 0.31 m to 0.85 m and the span-wise Δy from 0.40 m to 0.70 m depending on the local steepness of the hill. Outside the hill part, grid sizes in both horizontal directions are relaxed up to 4 m near the domain boundaries. The vertical grid spacing Δz also varies gradually from 0.15 m to 4.97 m at the top boundary. In this way, the entire computational grid consists of 940×428×100(≈40×106) hexahedron cells in each LES calculation.


### 3.3. Boundary Conditions

Reproducing the correct inflow velocity conditions is a challenge in LES. The inflow conditions must contain a realistic distribution of turbulence level throughout the boundary layer. Previously in several LES calculations [10, 25–27, 35], a separate precursor simulation over a flat and homogeneous surface has been done to generate inlet conditions for complex terrain simulations. This means that the simulation must be done in two parts (for two domains or grids), and then the data from the first (precursor) simulation is used in the second (main) simulation, that is, over complex terrain. As a consequence, the entire process (precursor and main simulations) takes a longer clock time for finalizing the simulation. Further, it requires large memory storage to store the data at each time step. Even more importantly, if one wants to make any changes in the main computational domain or grid, then the precursor simulation needs to be run again after introducing the new parameters. However, on the positive side, the precursor method may be the best option when several simulations are planned to be run using the same inflow conditions, as one can utilize the stored data in all those simulations. For the present case, the fully developed inflow boundary condition is obtained using the recycling (or mapping) method; see, for example, Chaudhari [18]. With the recycling technique, the precursor simulation is not required, and thus the entire simulation is performed in a single step; however, the user have to bear the extra CPU cost of the recycling section every time.Using the recycling technique, the flow variables are sampled on a span-wise plane (recycling plane in Figure2), which is sufficiently downstream of the inflow plane. Then, the sampled data are recycled back to the inflow plane at each time step, as explained in Figure 2. In this way, a recycling section is created between the inflow and the recycling planes, in which the flow becomes fully developed. Recently, Chaudhari [18] tested the suitability of the method under different recycling lengths Lrcy, that is, the distance between the inflow and the recycling planes. It was concluded that Lrcy should be at least three times the boundary-layer depth, that is, Lrcy≥3δ, in order to avoid artificial turbulence structures within the recycling section [18]. Hence, in this case, Lrcy=480m=4δ=40H was used. Furthermore, the average velocity flux is fixed at the inflow boundary in order to maintain the same amount of volumetric flow throughout the simulation. The static pressure is fixed to a constant value on the outlet plane and the zero-gradient (Neumann) [40] boundary condition is used for the rest of the flow variables. The slip boundary condition is used for all the flow variables at the top boundary. The periodic boundary conditions are set in the cross-wind direction. A roughness-length z0-based wall function is utilized to determine the ground surface fluxes on the lower boundary. The implementation and the validation of the wall function can be found from [18, 33]. Following Berg et al. [29], two types of the roughness-length values, z0=0.0003m and 0.015m, are used in the simulations for the water and the hill surfaces, respectively.In this work, two LES calculations for the two inflow wind directions (270∘ and 239∘) are carried out. The simulations are run until the physical time of t=2700 s, in which all the flow quantities are time-averaged over the last 600 s for calculating the flow statistics. The wall-clock time for each of the calculations was about 40 days using 512 processors for simulating the physical time of 2700 s. The calculations were carried out using Finland’s most powerful supercomputer called Sisu (Cray XC30). Due to the high CPU requirement, a grid independent study of the present case is subject to large computational resources. We therefore rely on the agreement of the results (see Table 1) between the present LES and the previous numerical and experimental cases in order to evaluate the quality and the reliability of the present simulations. For the same reason (high CPU resources), the current study is based on the 10 min time-averaged results. However, the averaging time taken here is more than 6.5 advection times along the entire computational domain, or 32 advection times along the hill length, and importantly it is 3 times as long as that used in the previous LES by Diebold et al. [35].Table 1
Mean absolute errors of speedup and TKE increase from the literature and in the present LES. The values in brackets show the best performance of the particular model type.
Model Mean error ofΔS (best) Mean error ofΔk (best) Wind direction(s) Reference RANS 2 eq. 15.1 (11.4) 47.0 (29.9) 27 0 °  and  239° Bechmann et al. [30] RANS 1 eq. 17.2 (13.8) 44.7 (42.7) 27 0 °  and  239° Bechmann et al. [30] Experiment 14.7 (13.3) 61.4 (59.4) 27 0 °  and  239° Bechmann et al. [30] LES 17.3 (14.1) 48.0 (41.6) 27 0 °  and  239° Bechmann et al. [30] Linearized 23.7 (20.6) 76.7 (71.4) 27 0 °  and  239° Bechmann et al. [30] RANS 2 eq. 10.3 — Only270° Prospathopoulos et al. [6] Wind-tunnel 12.4 42.2 Only270° Conan [36] LES-EPFL 9.6 — 27 0 °  &  239° Diebold et al. [35] Wind tunnel 17 56.3 Only270° Yeow et al. [37, 38] Present LES 10.3 24.1 27 0 °  and  239° —


## 3.1. LES Model and the Code

In LES, the turbulent eddies larger than the grid size, which are the most energetic ones and responsible for the majority of turbulent transport, are resolved directly, whereas the eddies smaller than the grid size are more isotropic and are modelled using a subgrid scale (SGS) model. The governing equations for LES can be obtained by filtering the original continuity and Navier-Stokes equations. The filtered governing equations for incompressible flow can be written as(1)∂u¯i∂xi=0,∂∂tu¯i+∂∂xju¯iu¯j=ν∂∂xj∂u¯i∂xj-1ρ∂p¯∂xi-∂τij∂xj,where u¯i and p¯ are the grid-filtered values of the velocity and the pressure, respectively, ν is the kinematic viscosity, ρ is the fluid density, and τij is the subgrid scale stress defined as follows:(2)τij=uiuj--u-iu-j.Since the smaller eddies are not resolved in the simulations, their effects are modelled using the one-equation eddy viscosity SGS model proposed by Yoshizawa [39]:(3)∂kSGS∂t+∂∂xju¯jkSGS=∂∂xjν+νSGS∂kSGS∂xj+PkSGS-CϵkSGS3/2Δ,where PkSGS=2νSGSS¯ijS¯ij, S¯ij is the strain rate tensor, and νSGS=CkkSGS1/2Δ, Δ=V1/3, where V is the volume of the computational cell. Ck=0.094 and Cϵ=1.048.In this work, the unstructured finite-volume method based code, OpenFOAM [40], is employed to solve the governing equations. The LES calculations are carried out by using our own in-house OpenFOAM-based LES solver [18, 33], recently implemented into OpenFOAM. The solver uses the classical fourth-order Runge–Kutta (RK4) scheme for time integration and a projection method where the velocity field is corrected using the pressure gradient in-between the RK4 substeps. The pressure is then obtained by solving the Poisson equation. More information on the LES solver, including its validation, can be found in [18]. The one-equation model (3) is time integrated numerically using the second-order backward implicit method, and it is solved after the fourth (last) RK4 substep. All the governing equations are discretized in space using the central differencing scheme. Following the recommendations of Berg et al. [29] and Bechmann et al. [30], the present simulations do not include the effects of the Coriolis force, and the atmosphere was considered to be neutrally stratified.


## 3.2. Numerical Domain and Grid Resolution

Figure2 shows the computational setup of the simulations, including an explanation of the recycling approach, and Figure 3 shows the numerical model of the hill with a surface grid. As discussed previously, the simulations use the recycling technique to generate the time-dependent inlet boundary conditions at the inflow plane (see Figure 2). For this reason, the length of the computational domain must be extended on the upstream side in order to ensure the fully developed flow before the reference mast M0. In the stream-wise (x) direction, the domain ranges from x=-810 m to 400 m. The location of the outlet boundary was followed from the previous RANS study [6]. The vertical height (z) of the domain is set to 120 m, and the boundary-layer depth δ is also 120 m, which is approximately 10 times higher than the height of Bolund Hill H, that is, δ≈10H. The computational domain is also extended from −205 m to 185 m in the span-wise (y) direction to avoid the periodic effects on the hill part, as the periodic boundary conditions are applied in that direction. The center of the hill is located at x=y=0. In this way, the LES computational domain is of the size 1210×390×120m3 in x, y, and z directions, respectively.Figure 2
Picture of the computational domain with an explanation of the recycling technique used in the simulation of the270∘ wind direction. A similar setting is used for the 239∘ wind direction, but with the x-axis aligned with that direction.Figure 3
Numerical model of Bolund Hill with a surface grid.The computational domains are then discretized using a block-structured hexahedron mesh with finer (order of cm) grid resolution on the hill surface. To select the grid resolutions in the present case, some guidelines were followed from the previous work [6, 35] dealing with the same hill. The computational domains were meshed using the commercial grid generator called Pointwise [41]. In both LES cases, the grid created was nonisotropic with the finest grid spacing of Δx=0.31m and Δy=0.40m being kept on the vertical escarpment (see Figure 3). On the rest of the hill part, the stream-wise resolution Δx varies from 0.31 m to 0.85 m and the span-wise Δy from 0.40 m to 0.70 m depending on the local steepness of the hill. Outside the hill part, grid sizes in both horizontal directions are relaxed up to 4 m near the domain boundaries. The vertical grid spacing Δz also varies gradually from 0.15 m to 4.97 m at the top boundary. In this way, the entire computational grid consists of 940×428×100(≈40×106) hexahedron cells in each LES calculation.


## 3.3. Boundary Conditions

Reproducing the correct inflow velocity conditions is a challenge in LES. The inflow conditions must contain a realistic distribution of turbulence level throughout the boundary layer. Previously in several LES calculations [10, 25–27, 35], a separate precursor simulation over a flat and homogeneous surface has been done to generate inlet conditions for complex terrain simulations. This means that the simulation must be done in two parts (for two domains or grids), and then the data from the first (precursor) simulation is used in the second (main) simulation, that is, over complex terrain. As a consequence, the entire process (precursor and main simulations) takes a longer clock time for finalizing the simulation. Further, it requires large memory storage to store the data at each time step. Even more importantly, if one wants to make any changes in the main computational domain or grid, then the precursor simulation needs to be run again after introducing the new parameters. However, on the positive side, the precursor method may be the best option when several simulations are planned to be run using the same inflow conditions, as one can utilize the stored data in all those simulations. For the present case, the fully developed inflow boundary condition is obtained using the recycling (or mapping) method; see, for example, Chaudhari [18]. With the recycling technique, the precursor simulation is not required, and thus the entire simulation is performed in a single step; however, the user have to bear the extra CPU cost of the recycling section every time.Using the recycling technique, the flow variables are sampled on a span-wise plane (recycling plane in Figure2), which is sufficiently downstream of the inflow plane. Then, the sampled data are recycled back to the inflow plane at each time step, as explained in Figure 2. In this way, a recycling section is created between the inflow and the recycling planes, in which the flow becomes fully developed. Recently, Chaudhari [18] tested the suitability of the method under different recycling lengths Lrcy, that is, the distance between the inflow and the recycling planes. It was concluded that Lrcy should be at least three times the boundary-layer depth, that is, Lrcy≥3δ, in order to avoid artificial turbulence structures within the recycling section [18]. Hence, in this case, Lrcy=480m=4δ=40H was used. Furthermore, the average velocity flux is fixed at the inflow boundary in order to maintain the same amount of volumetric flow throughout the simulation. The static pressure is fixed to a constant value on the outlet plane and the zero-gradient (Neumann) [40] boundary condition is used for the rest of the flow variables. The slip boundary condition is used for all the flow variables at the top boundary. The periodic boundary conditions are set in the cross-wind direction. A roughness-length z0-based wall function is utilized to determine the ground surface fluxes on the lower boundary. The implementation and the validation of the wall function can be found from [18, 33]. Following Berg et al. [29], two types of the roughness-length values, z0=0.0003m and 0.015m, are used in the simulations for the water and the hill surfaces, respectively.In this work, two LES calculations for the two inflow wind directions (270∘ and 239∘) are carried out. The simulations are run until the physical time of t=2700 s, in which all the flow quantities are time-averaged over the last 600 s for calculating the flow statistics. The wall-clock time for each of the calculations was about 40 days using 512 processors for simulating the physical time of 2700 s. The calculations were carried out using Finland’s most powerful supercomputer called Sisu (Cray XC30). Due to the high CPU requirement, a grid independent study of the present case is subject to large computational resources. We therefore rely on the agreement of the results (see Table 1) between the present LES and the previous numerical and experimental cases in order to evaluate the quality and the reliability of the present simulations. For the same reason (high CPU resources), the current study is based on the 10 min time-averaged results. However, the averaging time taken here is more than 6.5 advection times along the entire computational domain, or 32 advection times along the hill length, and importantly it is 3 times as long as that used in the previous LES by Diebold et al. [35].Table 1
Mean absolute errors of speedup and TKE increase from the literature and in the present LES. The values in brackets show the best performance of the particular model type.
Model Mean error ofΔS (best) Mean error ofΔk (best) Wind direction(s) Reference RANS 2 eq. 15.1 (11.4) 47.0 (29.9) 27 0 °  and  239° Bechmann et al. [30] RANS 1 eq. 17.2 (13.8) 44.7 (42.7) 27 0 °  and  239° Bechmann et al. [30] Experiment 14.7 (13.3) 61.4 (59.4) 27 0 °  and  239° Bechmann et al. [30] LES 17.3 (14.1) 48.0 (41.6) 27 0 °  and  239° Bechmann et al. [30] Linearized 23.7 (20.6) 76.7 (71.4) 27 0 °  and  239° Bechmann et al. [30] RANS 2 eq. 10.3 — Only270° Prospathopoulos et al. [6] Wind-tunnel 12.4 42.2 Only270° Conan [36] LES-EPFL 9.6 — 27 0 °  &  239° Diebold et al. [35] Wind tunnel 17 56.3 Only270° Yeow et al. [37, 38] Present LES 10.3 24.1 27 0 °  and  239° —


## 4. Results and Discussions

We begin this section by discussing the results at the reference mast M0. After this, the LES results at the rest of the masts (M1 to M8) are presented and discussed. The results from both wind directions are discussed simultaneously in order to illuminate the differences between the results from the two directions. Finally, simulation error is calculated for both wind directions, and the errors are compared against each other as well as against the literature values.


### 4.1. Results from the Reference Location

Figure4 compares the LES profiles of the mean velocity magnitude S and the turbulent kinetic energy k (TKE) with the field measurements at the reference location M0. The results are normalized by the frictional velocity u∗0 at M0. Here, the mean flow profile is very well reproduced by the recycling approach, as the simulated profile agrees well with the measurements at all heights. The TKE is also well predicted by the simulation at 5 m, but it is underestimated by about 25% at the height of 12 m. The total (modelled and resolved) TKE also underestimates the measured TKE at 12 m. The measurements show almost constant level of the TKE at M0, but in the simulation, the TKE decreases with height. Earlier, the LES runs submitted in the blind comparison also did not predict the constant TKE profile. This indicates that the deficiency could be very typical for LES modelling in general and thus occurs in the preset model, too. Since the modelled part of the TKE is significant in the area very close to the surface, we present the total TKE in all the results (Figures 6–13) hereafter. The error bars in the measurements show the standard deviations of the respective mean values. The standard deviation values were provided together with other statistics in the field data [29, 30].Figure 4
The LES profiles of the mean velocity magnitude (a) and the TKE (b) compared with the field measurements for the270∘ wind direction.
(a) (b)


### 4.2. Results on the Hill

Here, we first visualize the LES results on the instantaneous and mean flow structures around the hill. Figure5 depicts the instantaneous velocity fields as well as the 10-min time-averaged ones. It can clearly be seen that the local flow field is strongly affected by the hill. In particular, the cliff is noted to enhance the turbulent fluctuations. The complexity occurring in the flow due to the cliff is also illustrated by the streamlines which highlight the mean (10-min) flow patterns around the hill (Figure 5(b)). As expected, the flow is found to be separated at several places just after the cliff as well as on the lee sides. Particularly, the cliff is forming the small zone of the reversed flow, near the ground and along the entire cliff. Masts M2 and M6 are highly influenced with this reversed flow. However, the previous studies [6, 35] including those wind tunnels [34, 36–38] did not predict this flow separation nearby the cliff.Figure 5
LES-predicted instantaneous stream-wise velocity fields on the vertical plane (a), and the 10 min time-averaged streamlines along the hill (b). The results are from the270∘ wind direction.
(a) (b)Figure 6
Scatter plots of field data against LES results; solid lines: interval of±10% around perfect fitting (dashed line).
(a)
Mean velocity (b)
TKEFigure 7
Ratio of LES to field data for the normalized velocity magnitudeS/u∗0. The results from both wind directions are sorted with respect to the height above the ground level zagl (a) and the mast measurement locations (b).
(a) (b)Figure 8
Ratio of LES to field data for the normalized TKEk/u∗02. The results from both wind directions are sorted with respect to the height above the ground level zagl (a) and the mast measurement locations (b).
(a) (b)Figure 9
Vertical profiles of the velocity speedup (ΔS) at different masts (M1–M8) in the case of 270∘ (a) and 239∘ (b) wind directions.
(a) (b)Figure 10
Vertical profiles of the TKE increase (Δk) at different masts (M1–M8) in the case of 270∘ (a) and 239∘ (b) wind directions.
(a) (b)Figure 11
Simulation errors:RS for speedup and Rk for TKE increase errors in the case of 270∘ (a) and 239∘ (b) wind directions.
(a) (b)Figure 12
Simulation errors:RS for speedup and Rk for TKE increase errors at different measurement locations. The errors (RS and Rk) are averaged over the height for each location.Figure 13
Simulation errorsRS for speedup and Rk for TKE increase at different measurement heights. The errors (RS and Rk) are averaged over all the locations (masts) having the same height.Next, the results comparison between the simulation and the field data at different heights of all 10 masts (total 37-38 locations) is discussed. In order to see the correlation of the results, the scatter plots of all the available field data against the LES predictions for the mean velocity magnitude and the TKE are shown in Figures6(a) and 6(b), respectively. In the plots, an error interval of ±10% is also drawn for reference purpose. In Figure 6(a), it can be observed that 82% of the LES-predicted mean velocity in the case of 270∘ is within the interval of ±10%. The velocity prediction is almost the same, 81%, in the case of 239∘. However, a very small deviation (about 9%) is observed between the two cases when examining the turbulence part. For the TKE, 92% and 83% of the predictions appear within the interval (±10%) in the cases of 270∘ and 239∘, respectively. When one considers the interval of ±5%, the TKE prediction is still good, 80%, in both wind directions.To describe the discrepancies between the simulations and the measurements in detail, more scatter plots on the velocity and TKE predictions are presented in Figures7 and 8, respectively. In these figures, the results are sorted with respect to the measurement height above the ground level zagl (Figures 7(a) and 8(a)) and with respect to the masts (Figures 7(b) and 8(b)). The locations of the masts can be found from Figure 1. According to Figure 7(a), the predicted mean velocity is scattered for the height below 2 m, reaching a maximum error of almost 128% at M2 in the case of 270∘. On the other hand, the largest error (at 2 m) is smaller by about 35% in the case of 239∘. Generally, in both cases at zagl=2, masts M2, M6, and M8 have larger errors than the rest of the masts, as can also be seen in Figure 7(b). The large errors are due to the mast locations in which the flows are highly turbulent with flow separations (see Figure 5(b)). In both cases, the maximum error is reduced by about 100% within a distance of 3 m, from zagl=2m to 5 m. Above the height of 5 m, the prediction is greatly improved, as the maximum error is only 11%, while the majority of the error is below 5%.The overall prediction of the TKE is also scattered in both cases (270∘ and 239∘), although the level of maximum error is reduced by about 90% (Figure 8(a)) as compared to that of the velocity. The largest discrepancy is again found at M2, which is 70% (at 3.6 m) in the case of 270∘ and 65% (at 2 m) in the case of 239∘. However, unlike the velocity, the scatter is relatively independent of the mast locations and measurement heights (Figure 8(b)). Despite having less error (as compared to velocity), the TKE prediction does not seem to be improved with the height from 2 m to 5 m (Figure 8(a)). Overall, the velocity prediction is better in the case of 239∘, whereas the TKE prediction is better in the case of 270∘. More quantitative comparison between the two-wind-direction results is presented in Section 4.3.Following Bechmann et al. [30], the numerical results are further analyzed in terms of velocity speedup ΔS and TKE increase Δk. The speedup and the TKE increase are defined as follows [30]: (4)ΔS=Szagl-S0zaglS0zagl,Δk=kS02zagl-k0S02zagl,where the subscript 0 indicates the values (S0 or k0) from the reference mast M0. Vertical profiles of the speedup and the TKE increase at all the masts are shown in Figures 9 and 10, respectively.Here, it can be observed that the LES prediction of the speedup is generally acceptable (error less than 10%; Figure11(a)) at most (70%) of the places for both wind directions. Nevertheless, for the 270∘ wind, the simulation overestimates the measurements at the lowest heights of M2, M5, and M6, while it underestimates them at zagl=5 of M8. The maximum speedup error at the lowest height (2 m) is 53% and it decreases to 25% at 5 m in the case of 270∘ wind. Concerning the lowest heights of M2, M5, and M6, the speedup is predicted better in the case of 239∘, as the maximum error is only about 34%, which is at 2 m height of M2 (see Figure 11(b)). We note that the vertical profiles of the predicted speedup at M5 and M6 (below 5 m) differ significantly between the two directions. At the rest of the masts, there is no such a significant difference with respect to the wind direction. This means the flows in those regions (M5 and M6) are highly sensitive to the wind direction. According to the measurements, the speedup prediction in those areas is better in the case of 239∘.The large discrepancies between the simulations and the measurements at M2, M5, M6, and M8 are justified by the complex nature of the flows locally, that is, highly turbulent and reversing with several flow separations generated by the steepness of the hill. Particularly, masts M2 and M6 were placed about 10 m away from the cliff and were strongly influenced by the very small flow separation (due to the cliff) along the edge of the cliff, as it can be seen in Figure5(b). This particular separation is indeed much challenging to be even captured by numerical simulations because of its small size and location (near the cliff). Consequently, none of the previous attempts, including the wind-tunnel tests, have predicted this flow separation. Masts M8 and M4 were located on the lee side of the hill for 270∘ and 239∘ winds, respectively, and were also continuously influenced by flow separations. The simulation does predict most of the flow separations there, but the depth of the flow separation is somewhat compromised, although more at M8 than at M4. Previous studies [6, 30, 35–38] also reported large discrepancies at several masts, except for M1 and M7.Due to the flow separations, the turbulence level is significantly increased in a few locations. As a result, a clear peak in the TKE increase profile is visible at M2 and M6 for both winds and at M5 for the270∘ wind (see Figures 10(a) and 10(b)). Importantly, the simulation appears to accurately predict the extreme TKE vertical gradient, especially in the case of 270∘, including the lower heights (Figure 10(a)), even though the velocity speedup is highly overestimated at the same locations (i.e., M2 and M6). With the exception of M2 and M6, the TKE increase is also reproduced satisfactorily in the case of 239∘ direction. At 2 m height of M2 (for 239∘), the simulation shows the largest error of 140% for predicting the TKE increase (Figure 11(b)). Most of the discrepancies between the simulation and the measurements, visible in Figures 9 and 10, and especially at M2, M6, and M8, may also be justified by the uncertainties in the measured statistics. The high standard deviations at the lowest heights of M2 and M6 indicate high variabilities of the measured data at those places.


### 4.3. Simulation Error and Validation Metric

To evaluate the performance of the LES model in a quantitative way, the simulation (modelling or prediction) error from the LES results is calculated by considering the field data as a reference. Originally, Bechmann et al. [30] used the following equations to calculate the simulation errors for speedup and TKE increase. The error of the simulated speedup, RS, is defined as follows:(5)RS=100ΔSs-ΔSm,and the error of the simulated TKE, Rk, is(6)Rk=100Is-I0s-Im-I0mI0m,where ΔS is the speedup, I=k/S0, subscript zero denotes that the value is taken at the reference locations, and subscripts m and s denote the measured and simulated values, respectively. Figure 11 shows the simulation errors for the speedup and the TKE increase at each individual mast as well as at each possible height for the both wind directions. In the same way, Figure 12 displays the simulation errors, but they are averaged over the heights of each mast. The individual error from each height and/or location in Figure 11 may also be used for comparison with the results shown in the previous figures.In the following, we turn our focus on the averaged error to discuss the results. According to Figure12, the speedup error RS is less than 10% for half (50%) of the masts in both cases. The maximum error in the case of 270∘ is 27%, which is found at M6. On the other hand, the maximum error for the 239∘ wind is found at M5 and it is 18%. For the TKE increase, the error (Rk) is much higher than the speedup error (RS), reaching 59% (at M2) and 80% (at M6) in the cases of 270∘ and 239∘ winds, respectively. The high error for the TKE increase is obviously because the second-order statistic, for example, TKE or the Reynolds stress, is challenging to predict with enough accuracy as compared to mean flow (first-order statistic) in general. An interesting fact noted here is that the TKE increase error is almost the same at M2 for both winds, but a significant difference (75%) is observed at M6 due to the change in wind direction. Despite the difficulty of predicting the TKE compared to the mean wind, the Rk error is found to be smaller than RS at M6 and M8 in the case of 270∘ and at M5 in the case of 239∘. The high error at certain masts, such as M2 and M6, is due to their locations (near the cliff), where the aforementioned extremely complex flow behaviors occur. Indeed, virtually all previously reported models, including those wind tunnels, have failed to predict the flow accurately enough at masts M2 and M6. When excluding both M2 and M6, the averaged speedup error changes from 10.3% to 10.2%, which is insignificant. However, a high reduction from 24.1% to 13.5% is observed for the turbulence quantity. In addition, among all the masts, the standard deviation of the measured TKE is also at its highest at M2 and M6, reaching up to 24% in the case of 239∘, as could be expected from the error bars in Figures 9 and 10. The standard deviation or any other uncertainties of the measurements are, of course, not taken into account in the calculations of the simulation error.Next, Figure13 shows the (averaged) simulation error at three different heights. Here, the error is averaged over the measurement locations having the same height. For the 270∘ wind direction, it can be seen that the errors for the speedup and for the TKE increase at the height of 2 m are nearly the same (21%). This indicates that, at the lowest height, the model performs equally for both the speedup and the TKE increase. Above 2 m, the speedup error decreases rapidly down to 2.5% at 9 m. However, this does not happen with the TKE increase. The TKE increase error remains in the same order from 2 m to 5 m, and it decreases nominally, only by 4%, at 9 m.For the239∘ wind direction, the speedup error also decreases with increasing height, but not in the same trend as in the 270∘ case. It reduces from 13.5% at 2 m to 5.5% at 9 m. The TKE error is much higher, about 41% at 2 m (Figure 13). At this height, the TKE error is almost double compared to that of the 270∘ case. The high error is mostly due to the poor predictions at M2 and M6. At 5 m, the TKE error reduces greatly to 21%, about the same level as in the 270∘ case; however, it does not reduce further at 9 m.It should be emphasized that the measurements of the TKE at 2 m and 9 m heights of the reference location (M0) were missing in the field data (see Figure4(b)). Therefore, by following [30], a constant value of the TKE is considered at all three heights (2 m, 5 m, and 9 m) to calculate the measured TKE increase. This assumption (constant TKE with height) can highly bias the errors especially at 2 m, as the LES results follow the height-dependent TKE values from the simulation. Moreover, it should be remembered that the flow within the surface layer, for example, zagl≤2m in the Bolund case, may be of little interest for a wind-farm developer, as one is considering the wind condition within the rotor diameter. With the average over both the wind directions, the speedup error at 5 m is 9.3% and is reduced to 3.9% at 9 m. The error is not completely negligible, but the prediction could be considered accurate (less than 10% error) [30, 35].Table1 shows the mean errors from the two wind directions, and it compares the errors from the present case with those reported previously in the literature. The mean simulation error of the present LES is 10.3% and 24.1% for speedup and TKE increase, respectively. The LES prediction of the speedup is among the best results, as the mean speedup error from our case is only 0.7% higher than that of the best results [35]. Concerning the turbulence aspect, the present results show far smaller error in comparison with all previous results. Indeed, all prior models, except a RANS two-equation model in [30], are suffering a high TKE increase error. A couple of studies [6, 35] did not report the TKE increase error at all. The TKE increase error of 24.1% from this case is about 19% smaller than that of any other model. To our knowledge, this is thus the best prediction of the TKE over the Bolund Hill among all models employed over the hill so far.Finally, the present LES model is evaluated here using three different validation metrics. First, the validation metric called the fraction of two (FAC2) shows the fraction of data that is within the interval:(7)0.5≤PO≤2,where P is the predicted and O is the observed value. For an ideal model prediction, one would have FAC2 = 1 [43]. Second, the Fractional Bias (FB) can be defined as(8)FB=2·O-PO+P,where the angular bracket 〈·〉 represents the average value over the data samples. FB indicates the relative amount of bias in the modelled results such that FB =0 means no bias at all and FB >0 means systematic underestimation [44]. Finally, a third metric, known as the Normalized Mean Square Error (NMSE), is defined as follows:(9)NMSE=O-P2OP.NMSE is sensitive to large disagreements, in particular those occurring on the highest values [44]. The validation-metric values from the present results for both wind directions are listed in Table 2. According to the table, the present LES model passes all three validation tests, using the acceptance criteria for the model evaluation (in a rural area) proposed by Hanna and Chang [42]. In particular, the FAC2 value obtained by the present results is in general a notably high value in the evaluation of atmospheric flow models.Table 2
Calculated values of the validation metrics from the LES results. The acceptance criteria for the model evaluation are adopted from [42].
Validation metric Present model value Acceptance criteria Satisfy? FAC2 0.9194 FAC2 ≥ 0.50 ✓ FB 0.1130 F B ≤ 0.3 ✓ NMSE 0.2180 NMSE ≤ 3 ✓


## 4.1. Results from the Reference Location

Figure4 compares the LES profiles of the mean velocity magnitude S and the turbulent kinetic energy k (TKE) with the field measurements at the reference location M0. The results are normalized by the frictional velocity u∗0 at M0. Here, the mean flow profile is very well reproduced by the recycling approach, as the simulated profile agrees well with the measurements at all heights. The TKE is also well predicted by the simulation at 5 m, but it is underestimated by about 25% at the height of 12 m. The total (modelled and resolved) TKE also underestimates the measured TKE at 12 m. The measurements show almost constant level of the TKE at M0, but in the simulation, the TKE decreases with height. Earlier, the LES runs submitted in the blind comparison also did not predict the constant TKE profile. This indicates that the deficiency could be very typical for LES modelling in general and thus occurs in the preset model, too. Since the modelled part of the TKE is significant in the area very close to the surface, we present the total TKE in all the results (Figures 6–13) hereafter. The error bars in the measurements show the standard deviations of the respective mean values. The standard deviation values were provided together with other statistics in the field data [29, 30].Figure 4
The LES profiles of the mean velocity magnitude (a) and the TKE (b) compared with the field measurements for the270∘ wind direction.
(a) (b)


## 4.2. Results on the Hill

Here, we first visualize the LES results on the instantaneous and mean flow structures around the hill. Figure5 depicts the instantaneous velocity fields as well as the 10-min time-averaged ones. It can clearly be seen that the local flow field is strongly affected by the hill. In particular, the cliff is noted to enhance the turbulent fluctuations. The complexity occurring in the flow due to the cliff is also illustrated by the streamlines which highlight the mean (10-min) flow patterns around the hill (Figure 5(b)). As expected, the flow is found to be separated at several places just after the cliff as well as on the lee sides. Particularly, the cliff is forming the small zone of the reversed flow, near the ground and along the entire cliff. Masts M2 and M6 are highly influenced with this reversed flow. However, the previous studies [6, 35] including those wind tunnels [34, 36–38] did not predict this flow separation nearby the cliff.Figure 5
LES-predicted instantaneous stream-wise velocity fields on the vertical plane (a), and the 10 min time-averaged streamlines along the hill (b). The results are from the270∘ wind direction.
(a) (b)Figure 6
Scatter plots of field data against LES results; solid lines: interval of±10% around perfect fitting (dashed line).
(a)
Mean velocity (b)
TKEFigure 7
Ratio of LES to field data for the normalized velocity magnitudeS/u∗0. The results from both wind directions are sorted with respect to the height above the ground level zagl (a) and the mast measurement locations (b).
(a) (b)Figure 8
Ratio of LES to field data for the normalized TKEk/u∗02. The results from both wind directions are sorted with respect to the height above the ground level zagl (a) and the mast measurement locations (b).
(a) (b)Figure 9
Vertical profiles of the velocity speedup (ΔS) at different masts (M1–M8) in the case of 270∘ (a) and 239∘ (b) wind directions.
(a) (b)Figure 10
Vertical profiles of the TKE increase (Δk) at different masts (M1–M8) in the case of 270∘ (a) and 239∘ (b) wind directions.
(a) (b)Figure 11
Simulation errors:RS for speedup and Rk for TKE increase errors in the case of 270∘ (a) and 239∘ (b) wind directions.
(a) (b)Figure 12
Simulation errors:RS for speedup and Rk for TKE increase errors at different measurement locations. The errors (RS and Rk) are averaged over the height for each location.Figure 13
Simulation errorsRS for speedup and Rk for TKE increase at different measurement heights. The errors (RS and Rk) are averaged over all the locations (masts) having the same height.Next, the results comparison between the simulation and the field data at different heights of all 10 masts (total 37-38 locations) is discussed. In order to see the correlation of the results, the scatter plots of all the available field data against the LES predictions for the mean velocity magnitude and the TKE are shown in Figures6(a) and 6(b), respectively. In the plots, an error interval of ±10% is also drawn for reference purpose. In Figure 6(a), it can be observed that 82% of the LES-predicted mean velocity in the case of 270∘ is within the interval of ±10%. The velocity prediction is almost the same, 81%, in the case of 239∘. However, a very small deviation (about 9%) is observed between the two cases when examining the turbulence part. For the TKE, 92% and 83% of the predictions appear within the interval (±10%) in the cases of 270∘ and 239∘, respectively. When one considers the interval of ±5%, the TKE prediction is still good, 80%, in both wind directions.To describe the discrepancies between the simulations and the measurements in detail, more scatter plots on the velocity and TKE predictions are presented in Figures7 and 8, respectively. In these figures, the results are sorted with respect to the measurement height above the ground level zagl (Figures 7(a) and 8(a)) and with respect to the masts (Figures 7(b) and 8(b)). The locations of the masts can be found from Figure 1. According to Figure 7(a), the predicted mean velocity is scattered for the height below 2 m, reaching a maximum error of almost 128% at M2 in the case of 270∘. On the other hand, the largest error (at 2 m) is smaller by about 35% in the case of 239∘. Generally, in both cases at zagl=2, masts M2, M6, and M8 have larger errors than the rest of the masts, as can also be seen in Figure 7(b). The large errors are due to the mast locations in which the flows are highly turbulent with flow separations (see Figure 5(b)). In both cases, the maximum error is reduced by about 100% within a distance of 3 m, from zagl=2m to 5 m. Above the height of 5 m, the prediction is greatly improved, as the maximum error is only 11%, while the majority of the error is below 5%.The overall prediction of the TKE is also scattered in both cases (270∘ and 239∘), although the level of maximum error is reduced by about 90% (Figure 8(a)) as compared to that of the velocity. The largest discrepancy is again found at M2, which is 70% (at 3.6 m) in the case of 270∘ and 65% (at 2 m) in the case of 239∘. However, unlike the velocity, the scatter is relatively independent of the mast locations and measurement heights (Figure 8(b)). Despite having less error (as compared to velocity), the TKE prediction does not seem to be improved with the height from 2 m to 5 m (Figure 8(a)). Overall, the velocity prediction is better in the case of 239∘, whereas the TKE prediction is better in the case of 270∘. More quantitative comparison between the two-wind-direction results is presented in Section 4.3.Following Bechmann et al. [30], the numerical results are further analyzed in terms of velocity speedup ΔS and TKE increase Δk. The speedup and the TKE increase are defined as follows [30]: (4)ΔS=Szagl-S0zaglS0zagl,Δk=kS02zagl-k0S02zagl,where the subscript 0 indicates the values (S0 or k0) from the reference mast M0. Vertical profiles of the speedup and the TKE increase at all the masts are shown in Figures 9 and 10, respectively.Here, it can be observed that the LES prediction of the speedup is generally acceptable (error less than 10%; Figure11(a)) at most (70%) of the places for both wind directions. Nevertheless, for the 270∘ wind, the simulation overestimates the measurements at the lowest heights of M2, M5, and M6, while it underestimates them at zagl=5 of M8. The maximum speedup error at the lowest height (2 m) is 53% and it decreases to 25% at 5 m in the case of 270∘ wind. Concerning the lowest heights of M2, M5, and M6, the speedup is predicted better in the case of 239∘, as the maximum error is only about 34%, which is at 2 m height of M2 (see Figure 11(b)). We note that the vertical profiles of the predicted speedup at M5 and M6 (below 5 m) differ significantly between the two directions. At the rest of the masts, there is no such a significant difference with respect to the wind direction. This means the flows in those regions (M5 and M6) are highly sensitive to the wind direction. According to the measurements, the speedup prediction in those areas is better in the case of 239∘.The large discrepancies between the simulations and the measurements at M2, M5, M6, and M8 are justified by the complex nature of the flows locally, that is, highly turbulent and reversing with several flow separations generated by the steepness of the hill. Particularly, masts M2 and M6 were placed about 10 m away from the cliff and were strongly influenced by the very small flow separation (due to the cliff) along the edge of the cliff, as it can be seen in Figure5(b). This particular separation is indeed much challenging to be even captured by numerical simulations because of its small size and location (near the cliff). Consequently, none of the previous attempts, including the wind-tunnel tests, have predicted this flow separation. Masts M8 and M4 were located on the lee side of the hill for 270∘ and 239∘ winds, respectively, and were also continuously influenced by flow separations. The simulation does predict most of the flow separations there, but the depth of the flow separation is somewhat compromised, although more at M8 than at M4. Previous studies [6, 30, 35–38] also reported large discrepancies at several masts, except for M1 and M7.Due to the flow separations, the turbulence level is significantly increased in a few locations. As a result, a clear peak in the TKE increase profile is visible at M2 and M6 for both winds and at M5 for the270∘ wind (see Figures 10(a) and 10(b)). Importantly, the simulation appears to accurately predict the extreme TKE vertical gradient, especially in the case of 270∘, including the lower heights (Figure 10(a)), even though the velocity speedup is highly overestimated at the same locations (i.e., M2 and M6). With the exception of M2 and M6, the TKE increase is also reproduced satisfactorily in the case of 239∘ direction. At 2 m height of M2 (for 239∘), the simulation shows the largest error of 140% for predicting the TKE increase (Figure 11(b)). Most of the discrepancies between the simulation and the measurements, visible in Figures 9 and 10, and especially at M2, M6, and M8, may also be justified by the uncertainties in the measured statistics. The high standard deviations at the lowest heights of M2 and M6 indicate high variabilities of the measured data at those places.


## 4.3. Simulation Error and Validation Metric

To evaluate the performance of the LES model in a quantitative way, the simulation (modelling or prediction) error from the LES results is calculated by considering the field data as a reference. Originally, Bechmann et al. [30] used the following equations to calculate the simulation errors for speedup and TKE increase. The error of the simulated speedup, RS, is defined as follows:(5)RS=100ΔSs-ΔSm,and the error of the simulated TKE, Rk, is(6)Rk=100Is-I0s-Im-I0mI0m,where ΔS is the speedup, I=k/S0, subscript zero denotes that the value is taken at the reference locations, and subscripts m and s denote the measured and simulated values, respectively. Figure 11 shows the simulation errors for the speedup and the TKE increase at each individual mast as well as at each possible height for the both wind directions. In the same way, Figure 12 displays the simulation errors, but they are averaged over the heights of each mast. The individual error from each height and/or location in Figure 11 may also be used for comparison with the results shown in the previous figures.In the following, we turn our focus on the averaged error to discuss the results. According to Figure12, the speedup error RS is less than 10% for half (50%) of the masts in both cases. The maximum error in the case of 270∘ is 27%, which is found at M6. On the other hand, the maximum error for the 239∘ wind is found at M5 and it is 18%. For the TKE increase, the error (Rk) is much higher than the speedup error (RS), reaching 59% (at M2) and 80% (at M6) in the cases of 270∘ and 239∘ winds, respectively. The high error for the TKE increase is obviously because the second-order statistic, for example, TKE or the Reynolds stress, is challenging to predict with enough accuracy as compared to mean flow (first-order statistic) in general. An interesting fact noted here is that the TKE increase error is almost the same at M2 for both winds, but a significant difference (75%) is observed at M6 due to the change in wind direction. Despite the difficulty of predicting the TKE compared to the mean wind, the Rk error is found to be smaller than RS at M6 and M8 in the case of 270∘ and at M5 in the case of 239∘. The high error at certain masts, such as M2 and M6, is due to their locations (near the cliff), where the aforementioned extremely complex flow behaviors occur. Indeed, virtually all previously reported models, including those wind tunnels, have failed to predict the flow accurately enough at masts M2 and M6. When excluding both M2 and M6, the averaged speedup error changes from 10.3% to 10.2%, which is insignificant. However, a high reduction from 24.1% to 13.5% is observed for the turbulence quantity. In addition, among all the masts, the standard deviation of the measured TKE is also at its highest at M2 and M6, reaching up to 24% in the case of 239∘, as could be expected from the error bars in Figures 9 and 10. The standard deviation or any other uncertainties of the measurements are, of course, not taken into account in the calculations of the simulation error.Next, Figure13 shows the (averaged) simulation error at three different heights. Here, the error is averaged over the measurement locations having the same height. For the 270∘ wind direction, it can be seen that the errors for the speedup and for the TKE increase at the height of 2 m are nearly the same (21%). This indicates that, at the lowest height, the model performs equally for both the speedup and the TKE increase. Above 2 m, the speedup error decreases rapidly down to 2.5% at 9 m. However, this does not happen with the TKE increase. The TKE increase error remains in the same order from 2 m to 5 m, and it decreases nominally, only by 4%, at 9 m.For the239∘ wind direction, the speedup error also decreases with increasing height, but not in the same trend as in the 270∘ case. It reduces from 13.5% at 2 m to 5.5% at 9 m. The TKE error is much higher, about 41% at 2 m (Figure 13). At this height, the TKE error is almost double compared to that of the 270∘ case. The high error is mostly due to the poor predictions at M2 and M6. At 5 m, the TKE error reduces greatly to 21%, about the same level as in the 270∘ case; however, it does not reduce further at 9 m.It should be emphasized that the measurements of the TKE at 2 m and 9 m heights of the reference location (M0) were missing in the field data (see Figure4(b)). Therefore, by following [30], a constant value of the TKE is considered at all three heights (2 m, 5 m, and 9 m) to calculate the measured TKE increase. This assumption (constant TKE with height) can highly bias the errors especially at 2 m, as the LES results follow the height-dependent TKE values from the simulation. Moreover, it should be remembered that the flow within the surface layer, for example, zagl≤2m in the Bolund case, may be of little interest for a wind-farm developer, as one is considering the wind condition within the rotor diameter. With the average over both the wind directions, the speedup error at 5 m is 9.3% and is reduced to 3.9% at 9 m. The error is not completely negligible, but the prediction could be considered accurate (less than 10% error) [30, 35].Table1 shows the mean errors from the two wind directions, and it compares the errors from the present case with those reported previously in the literature. The mean simulation error of the present LES is 10.3% and 24.1% for speedup and TKE increase, respectively. The LES prediction of the speedup is among the best results, as the mean speedup error from our case is only 0.7% higher than that of the best results [35]. Concerning the turbulence aspect, the present results show far smaller error in comparison with all previous results. Indeed, all prior models, except a RANS two-equation model in [30], are suffering a high TKE increase error. A couple of studies [6, 35] did not report the TKE increase error at all. The TKE increase error of 24.1% from this case is about 19% smaller than that of any other model. To our knowledge, this is thus the best prediction of the TKE over the Bolund Hill among all models employed over the hill so far.Finally, the present LES model is evaluated here using three different validation metrics. First, the validation metric called the fraction of two (FAC2) shows the fraction of data that is within the interval:(7)0.5≤PO≤2,where P is the predicted and O is the observed value. For an ideal model prediction, one would have FAC2 = 1 [43]. Second, the Fractional Bias (FB) can be defined as(8)FB=2·O-PO+P,where the angular bracket 〈·〉 represents the average value over the data samples. FB indicates the relative amount of bias in the modelled results such that FB =0 means no bias at all and FB >0 means systematic underestimation [44]. Finally, a third metric, known as the Normalized Mean Square Error (NMSE), is defined as follows:(9)NMSE=O-P2OP.NMSE is sensitive to large disagreements, in particular those occurring on the highest values [44]. The validation-metric values from the present results for both wind directions are listed in Table 2. According to the table, the present LES model passes all three validation tests, using the acceptance criteria for the model evaluation (in a rural area) proposed by Hanna and Chang [42]. In particular, the FAC2 value obtained by the present results is in general a notably high value in the evaluation of atmospheric flow models.Table 2
Calculated values of the validation metrics from the LES results. The acceptance criteria for the model evaluation are adopted from [42].
Validation metric Present model value Acceptance criteria Satisfy? FAC2 0.9194 FAC2 ≥ 0.50 ✓ FB 0.1130 F B ≤ 0.3 ✓ NMSE 0.2180 NMSE ≤ 3 ✓


## 5. Conclusions

This paper has focused on simulating the atmospheric flows over natural complex terrain and on validating the LES model against full-scale measurements. Previously, several researchers utilized LES-based numerical models to validate their simulations, especially over Askervein Hill, which is characterized by a smooth terrain. However, in this paper, we have selected much challenging complex terrain, the Bolund Hill, for validating the LES methodology. Unlike most of the previous studies on the Bolund case, we have carried out LES calculations for two different wind directions. The results on the mean and turbulence (e.g., TKE) quantities from the two wind directions have been analyzed equally and reported in detail herein. The results from both wind directions have been presented simultaneously and discussed quantitatively in order to learn the differences and the similarities in the results due to the change in wind direction. Using the simulation error and the model-validation criteria, the paper has evaluated and compared the performance of the LES model with other models previously employed over the hill.The recycling method reproduced the correct mean velocity profile at the reference location (M0) but showed some difficulties in the prediction of especially constant TKE profile. Unlike the field data, the LES-predicted TKE profile decreases gradually abovezagl=2.5m; however, it reproduced the correct (measured) value at 5 m. Bechmann et al. [30] mentioned that all the LES models that participated in the blind test did not reproduce the correct inflow (velocity and TKE) profiles. This indicates that the reproduction of the correct inflow profiles, and especially the constant TKE, in the Bolund case is a global challenge to LES modelling in general. However, in comparison with other prior LES models [30, 35], the prediction of the inflow conditions is better in this study thanks to the recycling technique.Generally speaking, the overall predictions of the speedup and the TKE increase over the hill are in close agreement with the measurements at all the positions, with the averaged difference of 10.2% for the speedup and 13.5% for the TKE increase, excluding M2 and M6. Due to the complex locations of M2 and M6, the simulation suffers from a maximum error of 140% to predict the measured values at 2 m height. Indeed, the measurements also show high standard deviations (up to 24%) at the same places, and this indicates that the measured data have high variabilities at those places (M2 and M6). Above the height of 5 m, the LES prediction of the mean flow is convincingly accurate, as the maximum difference between the simulations and the measurements is only 11%. In particular, the successful prediction of the reversed flow at the escarpment edge is believed to be an excellent achievement of this particular LES model, as no other previous numerical simulations had predicted the separation. For the turbulence part, the simulation reports the best prediction of the TKE, with the averaged error of 24.1%, which is 19% smaller than any other models, including the wind-tunnel experiments, reported in the past. Hence, it can be concluded that the present LES model, using the recycling technique on the inlet side, is able to reproduce the complex turbulent structures of the local wind flow over a complicated terrain such as Bolund Hill. Thus, it is possible to employ the same LES methodology to analyze the wind structures over other real terrains.Currently, this work is being extended with the inclusion of the full depth of the atmospheric boundary layer with capping inversion and the Earth’s rotation into the simulations. In the near future, the aim is to study the influence of capping inversion and the Coriolis force on the local atmospheric flows over this challenging Bolund site. The new results will then be compared with the existing results presented in this paper.



---

*Source: 9232759-2016-10-24.xml*