Turbulent Inflow Generation for Large-Eddy Simulation of Winds around Complex Terrain

Abstract

Accurate turbulent inflow conditions are needed to broaden the application of the large-eddy simulation technique to predict winds around arbitrarily complex terrain. We investigate the concept of buoyancy perturbations with colored noise to trigger turbulence in upstream flows approaching complex terrain regions. Random perturbations are imposed on the source term in the pseudo-temperature transport equation. These perturbations are effective within three-dimensional boxes and scaled using a bulk Richardson number defined for each box. We apply the turbulent inflow generation technique to predict winds around the Askervein and Bolund Hills under neutrally stratified conditions. We find that a common value for the bulk Richardson number works well for a variety of flow problems. Additionally, we show that the height of the perturbation box plays an important role in the accuracy of the predictions around complex terrain. We consistently obtained good results for both simulation cases when the perturbation box height was made a fraction of the Obukhov length scale.

1. Introduction

The large-eddy simulation (LES) is increasingly applied to simulate wind farm power production, wind turbine wake dynamics, and various micro-meteorological applications [1]. LES holds great promise to improve the estimation of optical turbulence parameters over regions that are relevant to the astronomy field as well [2]. The fundamental idea behind LES is to resolve turbulent eddies larger than a particular length scale while modeling the effect of eddies smaller than the said length scale on the larger eddies. The resolved flow field is then expected to exhibit the expected statistics of turbulence. LES needs lateral boundary conditions that admit arbitrary flow conditions to be able to simulate complex flow scenarios. Inflow flow must contain realistic turbulent correlations, but achieving this goal in a simulation is not straightforward. Turbulent inflow conditions have been the focus of numerous studies over the past two decades as evidenced by several review articles on the subject [3,4,5].

Wu [5] highlights the need for robust inflow generation techniques for LES of atmospheric flows that take into account surface roughness and the three-dimensionality of the mean wind flow, i.e., time-varying mean flow directions. Turbulent inflow generation techniques generally fall into two categories: precursor methods and synthetic methods [4]. In precursor methods, a separate auxiliary boundary-layer simulation feeds fully developed turbulent flow information into a target simulation at its inflow boundaries.

Precursor methods have been adopted for LES of atmospheric flows in several studies [6,7,8,9]. In a precursor method, accurate turbulent flow correlations can be introduced into the target simulation. The auxiliary simulation in the precursor method consists of obtaining a converged solution in a doubly periodic domain driven by either a constant mass flow rate [6] or a constant mean-pressure gradient [7]. Often, the auxiliary domain is stored as an independent database that can be accessed by the target simulation. The drawback of the precursor method, however, is having to perform a computationally expensive auxiliary simulation. Simulation turnaround time also slows down because of frequent data input and output.

Synthetic methods avoid the need for an auxiliary simulation and seek to generate realistic turbulence in the incoming flow as part of the target simulation. However, as turbulent flow correlations are typically not known in advance for general flows, synthetic methods require a development fetch before realistic turbulent flow conditions can be observed within the simulation domain. The promise of synthetic methods is that a development fetch shorter than a separate auxiliary domain would greatly reduce the computational overhead of generating turbulence.

As a subset of synthetic methods, recycling methods also avoid the creation of inflow databases. Lund et al. [10] proposed a recycling method for turbulent-boundary layers observed in the field of aerodynamics. In this particular method, the flow is sampled from several boundary-layer thicknesses downstream of the inlet and then re-introduced at the inlet after being scaled by inner and outer layer similarity laws to account for the effect of a growing turbulent-boundary layer. Mayor et al. [11] applied the recycling idea to develop homogeneous turbulence, which would then become the inflow condition for an inhomogeneous-flow field regarding meso-scale internal boundary layers. The original method of Lund et al. [10] used empirical relationships when applying the scaling arguments between the recycling location and the inlet. Araya et al. [12] proposed using a test plane to determine flow parameters dynamically. Recycling methods, however, do not readily extend to complex-surface boundaries, because similarity laws for such problems are not available.

Recent wind-farm simulation studies have found that applying mean shear and turbulent fluctuations independently can lead to erroneous results [1,13]. The synthetic method of Mann [14] constructs a spectral tensor by applying the rapid-distortion theory to linearized Navier–Stokes equations. While a turbulent-like state can be generated, it may not necessarily represent a realistic turbulence state. Regardless, the synthetic method of Mann has been used in a variety of ABL-related studies. Xie and Castro [15] proposed generating a two-dimensional turbulent plane using digital filters on unit-variance random data with a zero mean. Based on the observation that correlation functions in shear flows follow a decaying exponential, Xie and Castro correlated each two-dimensional plane with the plane from a previous time step using a decaying exponential weighting function. They were able to obtain good results for the simulation of flow over urban canopies.

Most synthetic methods rely on the information on the Reynolds stress tensor or turbulent-spectral content at the inlet domain in advance of the simulations. For instance, in their synthetic-eddy method (SEM), Jarrin et al. [16] defined the Reynolds stresses and length scales for turbulence by extracting that information from a separate Reynolds-averaged Navier–Stokes (RANS) simulation. Poletto et al. [17] improved upon the SEM by ensuring divergence-free conditions at the inlet. In their revised version of the SEM, Poletto et al. used an existing direct numerical simulation (DNS) database to provide the Reynolds stresses at the inlet. This issue of providing turbulence information at the inlet is important. A common metric when comparing different turbulent inflow generation methods is the development length to achieve a realistic turbulent state in the domain. A short development length is desirable, but it is important to keep in mind the initial conditions for each method when comparing different turbulent-inflow methods. When the initial state at the inlet is obtained from a separate turbulence simulation, the development length is typically short. However, we should keep in mind that detailed statistics of turbulence at the inlet may not be available for a general flow problem.

A synthetic turbulent-inflow method that does not require detailed information on turbulence statistics at the inlet can be achieved through artificial perturbations imposed on a pseudo-temperature field that is only coupled to the momentum equations over a short buffer zone placed downstream on the inflow boundary. This original idea was introduced by Munoz-Esparza et al. [18]. In their so-called cell perturbation method, white-noise perturbations, scaled by an Eckert ($Ec$) number argument, were added onto the potential temperature field within a buffer region near the inlet. These perturbations were updated at a specified frequency. These perturbations in the potential temperature field then modified the momentum field through the action of buoyancy to create three-dimensional vortices that broke down into turbulence after a development fetch. Eckert number scaling of the perturbations and a time-scale for the perturbation update frequency were calibrated in Munoz-Esparza et al. [19]. The cell perturbation method is easy to implement and does not require knowledge of turbulence statistics at the inlet. The method readily extends to problems with time-varying mean-flow directions.

Umphrey and Senocak [20] applied the cell perturbation method to a canonical incompressible turbulent-channel flow case. Umphrey and Senocak found that the Eckert number guidance does not extend to the turbulent channel flow case, which had a Reynolds number that was much lower than the corresponding value for atmospheric flows. Subsequently, Umphrey and Senocak modified the cell perturbation method by generating colored noise in the pseudo-temperature field. They also imposed perturbations in three-dimensional volumes instead of planar cells as done in the cell perturbation method. Deleon et al. [21] focused on eliminating the ad hoc nature of determining the maximum amplitude of temperature perturbations in the work of Umphrey and Senocak [20]. Deleon et al. [21] proposed a bulk Richardson number scaling argument to calculate the maximum perturbation amplitude. They also took into account the variation of the velocity with distance from the wall. Deleon et al. [21] applied their proposition to the canonical cases of turbulent-channel flow and the backward-facing step. In the backward-facing step case, they also simulated the spatially evolving boundary layer upstream of the step. They obtained excellent predictions and demonstrated that a common set of parameters can be used for different flow problems at different Reynolds numbers. Deleon et al. [21] referred to their approach as the box perturbation method to differentiate the differences from the original cell perturbation method [19].

Prediction and forecasting of winds around complex terrain with spatial resolutions that are on the orders of 10 m are desirable for the wind energy field in addition to several other applications. Stand-alone micro-scale wind simulation codes provide many advantages, but such models need to be informed by meso-scale numerical weather prediction models for accurate predictions. The complexity of underlying terrain, meandering wind directions, and high Reynolds number characteristics of atmospheric flows call for turbulent inflow generation techniques that are versatile and computationally inexpensive. Despite the promising results presented in the work of Deleon et al. [21], the box perturbation method has not been tested for turbulent inflow generations for atmospheric flows around complex terrain. In the present investigation, we apply and refine the box perturbation method to simulate wind over complex terrain at realistic Reynolds numbers using a hybrid LES-RANS turbulence model.

2. Numerical Formulation

We solve the filtered form of the incompressible Navier–Stokes equations

$$\begin{array}{c}\hfill \frac{\partial {\overline{u}}_j}{\partial x_j}=0,\end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial }{\partial x_j}\left({\overline{u}}_i{\overline{u}}_j\right)=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_i}& +\frac{\partial }{\partial x_j}\left(2\nu {\overline{S}}_{ij}-{\tau }_{ij}\right),\hfill \end{array}$$

(2)

where the ${\overline{S}}_{ij}$ is the strain rate tensor and ${\tau }_{ij}$ are the filtered Reynolds stresses, which are modeled with an eddy-viscosity assumption as follows

$$\begin{array}{cc}\hfill {\tau }_{ij}=& {\nu }_t{\overline{S}}_{ij},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\nu }_t=& l_{mix}^2\left|\overline{S}\right|.\hfill \end{array}$$

(4)

For the mixing length scale $l_{mix}$, we use the hybrid LES-RANS turbulence model described in Senocak et al. [22]. In this model, the length scale used in the LES eddy-viscosity model is mathematically blended with Prandtl’s mixing-length scale used in the RANS model

$$\begin{array}{c}\hfill l_{mix}=\left(1-exp\left(-h/h_B\right)\right)C_S\Delta +exp\left(-h/h_B\right)\kappa h,\end{array}$$

(5)

where h is the shortest distance to the surface and is not a constant value, $h_B$ represents the RANS-LES blending height, and $\Delta$ is the LES filter width.

As part of the box perturbation method [21] for turbulent inflow generation, we solve the filtered form of the temperature transport equation

$$\begin{array}{c}\hfill \frac{\partial \overline{\varphi }}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\varphi }{\overline{u}}_j\right)=\frac{\partial }{\partial x_j}\left(\Gamma \frac{\partial \overline{\varphi }}{\partial x_j}+{\tau }_{\varphi j}\right)+{\mathsf{\Phi }}_{PB},\end{array}$$

(6)

where ${\tau }_{\varphi j}$ is the filtered turbulent-heat flux, $\Gamma$ is the thermal diffusivity, and ${\mathsf{\Phi }}_{PB}$ is the source term to generate perturbations in the temperature field. The turbulent-heat flux is modeled by

$$\begin{array}{c}\hfill {\tau }_{\varphi j}=\frac{{\nu }_t}{P{r}_t}\frac{\partial \overline{\varphi }}{\partial x_j},\end{array}$$

(7)

where ${\nu }_t$ is the kinematic turbulent-eddy viscosity and $P{r}_t$ is the turbulent Prandtl number, which we set to 0.85.

Deleon et al. [21] proposed a bulk Richardson number ($Ri$) scaling for the source term that appears in Equation (6)

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_{PB}=\frac{R{i}_{PB}{〈U\left(z\right)〉}_{PB}^3}{g\beta L_{PB}{H}_{PB}},\end{array}$$

(8)

where g is the gravitational constant, $\beta$ is the thermal expansion coefficient for air, and $L_{PB}$ and $H_{PB}$ are the length and height of the perturbation box, respectively. ${\langle U\left(z\right)\rangle }_{PB}^2$ is the averaged incoming velocity calculated at the height of the perturbation box center measured from the surface. The bulk Richardson number $R{i}_{PB}$ for each perturbation box is prescribed as an empirical constant. We investigate the effect of $R{i}_{P}B$ on the wind field in the present study. For the derivation and reasoning behind Equation (8), we refer the reader to Deleon et al. [21].

Pseudo-random source term perturbations that are uniformly distributed in the [$-{\mathsf{\Phi }}_{PB},+{\mathsf{\Phi }}_{PB}$] are then applied within each perturbation box in the buffer zone placed downstream of the inlet face of the computational domain. We note that perturbing the source term, as opposed to directly perturbing the pseudo temperature field with white noise as done in the cell perturbation method, leads to colored noise in the velocity field as the pseudo temperature and velocity field are transported and additionally filtered under the action of the eddy viscosity. To ensure thermally neutral conditions within the simulation domain, all temperature perturbations are uniformly shifted every time step such that integration of the buoyant forces leads to a net-zero thermal-energy input. Source term perturbations on each box, which encompasses several computational cells in a three-dimensional volume, are kept constant over a time period $t_P$ determined based on the following guidance:

$$\frac{t_p{\langle U\left(z\right)\rangle }_{PB}}{L_{PB}}=1.$$

(9)

Each box is then perturbed according to Equation (8) at every $t_P$ throughout the entire simulation.

3. Computational Setup

The above formulation with the immersed-boundary method to represent complex terrain is described and validated in the works of Umphrey et al. [23] and DeLeon et al. [24]. We refer the reader to those studies for further details. The finite-difference method with a second-order accurate Adams–Bashforth scheme for time advancement and a second-order accurate central-difference scheme for spatial terms are adopted in the flow solver. The three-dimensional incompressible flow solver has been parallelized for fast computations on clusters of graphics processing units [25,26,27]. The pressure Poisson’s equation, which is the most time-consuming element of incompressible flow computations, is solved iteratively with an amalgamated, parallel, geometric-multigrid technique [28].

In the following simulation cases, the top boundary condition is a zero-shear-stress wall and the terrain surface boundary conditions are calculated with the immersed-boundary method presented in the study by DeLeon et al. [24]. In the simulations where periodicity is adopted in the streamwise direction, it is enforced using the shifted periodic boundary conditions [29], whereas regular periodicity is adopted in the spanwise lateral direction. For simulations with the inflow turbulence generation technique, a convective type outflow condition for the momentum field and a constant pressure value are adopted as boundary conditions at the outlet face of the computational domain.

4. Results

We consider the Askervein Hill and Bolund Hill cases to evaluate the box perturbation method as an inflow generation method for complex-terrain flows. The Askervein Hill case was the focal point of a field measurement campaign in the 1980s [30,31,32,33]. Bolund Hill is relatively a new case. It is an isolated hill with a steep escarpment surrounded by shallow waters. It is located in Denmark. It was the focus of a field study [34] and a blind evaluation study of different computer models [35]. DeLeon et al. [24] validated the current flow solver with shifted periodic boundary conditions [29] in the streamwise direction for both hills. In what follows, we adopt the same bulk Richardson number guidance as in the work of Deleon et al. [21], which was tested only for canonical non-atmospheric flow cases, and assess the effect of perturbation box height on predictions and develop guidance for it using the aforementioned complex-terrain cases.

4.1. Askervein Hill

The topography of the Askervein Hill case is shown in Figure 1. The domain size for the Askervein simulation is 8218 m × 5800 m × 1000 m in the streamwise, horizontal spanwise, and vertical directions with a grid resolution of 641 × 513 × 257 in the same respective directions. We use the Askervein terrain map that does not include neighboring hills [33] and immersed it such that the surrounding flat surface is located 14.7 m above the bottom of the domain. As suggested by Taylor and Teunissen [30], we set the aerodynamic roughness length $z_0$ = 0.03 m. We use a shift of 1000 m for the shifted periodic boundary conditions and drive the flow by dynamically calculating the mean-pressure gradient to keep a constant mass flow rate of approximately $9\times {10}^7\mathrm{kg}{\mathrm{s}}^{-1}$ using the technique described by Benocci and Pinelli [36]. The RANS-LES blending height, $h_B$, is set to 35 m.

We use the simulation with shifted periodic boundary conditions in the streamwise direction as a reference and compare the box perturbation method results against it.

In the box perturbation method, we prescribe a mean-velocity profile at the inlet based on the logarithmic law using $u_{*}$ = 0.68 $\mathrm{m}{\mathrm{s}}^{-1}$ for neutral stratification. The inlet profile was set to a constant value beyond 850 m above the ground level. A convective-outflow boundary condition is used at the domain outlet. The free-stream temperature of the air ${\varphi }_0$ was assumed to be 293 K and we set $\beta$ = ${\varphi }_0^{-1}$ in Equation (8).

In this study, we focus on the effect of perturbation box height on the predictions. The three box heights considered in our study are presented in

Table 1. Parameters of the box perturbation method used in the simulation of flow over Askervein Hill. $L_{PB}$, $W_{PB}$, and $H_{PB}$ are the length, width, and height of the perturbation boxes, respectively. Depth is the number of perturbation boxes in the streamwise direction.

Case	${\mathit{H}}_{\mathit{PB}}$ (m)	${\mathit{L}}_{\mathit{PB}}$ (m)	${\mathit{W}}_{\mathit{PB}}$ (m)	PB Depth
Ask1	12	23	23	3
Ask2	38	76	76	3
Ask3	62	123	123	3

. Based on the work of Deleon et al. [21] on canonical turbulent-channel and backward-facing step flows, we use the same $R{i}_{PB}$ and set it to 0.042 for all the simulations presented in this study, which worked well. Deleon et al. [21] recommended a perturbation box height, $H_{PB}$, of approximately one-eighth of the turbulent-boundary layer height. The perturbation box length, $L_{PB}$, and width, $W_{PB}$, were then selected as twice the $H_{PB}$.

For atmospheric flows, setting $H_{PB}$ to one-eighth of the boundary-layer height would result in large perturbation boxes. The perturbations at the inlet would then take a significant development fetch to break down into fine-scale turbulence. Therefore, we seek different guidance to determine the perturbation box sizes for complex-terrain flows. In Lopes et al. [6], an Obukhov length of 600 m was used in their analysis. Cases Ask2 and Ask3 in Table 1 are approximately one-sixteenth and one-eighth of that Obukhov length scale, respectively. Note that the actual box height used in the simulation is based on having a whole number of computational cells within that suggested height. As the turbulent stresses greatly affect the mean flow in the inner layer, we also consider a perturbation box height comparable to the height of the inner layer, expecting that a small buoyancy perturbation box will break down quickly. An inner-layer height of 11.6 m was used in the work of Lopes et al. [6]. In case Ask1, we use $H_{PB}=12$ m.

Turbulent-flow fields generated by the three perturbation boxes are shown in Figure 2. The visualization plane displays a wind speed of 100 m above ground level. As expected, the inflow with the smallest perturbation box (i.e., Figure 2b breaks down into a turbulent state over a short distance, whereas the inflows with larger perturbation boxes (i.e., Figure 2c,d) takes a longer development fetch to break down into a realistic turbulent state. Based on the overall domain size, perturbations in cases Ask1 and Ask2 break down into realistic turbulence long before the hill. In case Ask3, however, perturbations do not break down over such a short distance.

To assess the impact of the perturbation box height on the mean flow, we compare our results to the case of MF03-D taken between 1400–1700 BST (British Summer Time = UTC + 1 h) on 3 October 1983 [37]. We present the fractional speed-up ratio to compare our results to the field data. The fractional speed-up ratio is the horizontal velocity speed-up relative to an undisturbed mean-velocity profile,

$$\begin{array}{c}\hfill \Delta S=\frac{U\left(h_{agl}\right)-U_{ref}\left(z\right)}{U_{ref}\left(z\right)},\end{array}$$

(10)

where $h_{agl}$ and z are both vertical heights above ground level of the terrain, except z, which contains the global elevation datum.

Figure 3 shows the horizontal wind speed at the reference site (RS) and the fractional speed-up ratio at the hilltop (HT). The mean wind speed shown in Figure 3a, produced by the box perturbation method, agrees reasonably well with both the field data and the logarithmic-law wind profile. The box perturbation method produces better agreement with the logarithmic-law profile than it does with the periodic simulation. This is expected because the logarithmic-law profile is prescribed at the inlet. Upon close inspection, we see that predictions below $h_{agl}\approx 60$ m are sensitive to the perturbation box height choice, whereas above 60 m, the profiles collapse onto each other.

The fractional speed-up ratio in Figure 3b shows the effect of perturbation box height on how the flow accelerates over the hill. The two larger perturbation boxes show the best agreement with the field data above $h_{agl}\approx 15$ m. The smallest perturbation box is drastically lower in fractional speed-up ratio, which we also observe consistently in the case of Bolund Hill, as we will show in the next section. The flow structures are induced by synthetic buoyancy perturbations. The size of the flow structures appears to have a significant effect on the acceleration over the hill. A plausible explanation for this effect could be the reduction in momentum transfer from the free-stream flow toward the flow near the surface because of the absence of larger eddies in the free-stream flow. We should mention that we do not expect a good agreement with the field date very close to the surface, because we lack the vertical resolution needed to capture the variation of the flow as observed in the field experiment.

Figure 4 and Figure 5 show the fractional speed-up ratio along the AA and A lines at 10 m above ground level. On the windward side of the hill, both the two largest perturbation boxes and the periodic simulation show a good collapse of the predictions. On the other hand, the smallest perturbation box, which developed into realistic turbulence well before approaching the hill, shows a significantly lower fractional speed-up ratio. Again, we attribute this behavior to the lack of large-scale eddies in the simulation when a small perturbation box is adopted at the inlet. The periodic reference simulation and the two largest perturbation boxes do not show much sensitivity on the leeward side of the hill, as seen in Figure 4. The separation behind the hill along the AA line is not captured well in general due to the deficiencies of the current turbulence model for separated flows, as discussed by DeLeon et al. [24]. The leeward side of the A line experiences the full shadow of the hill. Consequently, a stronger flow separation than the separation along the AA line exists there. This high sensitivity of the results on the A line was also reported by DeLeon et al. [24] for different simulation parameters.

4.2. Bolund Hill

In the study by Deleon et al. [21], the same set of parameters worked well for the canonical cases of turbulent-channel flow and backward-facing step. We have already shown that the same bulk Richardson number for the perturbation box suggested in Deleon et al. [21] produces satisfactory results for atmospheric flow around Askervien Hill. However, we had to adopt different guidance to determine the size of the perturbation boxes. The guidance to size the perturbation boxes relative to the turbulent-boundary layer observed in the canonical cases considered in Deleon et al. [21] does not readily extend to complex terrain flows. To this end, we consider Bolund Hill as an additional test case to evaluate the box perturbation method and the guidance to determine the perturbation box height for complex-terrain flows.

Simulation of winds around Bolund Hill is a challenging case to simulate, as evidenced by the large variation among the results submitted to the blind evaluation study [35]. The maximum terrain height is approximately 12 m. The hill is much smaller in size than the Askervein Hill. Therefore, we are able to adopt a spatial resolution that is much finer than the spatial resolution used in the Askervein Hill cases. We are able to test different perturbation box heights but they are derived from the same guidance used for Askervein Hill.

DeLeon et al. [24] simulated both the field-scale and the wind-tunnel-scale version of the Bolund Hill case with the current numerical solver using shifted periodic boundary conditions in the streamwise direction. The current investigation is built upon that study and the best solution from that study is used as a reference solution here.

We use a computational domain of 1533 m × 765 m × 127 m in the streamwise, spanwise, and vertical directions. The grid spacings are 2 m × 2 m × 0.5 m, in the x, y, and z directions, respectively. The RANS-LES blending height $h_B$ is set to 5.3 m. The boundary conditions are identical to the Askervein Hill setup, except that the shift for the periodic boundary conditions [29] is set to 150 m, and the incoming mass flow rate is approximately 1.26 $\times {10}^6$ $\mathrm{kg}{\mathrm{s}}^{-1}$. The aerodynamic roughness length $z_0$ is set to $6\times {10}^{-4}$ m.

The perturbation box sizes are given in

Table 2. Parameters of box perturbation method used in flow over the Bolund Hill case. $L_{PB}$, $W_{PB}$, and $H_{PB}$ are the length, width, and height of the perturbation boxes, respectively. Depth is the number of perturbation boxes in the streamwise direction.

Case	${\mathit{H}}_{\mathit{PB}}$ (m)	${\mathit{L}}_{\mathit{PB}}$ (m)	${\mathit{W}}_{\mathit{PB}}$ (m)	PB Depth
Bol1	2	4	4	3
Bol2	17	35	35	3

. Case Bol1 is based on the inner layer height of the Bolund Hill measurement campaign, reported as ≈2 m in Bechmann et al. [35]. Case Bol2 follows the Obukhov length recommendation that we discussed earlier in Section 4.1. As reported in Bechmann et al. [35], one of the criteria to determine if measurements contributed to the neutrally stratified database was to have an Obukhov length greater than 250 m. Here, as an approximate number, we will assume an Obukhov length of 300 m. With that assumption, the one-eighth recommendation would result in only three perturbation boxes in the vertical direction, which is not enough. Thus, we only pursue the one-sixteenth recommendation here, giving an $H_{PB}$ = 17 m. Again, $R{i}_{PB}$ is unchanged and set to 0.042. The working fluid is air at 293 K. We compare our predictions against the field measurements collected along line B in Berg et al. [34] for the 270-degree wind direction. The measurement locations relative to the terrain are shown in Figure 6.

Figure 7 presents the development of turbulence in the incoming flow as a result of buoyancy perturbations. Two perturbation box sizes are considered. The largest of the perturbation boxes is similar in size to the smallest perturbation box used in our Askervein Hill case simulation. The stochastic buoyancy perturbations develop into realistic-looking turbulence after an approximately 1000 m fetch, as shown in Figure 7b. Given that the development region before the hill is less than 1000 m, it is not surprising that the flow has not reached a fully turbulent state at the hill site. A fully turbulent state is eventually reached past the hill. The case with smaller perturbation boxes breaks down into small structures after a short-development fetch as in Figure 7c. The flow maintains its turbulence state for some significant distance before the hill.

Figure 8 shows the fractional speed-up ratios at the four masts, M7, M6, M3, and M8. The wind profiles corresponding to the three cases shown in Figure 7 are compared against each other and the field data. Overall, differences among the three cases are small, but differences among them grow slightly close to the surface. The taller perturbation box (i.e, $H_{PB}=17$ m) results in faster flow close to the surface relative to the short perturbation box. We observed a similar effect in the Askervein Hill case as well (i.e., taller perturbation boxes create faster flows close to the surface). With the exception of predictions at M6, results obtained with the shorter perturbation box agree better with the field data. It is likely that the large eddy viscosity at the escarpment produced by the turbulence model at hand is damping the small-scale turbulence in the case with the short perturbation boxes. Despite the differences among the predictions, the variation is at the same level of uncertainty that exists in the field measurements. Overall, we deem the predictions with the box perturbation method successful despite some dependence on the size of the perturbation boxes.

5. Discussion

Inducing stochastic buoyancy perturbations [18,19] to trigger turbulence in the inflow for large-eddy simulations is a promising concept. In the present study, we applied the box perturbation method [21], which was only tested for canonical fluid dynamics problems at low Reynolds numbers in the original study, to complex terrain flows. A unique aspect of the box perturbation method is that it introduces colored noise as volumetric perturbations in boxes that encompass several computational cells, where a bulk Richardson number is used to scale the amplitude of random source term perturbations in the pseudo-temperature transport equation. We showed that the same bulk Richardson number guidance as proposed in Deleon et al. [21] for the canonical cases of turbulent-channel flow and backward-facing step applies to complex-terrain flows as well, which is encouraging. However, the physics of complex-terrain flows is more complicated than the physics of canonical aerodynamic boundary-layer type flows. Consequently, we found that the guidance to determine the perturbation box sizes, developed for canonical turbulent-boundary layer type flows, does not readily extend to complex-terrain atmospheric flows and has to be revised.

We tested two options to decide on the perturbation box height: as a fraction of the Obukhov length and as the height of the inner layer. The shorter perturbation boxes are considered to develop into turbulence over a shorter distance than their taller counterparts. However, the fine-scale structures, created by the shorter boxes, resulted in low fractional speed-up ratios in the Askervein Hill case. We obtained better results with taller perturbation boxes, but a longer development fetch was needed for those perturbations to break into realistic turbulence. We observed a similar trend in the Bolund Hill case as well. Based on our numerical experiments, we recommend a perturbation box height of about one-eighth to one-sixteenth of the Obukhov length for the complex-terrain flow problem at hand. Our proposition should be tested for unstable and stable stratification as we have only considered neutral stratification in the present study. We also found that the bulk Richardson number guidance (i.e., $R{i}_{PB}=0.042$, as suggested in Deleon et al. [21]) for low Reynolds number canonical fluid dynamics cases works well for the complex terrain flows considered in the present study.

Our current guidance to use the Obukhov length to determine the perturbation box height for complex terrain flows is preliminary and limited to neutral stratification. However, we observe that the effect of perturbation box size on the first-order statistics is not substantial. Detailed information on the statistics of turbulence in the incoming flow and around complex terrain will be useful to develop robust guidance for determining the perturbation box sizes.

6. Conclusions

Accurate wind predictions over complex terrain benefit greatly from temporal and spatial resolutions that are much finer than the resolutions used in meso-scale numerical weather prediction models. To this end, the present micro-scale wind simulation approach is well-suited to accommodate such fine resolutions. The inflow generation technique that is investigated here is computationally inexpensive, and, therefore, ideal for forecasting applications as well.

In the current investigation, we simulated neutrally stratified flows around complex terrain. Buoyancy perturbations due to the turbulent inflow generation technique were coupled to the momentum field within a short buffer zone. Future work should investigate the impact of such perturbations on simulations with both stable and unstable stratification. The inflow generation technique could also benefit from the quantification of turbulent scales generated due to different perturbation box sizes and comparison with experimental data.