# Implicit Large eddy simulation of wind flow over rough terrain

Bofu Wang, Fei Cheng, Hui Xu, Xuerui Mao

► To cite this version:

Bofu Wang, Fei Cheng, Hui Xu, Xuerui Mao. Implicit Large eddy simulation of wind flow over rough terrain. 17th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery (ISROMAC2017), Dec 2017, Maui, United States. hal-02981011

HAL Id: hal-02981011

https://hal.science/hal-02981011

Submitted on 27 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Implicit Large eddy simulation of wind flow over
rough terrain

Bofu Wang<fnref n="1" />, Fei Chen<fnref n="2" />, Hui Xu<fnref n="3" />, Xuerui Mao<fnref n="1" />*

<figure><img src="image_2.png" /></figure>

ISROMAC 2017

International
Symposium on
Transport Phenomena
and
Dynamics of Rotating
Machinery

Maui, Hawaii

December 16-21, 2017

**Abstract**

Atmospheric boundary layer (ABL) flow over a rough terrain is studied numerically. A drag force model is employed to represent surface roughness. The spectral vanishing viscosity (SVV) method is adopted for implicit large eddy simulation (ILES). These techniques enable us to capture the dominant features of turbulent flow over non-smooth terrain with inhomogeneous roughness at high Reynolds number without resolving the very fine scale turbulence near the wall. We then run a test case of wind resource assessment on rough terrain, namely ABL flow over a hill with small-scale surface roughness at $Re = 10^6$ ($Re$ is defined based on hill height). The instantaneous flow, mean flow and Reynolds stress are analysed in order to interpret the obtained turbulent flow. The effect of roughness length and associated drag coefficient are also evaluated.

**Keywords**

Large eddy simulation — Rough terrain — Wind energy — Spectral element method — Spectral vanishing viscosity.

## 1. INTRODUCTION

Modelling atmospheric turbulent flow over a complex terrain with roughness is of great interest for many engineering, environmental, and geophysics applications, such as wind energy assessment; weather predictions; wind effects on agriculture and forest; and transport and dispersion of pollutants in the atmospheric environment. Specially, for wind energy application, accurate wind resource assessment is critical in the early stage development of wind farms. Inland wind turbines are usually sited on complex rough terrains with grass, forests, hills and mountains which make the evaluation of wind resources extremely challenging as the height of roughness is commensurate with the turbines. Moreover, the Reynolds number of turbulent Atmospheric boundary layer (ABL) is very high (up to $10^8$), imposing extra challenges to numerical assessment of wind resources and wind turbine performances. During the past decades, serious attempts have been made to predict the flow over complex terrain with roughness through conducting field experiments, employing physical models in wind-tunnel experiments, and carrying out numerical simulations. The most widely studied case is flow over forested hills.

The area with forested hills is commonly difficult to access, which prevents field measurements by either masts or LiDAR. A well-known hill experiment was conducted by an international group on the Askervein Hill in the UK [1, 2]. The slope of this hill is relatively smooth and the mean wind profiles are extensively used for research with particular interest in the so-called linearised model. This model does not

include detachment and recirculation under gentle slope assumption, and it reproduces quite satisfactorily the measured speed-up on Askervein at the hill top [3] but performs less successfully on the lee side. Another prominent measurement is the Bolund field campaign [4, 5, 6] over Bolund Hill in Denmark which provided new experimental datasets of the mean flow and the turbulence properties. The hill has a elatively smaller scale (12 m high) 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. Up to now, these two measurements have been widely used as complex-terrain benchmark cases for validating experimental and numerical implementations as well as evaluating turbulence models.

Regarding surface roughness, in early times, some model experiments were performed in wind-tunnel to obtain turbulent flow fields within and above a model plant canopy. Finnigan and Brunet [7] carried out a "Furry hill" experiment with a two-dimensional ridge using and provided a first view of the main features of canopy flow over a hill. Neff and Meroney [8] conducted single-hot film anemometer measurements to evaluate how hill and vegetation affect wind power availability. More recently, Poggi and Katul [9, 10, 11] conducted a series laser Doppler anemometry measurements over a train of gentle and narrow forested hills in a flume facility. Their results confirmed the theoretically predicted intermittent recirculation region downstream of a hill cov-

<footnote n="1">¹Faculty of Engineering, University of Nottingham, Nottingham, NG7 2RD, UK</footnote>

<footnote n="2">²Xinjiang Goldwind Science & Technology Co., Ltd., Beijing, 100176, China</footnote>

<footnote n="3">³Department of Aeronautics, Imperial College London, London SW7 2AZ, UK</footnote>

<footnote>*Corresponding author: xuerui.mao@nottingham.ac.uk</footnote>
ered by a dense canopy. They also observed sweep motions dominating momentum transfer within and near the canopy top all along the hill, while ejections dominating further above. Teunissen et al. used two wind tunnels and constant-temperature hot-wire anemometry to study the Askervein hill at three scales: 1/800, 1/1200, and 1/2500. Turbulence statistics gave more scatter compared to field observations with more than 50% error [12]. Discrepancies observed were higher in the wake and weaker at the location of the maximum speed-up. Recently, Yeow et al. [13, 14, 15] performed a series of wind-tunnel analyses with hot-wire anemometry and particle image velocimetry to study flow around the Bolund island with a 1/115 scale model and gave details on the unsteady flow behaviour. There are also some other experiments with complex topographies inducing complex flow phenomena, such as Chock and Cochran [16].

Measurement from meteorological towers and model experiments are currently indispensable in industrial practices, but not sufficient to analyse the intricate wind profile on complex terrains. Numerical models have the potential to offer insight with high fidelity information in this regard. Wind flow in complex terrains exhibits complex phenomena such as flow stagnation, rapid acceleration deceleration, recirculation-reattachment, vortex shedding, etc. Such micro-scale events occur at higher frequencies as well as at smaller spatial scales compared to meteorologic scale events, which imposes challenges to numerical simulation. Many efforts have been devoted to the behavior of the flow phenomena over terrain hills using Reynolds averaged models [17, 18], which predict the mean flow fields using far less computational resources. However, such models are based on a variety of assumptions and intrinsically ignore the flow unsteadiness; they are also known to be inaccurate in the flow separation region. In recent years, due to the development of numerical techniques and high-performance computing facilities, more advanced CFD methods, e.g. large-eddy simulation (LES) and direct numerical simulation (DNS) [19, 20] are employed in wind energy. It has been demonstrated that the LES techniques are able to reproduce many observed features of turbulent flow over flat terrain with homogeneous roughness [21] and downwind of forest leading edges [22]. There are investigations applying LES to flow over a rough hill (see [23] and references therein). Recently, developing LES models in the open source software OpenFOAM [24] becomes very popular for wind energy application [25, 26]. Three dimensional calculations for the Bolund hill case have been performed for validations. Improved comparison with measurements [4, 5] were obtained.

This paper is oriented towards applying LES for ABL flows over complex terrains at very large Reynolds number. The boundary layer flow over an idealistic hilly terrain is studied. This work paves the way for wind resource assessment over an area with a mixed scale of roughnesses.

## 2. METHODS

The governing equations for the current problem are the incompressible NS equations with a volume force term representing the force induced by the roughness,

$$ \partial_t u + (u \cdot \nabla)u = -\nabla p + Re^{-1} \nabla^2 u + F \quad \text{with} \quad \nabla \cdot u = 0, \quad (1) $$

where $u$, $p$ and $F$ are velocity, pressure and volume force, respectively. $Re$ is the Reynolds number defined based on the free-stream velocity and the height of the hill. $Re = 10^6$ is used throughout this work if not otherwise specified.

### 2.1 Spectral element method

A majority of numerical tools for wind flow simulation are based on finite volume methods. There are also some examples of successful applications of element-based Galerkin methods in numerical weather prediction [27]. In this study, we adopt the open-source software Nektar++ [28] which is designed to support the development of high-performance scalable solvers for partial differential equations using the spectral/hp element method. As a high order finite element method, it can deal with arbitrary geometric complexity, and are capable of local mesh adaption by either increasing the number of elements or increasing the polynomial order within elements. Moreover, it has excellent scalability for parallel simulation and is scalable to more than one million processor cores.

### 2.2 Roughness model

In LES applied in offshore wind energy, the roughness (waves) is considered to be sub-scale of the computational grid, and therefore can be modelled using wall turbulence models, which are commonly in the form of the boundary condition for shear stress. However regarding the deployment of wind turbines on rough terrains, the scale of the roughness elements, e.g. buildings and forest, is commensurate with the turbine and these elements cannot be considered as sub-grid. In the present work, a drag force approach is introduced to model roughness. The drag term can be expressed as:

$$ F_i = C_d A_f \sqrt{u_j u_j} u_i, \quad (2) $$

where $i=1, 2, 3$ refers to streamwise, lateral and vertical components respectively, $C_d$ is the mean drag coefficient of the roughness, $A_f$ is the density of roughness, and $u_i$ is the wind-velocity component. Furthermore, the dimensionless roughness length is defined as $z_0$. Both $C_d$ and $A_f$ are assumed to be independent on wind velocity. In our following simulations, we set $A_f = 1$, $C_d = 0.4$ and $z_0 = 0.1$ if not otherwise specified. Without loss of generality, the roughness is assumed to be with a uniform height.

### 2.3 Spectral vanishing viscosity method

The spectral vanishing viscosity (SVV) method is adopted to account the small-scale turbulence. The SVV define a viscosity in the spectral space that reaches a maximum at high wave numbers, but vanishes for wave numbers below a
resolution-dependent threshold [29]. In this approach, one adds an additional reaction term of the form to equations (1)

$$ \epsilon_{SVV} \frac{\partial}{\partial x} (\hat{Q} * \frac{\partial u}{\partial x}), \qquad (3) $$

where $\epsilon_{SVV}$ is a constant, * denotes the convolution operator, and $\hat{Q}$ is a kernel dictating which modes receive damping. This approach has been widely applied as a tool for implicit large eddy simulation (ILES). The key point in SVV filtering is that due to the shape of the kernel $\hat{Q}(k)$,

$$ \hat{Q}(k) = \begin{cases} \exp(-\frac{(N-k)^2}{(M-k)^2}), & M < k \le N \\ 0, & k \le M \end{cases} \qquad (4) $$

artificial viscosity for any mode number k is only applied above a cutoff mode M. For these damped modes, the total viscosity can thus be expressed as $1/Re + \epsilon_{SVV}$.

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic of the computational domain.</figcaption></figure>

### 3. DISCRETIZATION

The model of ABL flow over a small hill with uniform roughness near the wall is shown in Figure 1. The hill is quasi-three-dimensional (the spanwise direction is assumed to be homogeneous and discretised by Fourier transformation) and only a two-dimensional (2d) plane is shown here. The shape of the hill is sinusoidal with height $H = 1$ and width $W = 2$. At inflow, a power law velocity profile is applied. In order to trigger transition, a random perturbation with 5% turbulence intensity is applied at $y < 2$ in most of our calculations. A high order outflow condition is used at outlet [30]. At the bottom, no-slip wall condition is applied combined with a uniform roughness with length $z_0$. Far field condition is used at the top. The computational domain spans from -60 to 50, 0 to 50 and 0 to $\pi$ in the streamwise, vertical and spanwise directions, respectively.

Figure 2 shows the mesh generated by gmsh. There are 2815 quadrilateral elements in each 2d plane, and at least one element is set within the rough region near bottom. In each element, a spectral method is used to further decompose the element to a $(P+1) \times (P+1)$grid, where P represents a polynomial order. P = 6 is adopted in all our following simulations. In the homogeneous direction, a Fourier decomposition is adopted and 32 Fourier modes are computed.

<figure><img src="image_4.png" /><figcaption>Figure 2. Mesh in each x - y plane.</figcaption></figure>

### 4. RESULTS AND DISCUSSION

#### 4.1 Instantaneous flow field

<figure><img src="image_5.png" /><figcaption>Figure 3. Instantaneous contours of the streamwise velocity on horizontal cross sections (x-z planes) at (a) y=0.05, (b) y=0.1, (c) y=0.5, (d) y=1 and (e) y=2.</figcaption></figure>

We first evaluate the instantaneous flow fields. Figure 3 shows contours of streamwise velocity at various heights. Due to the roughness effect, upstream has already developed to turbulence. When this turbulent flow past the hill a turbulent wake region develops, the wind speed is reduced. The upstream streamwise velocity is characterised by longitudinal patterns which can be regarded as the signature of coherent structures induced by the bottom roughness. With increasing height, these structures increase in size and look like coherent streaks as usually observed in boundary layers. With further increasing height, fluctuations gradually decrease and the laminar ambient flow is recovered. The effect of the hill on instantaneous wind flow velocity is strong at the near wall region as shown in Figure 3 (a). When reaching the hill, the longitudinal pattern is almost eliminated in the wake and then regenerated further downstream as shown in Figure 3 (b), (c), (d). The effect of the hill diminishes with increasing height which can be seen from the contrast between upwind and downwind velocity. In Figure 3 (d) which shows the contour at y=1 (hill height), stronger fluctuation can be observed manifesting higher turbulence intensity there. At y=2, the effect of the hill on instantaneous wind flow is no longer observable (see Figure 3 e).
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Averaged contour plots for (a) streamwise velocity and (b) spanwise velocity.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Mean streamwise velocity profiles with inflow perturbation applied at $y < 4$.</figcaption></figure>

## 4.2 Mean flow

Figure 4 gives streamwise velocity profile sliced at different streamwise positions $x=-25, x=-20, x=-15, x=-10, x=-5$ and $x=3$ as well as the 1/7 power law velocity profile. We can see the velocity profile does not have much variation with upstream position (from $x=-25$ to $x=-5$), which is remarkably different from turbulent boundary layer flow. The extracted profiles deviate from the power law profile in the near wall region where the velocity decreases a little due to the surface roughness. The velocity profile indicates the boundary layer does not grow until near the hill side, a feature resulting from the uniform roughness. The profile at $x=3$ shows an inflection point near boundary which is caused by the inverse pressure gradient in that region. From this profile, it is also observed that the velocity above the hill is slightly larger than free stream velocity due to the speed up. Note that one important feature so called low-level jet [31] usually observed in ABL is not captured in this case. We then increase the perturbation region at inflow to $y < 4$. The streamwise velocity profiles

can be found in Figure 5. We can see that there is a small zone with concentrated stream which forms low-level jet like profile. In meteorological studies, low-level jets are caused by a lot of mechanisms which is beyond the scope of current study.

Figure 6 presents x-y slices of time- and z-averaged velocity components. From Figure 6 (a), the streamlines are distorted by the presence of the hill and a separated region on the lower lee side of the hill is generated as can be expected. The flow reacts as if the hill and its separation region were a single smooth obstacle. In Figure 6 (b), the average spanwise velocity has a relative small magnitude and shows coherent structure similar in instantaneous field. This distribution suggests that the flow in spanwise direction is dominated by small scale structures.

The average pressure field is presented in Figure 7 (a). A large-scale horizontal pressure gradient appears due to topography, inducing the distortion of the mean flow. The minimum pressure is located in wake region followed by an adverse pressure gradient. The surface pressure shown in Figure 7 (b) also has a maximum before the hill. A minimum in surface pressure occurs as the flow passes the hill summit, and large-magnitude oscillation is observed near hill top. The variation of surface pressure is qualitative in agreement with the wind tunnel measurement [32].

## 4.3 Reynolds stress

Figure 8 shows the distribution of spanwise averaged streamwise normal Reynolds stress, $\langle u'u' \rangle = \frac{1}{\pi} \frac{1}{T} \int_{0}^{\pi} \int_{0}^{T} u'u'dzdt$, where $u'$ is perturbation from the time averaged streamwise velocity. It can be seen that the wake turbulence exhibit plume-like structure [33]. The region of maximum production is some distance along the downstream of the hill. The structure extends downwind above the region of reversed
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Average pressure distribution, (b) Surface pressure.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Reynolds stress $<u'u'>$.</figcaption></figure>

mean flow. This structure indicate high turbulence level in this region and it is associated with strong shear layer formulated on the top of the recirculation zone.

Figure 9 shows wall shear stress at bottom surface with varies roughness length $z_0$ and drag coefficient $C_d$. Take $z_0 = 0.1$, $C_d = 0.4$ as a reference case, either decreasing $z_0$ or decreasing $C_d$ will lead to an increase of wall shear stress. A high peak in wall shear stress is observed at hill top.

## 5. SUMMARY

In this paper, we have carried out ILES to investigate the ABL flows over a small hill at $Re = 10^6$. It is found that under the effect of surface roughness, the flow upstream of the hill behaves like a fully developed turbulent boundary layer. Streak-like structures can clearly be observed from instantaneous flow fields. Low-level-jet like velocity profile can be captured due to roughness effect. As the turbulent flow past the hill, a wake region is developed on the lee side of the hill. The wake region is characterised by a recirculation zone, a strong elevated shear layer downstream from the summit and high turbulence levels. The wake turbulence is characterised by plume-like structure. The roughness length

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Wall shear stress.</figcaption></figure>

and drag coefficient can strongly increase the wall shear stress.

## REFERENCES

[1] Taylor P. and Teunissen H. W. Askervein'82: report on the september/october 1982 experiment to study boundary-layer flow over askervein, south uist. *Ontario technical report MSRB-83-8*. *Atmospheric Environment Service, Downsview*, 1983.

[2] Taylor P. and Teunissen H. W. The askervein hill project: report on the sept./oct. 1983, main field experiment. *Ontario technical report MSRB-84-6*. *Atmospheric Environment Service, Downsview*, 1984.

[3] Hunt J. C. R., Leibovich S., and Richards K. J. Turbulent shear flows over low hills. *Q. J. Roy. Meteorol. Soc.*, 114:1435, 1988.

[4] Berg J., Mann J., Bechmann A., Courtney M. S., and Jørgensen H. E. The bolund experiment, part i: flow over a steep, three-dimensional hill. *Boundary-Layer Meteorol.*, 141:219-243, 2011.

[5] Bechmann A., Sørensen N. N., Berg J., Mann J., and Rëthorë P. E. The bolund experiment, part ii: blind comparison of microscale flow model. *Boundary-Layer Meteorol.*, 141:245-271, 2011.

[6] Bechmann A., Berg J., Courtney M., Jørgensen H., Mann J., and Sørensen N. N. The bolund experiment: Overview and background. *Danmarks Tekniske Universitet Risø National Laboratory for Sustainable Energy*, 2009.

[7] Finnigan J. J. and Brunet Y. Turbulent airflow in forests on flat and hilly terrain. in *Wind and trees, Coutts MP, Grace J* (eds). *Cambridge University Press: Cambridge, UK.*, pages 3-40, 1995.

[8] Neff D. E. and Meroney R. N. Wind-tunnel modeling of hill and vegetation influence on wind-power availability. *J. Wind Eng. Ind. Aerodyn.*, 74-76:335-343, 1998.

[9] Poggi D. and Katul G. G. An experimental investigation of the mean momentum budget inside dense canopies on narrow gentle hilly terrain. *Agric. Forest Meteorol.*, 144:1-13, 2007.
[10] Poggi D. and Katul G. G. Turbulent flows on forested hilly terrain: The recirculation region. *Q. J. R. Meteorol.* Soc., 133:1027–1039, 2007.

[11] Poggi D., Porporato A., Ridolfi L., Albertson J. D., and Katul G. G. The effect of vegetation density on canopy sub-layer turbulence. *Boundary-Layer Meteorol.*, 111:565–578, 2004.

[12] Teunissen H. W., Shokr M., Bowen A., Wood C., and Green D. The askervein hill project: wind-tunnel simulations at three length scales. *Boundary-Layer Meteorol.*, 40:1–29, 1987.

[13] Yeow T. S., Cuerva A., Conan B, and Pérez Álvarez J. Wind tunnel analysis of the detachment bubble on bolund island. In: *Proceedings of the science of making torque from wind*, Oldenburg, Germany, 2012.

[14] Yeow T. S., Cuerva A., Conan B, and Pérez Álvarez J. Wind tunnel analysis of the detachment bubble on bolund island. *Journal of Physics: Conference Series*, 555:012021, 2014.

[15] Yeow T. S., Cuerva A., and Pérez Álvarez J. Reproducing the bolund experiment in wind tunnel. *Wind Energy*, 18:153–169, Wind Energy.

[16] Chock G. and Cochran L. Modeling of topographic wind speed effects in hawaii. *J. Wind Eng. Ind. Aerodyn.*, 93:623–638, 2005.

[17] Finardi S., Brusasca G., Morselli M., Trombetti F., and Tampieri F. Boundarylayer flow over analytical two-dimensional hills: A systematic comparison of different models with wind tunnel data. *Boundary Layer Meteorol.*, 63:259–291, 1993.

[18] Ying R. and Canuto V. Numerical simulation of flow over two-dimensional hills using a second-order turbulence closure model. *Boundary Layer Meteorol.*, 85:447–474, 1997.

[19] E. Kleusberg, R. F. Mikkelsen, P. Schlatter, S. Ivanell, and D. S. Henningson. High-order numerical simulations of wind turbine wakes. *J. Phys. Conf. Ser.*, 854:012025, 2017.

[20] A. Ducoin, M. S. Shadloo, and S. Roy. Direct numerical simulation of flow instabilities over savonius style wind turbine blades. *Renew. Energ.*, 105:374–385, 2017.

[21] Su H. B., Shaw R. H., Paw U. K. T., Moeng C. H., and Sullivan P. P. Turbulent statistics of neutrally stratified flow within and above a sparse forest from large-eddy simulation and field observations. *Boundary-Layer Meteorol.*, 88:363–397, 1998.

[22] Yang B., Raupach M., Shaw R. H., Paw U. K. T., and Morse A. P. Large eddy simulation of turbulent flow across a forest edge. part i: Flow statistics. *Boundary-Layer Meteorol.*, 120:377–412, 2006.

[23] Dupont S., Bruneta Y., and Finnigan J. J. Large-eddy simulation of turbulent flow over a forested hill: Validation and coherent structure identification. *Q. J. R. Meteorol.* Soc., 134:1911–1929, 2008.

[24] Opencfd ltd. user and programmer's guide. http://www.openfoam, 2012.

[25] Churchfield M. J., Lee S., Michalakes J., and Moriarty P. J. A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics. *J. Turbul.*, 13:N14, 2012.

[26] Vuorinen V., Chaudhari A., and Keskinen J. P. Large-eddy simulation in a complex hill terrain enabled by a compact fractional step openFOAM solver. *Adv. Eng. Soft.*, 79:70–80, 2015.

[27] Simone Marras, James F. Kelly, Margarida Moragues, Andreas Müller, Michal A. Kopera, Mariano Vázquez, Francis X. Giraldo, Guillaume Houzeaux, and Oriol Jorba. A review of element-based galerkin methods for numerical weather prediction: Finite elements, spectral elements, and discontinuous galerkin. *Arch. Comput. Method E.*, 23:673–722, 2016.

[28] Cantwell C.D., Moxey D., Comerford A., Bolis A., Rocco G., Mengaldo G., Grazia D. De, Yakovlev S., Lombard J.-E., Ekelschot D., Jordi B., Xu H., Mohamied Y., Eskilsson C., Nelson B., Vos P., Biotto C., Kirby R. M., and Sherwin S. J. Nektar++: An open-source spectral/hp element framework. *Comput. Phys. Comm.*, 192:205–219, 2015.

[29] Tadmor E. Convergence of spectral methods for nonlinear conservation laws. *SIAM J. Numer. Anal.*, 26:30–44, 1989.

[30] Dong S., Karniadakis G. E., and Chryssostomidis C. A robust and accurate outflow boundary condition for incompressible flow simulations on severely-truncated unbounded domains. *J. Comput. Phys.*, 261:95–136, 2014.

[31] Wagner A. Uber die tageswinde in der freien atmosphire. *Grrl. Britr. Grophys.*, 1939.

[32] Brunet Y., Finnigan J. J., and Raupach M. R. A wind tunnel study of air flow in waving wheat: Single-point velocity statistics. *Boundary Layer Meteorol.*, 70:95–132, 1994.

[33] Keith W. Ayotte, Peter P. Sullivan, and Edward G. Patton. Les and wind tunnel modelling over hills varying steepness and roughness. *The Fifth International Symposium on Computational Wind Engineering*, 2010.