Evaluating Terrain as a Turbulence Generation Method

: When driving microscale large-eddy simulations with mesoscale model solutions, turbulence will take space to develop, known as fetch , on the microscale domain. To reduce fetch, it is common to add perturbations near the boundaries to speed up turbulence development. However, when simulating domains over complex terrain, it is possible that the terrain itself can quickly generate turbulence within the boundary layer. It is shown here that rugged terrain is able to generate turbulence without the assistance of a perturbation strategy; however, the levels of turbulence generated are improved when adding perturbations at the inlet. Flow over smoothed, but not ﬂat, terrain fails to generate adequate turbulence throughout the boundary layer in all tests conducted herein. Sensitivities to the strength of the mean wind speed and boundary layer height are investigated and show that higher wind speeds produce turbulence over terrain features that slower wind speeds do not. Further, by increasing the height of the capping inversion, the effectiveness of topography alone to generate turbulence throughout the depth of the boundary is diminished. In all cases, the inclusion of a perturbation strategy improved simulation performance with respect to turbulence development.


Introduction
In the real atmosphere, no matter how far upstream one travels from a region of interest, there is turbulence in the atmospheric boundary layer (ABL), and it propagates and affects the flow downstream. However, when one models atmospheric flow, the computational domain is finite, and ends some distance upstream of the region of interest. This creates a challenge in performing non-idealized, non-periodic, turbulence-resolving, microscale large-eddy simulations (LES) of the ABL in situations such as those with complex terrain [1,2]. At the upstream boundary we have the options to impose: (i) turbulent flow that perfectly matches the flow situation further upstream, (ii) the mean flow only, or (iii) the mean flow plus some perturbation or approximation of turbulence. Inflow type (i) is an unrealistic one in complex situations because it requires a turbulence-resolving precursor simulation of the upstream, which in itself faces the same problem. This leaves the more general inflow types (ii) or (iii), which come with the disadvantage that it takes distance, also known as "fetch", for the ABL turbulence to develop and come to an adequately equilibrium state with no memory of the imperfect upstream boundary conditions. In this work, we examine whether or not terrain itself is an effective means of producing a realistic turbulent ABL within a reasonable fetch using LES with inflow boundary conditions of type (ii) or (iii).
The problem of turbulence generation and its evolution to equilibrium over complex terrain arises in many situations. Examples include LES of wind flow through wind farms or urban areas in complex terrain [3,4]. It is also a pertinent problem in mesoscale-coupled microscale simulations in complex terrain [5]. Mesoscale simulations are performed with numerical weather prediction (NWP) tools resolving regional-scale atmospheric processes over hundreds of kilometers. It would be too computationally expensive to use high enough resolution grids to resolve turbulence over the large domain extents of a typical mesoscale simulation. Thus, mesoscale simulation grid resolutions are typically larger than 1 km horizontally, and turbulence is parameterized by an ABL turbulence model [6]. In a mesoscale-coupled microscale simulation, the microscale domain is often O(1 − 10) kilometers in extent and lies somewhere within the mesoscale domain. Time-varying flow data, such as velocity and temperature, are extracted from the mesoscale model on the microscale boundaries and used as boundary conditions for the microscale simulation. Because the mesoscale model does not resolve ABL turbulence, we are left with the turbulence generation problem discussed above.
It has been shown that the turbulence generation process can be accelerated in idealized simulations over flat domains through the addition of perturbations to non-turbulent inflow boundary data [7][8][9]. Additional studies have been carried out simulating from synoptic scale flow to microscale flow using a domain with multiple nests and initial and boundary conditions provided by a global climate model or reanalysis data set [10,11]. Within these studies, due to the complexities of the flow, it is difficult to assess what process is generating the microscale turbulence. To the best of our knowledge, little has been done to show how turbulence is generated by complex terrain itself and what effect inflow perturbation techniques have to enhance this process. In the work of Mirocha et al. [7], small amplitude, sinusoidal perturbations to the surface topography were applied in order to assess their usefulness as a perturbation technique in the situation in which mesoscale inflow boundary data was driving a flat microscale LES. While these topographic perturbations did decrease the fetch within the microscale domain, it was found that the turbulence spectrum of the u-component wind field, wind speeds, resolved TKE, and friction velocity were poorly predicted. The terrain forcing utilized by Mirocha et al. [7] was idealized and of low ruggedness. Lee et al. [12] studied the ability of buildings to generate turbulence in an LES and find that buildings do not generate sufficient turbulence above the building canopy and recommend the use of a perturbation strategy to accurately simulate turbulent processes in urban simulations. Lee et al. [12] hypothesize that complex terrain will produce a similar result to urban simulations in which the terrain features will be unable to generate adequate turbulence throughout the ABL. In this work, we test this hypothesis by studying flow over realistic terrain over a range of ruggedness to better understand how ruggedness affects efficient turbulence generation. In most studies concerning turbulence formation acceleration over flat terrain, a key metric for method effectiveness is how the required fetch is decreased horizontally (see Mirocha et al. [2,7], Muñoz-Esparza et al. [8], Muñoz-Esparza et al. [13]). With terrain-generated turbulence, though, the vertical growth in the turbulent field must be considered, as turbulence is first generated through interactions with the surface. As non-turbulent flow passes over topographic features, turbulence may become "tripped" depending on how rugged the topographic feature is and the speed of the flow. Vertical growth of the turbulence will also depend on the steepness of the topographic feature (an aspect not considered in the urban simulations conducted within Lee et al. [12]), stability of the atmosphere, and the height and strength of the capping inversion at the top of the boundary layer.
In order to examine the effectiveness of terrain, itself, as a means of generating turbulence in non-periodic, non-idealized LES of the ABL, we use the terrain from a portion of the Rocky Mountains and their foothills as a case study, and we vary four different aspects of the flow: 1.
Terrain Ruggedness: The effect of ruggedness on turbulence generation, the main focus of this work, is studied through varying ruggedness of this terrain by applying different levels of smoothing to it. It is important to also keep in mind that flow from a different direction will interact differently with the terrain and the flow will effectively experience different levels of ruggedness depending on wind direction. However, to reduce overall complexity of this study, we only study winds from a single wind direction.

2.
Inflow Turbulence Seeding: As a control, we use a fully-turbulent, canonical LESgenerated inflow; second, we apply the non-turbulent average of that turbulent inflow; and third, we use that average non-turbulent inflow plus turbulent-triggering perturbations.

3.
Mean Wind Speed: We expect that varying mean wind speed affects turbulence evolution and fetch. For example, we conjecture that it is possible that a case with moderately complex terrain and slow wind speeds will not generate turbulence as effectively as when the wind speeds are higher because wind speed affects the flow separation process.

4.
Boundary Layer Height: We conjecture that it may require more fetch for taller boundary layers to become fully populated with turbulence.
In the following sections, we describe the details of our numerical experiment and then we present and discuss the results, and finally provide conclusions and future directions.

Outline of Numerical Experiment
We view this work as a numerical experiment in which we examine how turbulence develops within the boundary layer simulated in LES of flow over terrain where the domain has distinct inflow and outflow boundaries and, therefore, a fetch region where turbulence must develop and come to a fully developed state. The four parameters we adjust are terrain ruggedness, inflow turbulence seeding type, mean wind speed, and boundary layer height. Ideally, one wishes to minimize the fetch in one of these simulations to reduce computational cost.
To examine the effect of varying these four parameters, we performed a suite of cases described in Table 1. Two terrain ruggednesses are studied, one from the full-resolution shuttle terrain radar mission (SRTM) data, and one with a low-pass filtered version of that same terrain. We use three different types of inflow turbulence seeding: seeding with realistic turbulence from a "precursor LES", no seeding, and seeding through organized random temperature perturbations. We vary wind speed through the strength of the applied background pressure gradient. Capping inversion height is varied simply through the height at which we place the inversion in the inflow boundary conditions and initial conditions for the potential temperature field.
In the Base cases, a baseline boundary layer is produced with an initial boundary layer height of 550 m, and a specified density-normalized driving pressure gradient of 2.2 × 10 −4 m 2 s −2 . In the Fast cases, the driving pressure gradient is doubled in order to produce a faster mean wind. The Tall cases revert back to the original driving pressure gradient but double the initial boundary layer height to 1100 m. For all three cases, a suite of simulations is performed in which the terrain complexity is varied between complex and smoothed, and the inflow is varied between turbulent, non-turbulent, or containing cell perturbations. This results in 3 precursor simulations and 18 terrain simulations in total.

Numerical Method
The simulations herein are conducted using SOWFA (the Simulator For Wind Farm Applications) [14,15], which is built upon the open-source computational fluid dynamics software, OpenFOAM [16,17]. SOWFA is used to perform large-eddy simulations (LES) of atmospheric/wind-energy flows. It solves the filtered incompressible Navier-Stokes equations. To account for density stratification, the vertical momentum equation includes a Boussinesq buoyancy term, which means that a transport equation for filtered virtual potential temperature is also solved. A one-equation subgrid-scale (SGS) turbulent kinetic energy equation is solved, which is used to compute the SGS turbulent stresses and temperature flux [18]. These partial differential equations are discretized using a cell-centered, collocated-variable, finite-volume formulation on unstructured meshes. Although the code can handle arbitrary mesh cell shapes, we always use hexahedral cells. We can add regions of higher resolution through splitting of sets of cells, and the resultant mesh is still a single, un-nested mesh. For flow over terrain, our meshes are terrain following with vertical grid lines near the terrain as orthogonal to the terrain as is reasonably possible. Time-advancement uses a second-order accurate, backward-in-time, implicit scheme. The equation system is solved sequentially using the pressure-implicit splitting operation (PISO) [19,20] but with an outer loop to maintain the second-order accuracy in time. Because the code is designed to handle unstructured meshes, the linear systems solvers are all iterative Krylov-type solvers. For example, the pressure equation is solved using the conjugate-gradient method with a multigrid preconditioner.
All simulations use the same general boundary conditions. The inflow boundary conditions for velocity, potential temperature, and subgrid-scale turbulent kinetic energy are of the Dirichlet type while the outflow conditions are of the Neumann type. To assure global mass flow rate balance, the outflow velocity flux derived from a zero-normal-gradient outflow condition is adjusted to match the specified inlet velocity flux. Velocity and temperature are not specified on the lower surface, but rather Schumann's model for surface shear stress and heat flux are specified. Monin-Obukhov similarity is used to estimate friction velocity, and thus shear stress, at the surface. In flat horizontally homogeneous conditions, we would normally use Monin-Obukhov similarity in a horizontally averaged way, but because we are dealing with complex terrain, we use Monin-Obukhov similarity locally at each surface-adjacent grid cell. The upper surface is treated as a rigid, stress-free lid for velocity, and the gradient of potential temperature is specified to match the initial potential temperature profile gradient at that location. Because this is an incompressible solver and the pressure Poisson equation is derived from the continuity equation (which in finite-volume form is simply a balance of velocity flux over grid cell faces), at any location where velocity flux is specified, no boundary condition for pressure is required. Therefore, in this solver, there is no need to specify boundary conditions on pressure: the upper and lower boundaries specify zero velocity flux, the inlet boundary directly specifies the velocity flux through the Dirichlet specification of velocity, and the outflow boundary indirectly specifies mass flow rate through first a zero-gradient extrapolation from the interior followed by a correction to ensure that global outflow flux equals global inflow flux.

Idealizations
As stated in the introduction, we focus on real terrain from the Rocky Mountains. We consider this real terrain instead of an idealized topography to ensure that the topographic features are realistic and not under-or over-predicting the contributions to turbulence generation by terrain. While we recognize this inhibits the generality of the study, we feel that considering an extremely rugged region such as the Rocky Mountains offers insight into how turbulence is developed over some of the most rugged terrain regions and whether these topographic features are sufficient to generate adequate turbulence throughout the ABL. Our other option was to synthetically generate terrain, but we felt that this had more disadvantages in that there is simply no good replacement for real terrain.
In order to simplify this study and create flow that is predominantly west to east at all heights, Coriolis forcing is not included. By including Coriolis and the subsequent wind veer and flow from the southwest at lower levels, we would have turbulent-deficient areas along the southern edge of the domain that would change in distance along the domain and with height creating an increasingly small analysis domain. Additionally, as noted in Stull [21], the Coriolis term within turbulent kinetic energy equations cannot create turbulence but merely redistributes it. Thus, we do not anticipate the exclusion of Coriolis to significantly change the results other than from the aforementioned issues with its inclusion in the model setup.
Another simplification is the use of a single constant surface roughness. We choose medium roughness height indicative of long grasses, which reflects the real land cover a significant part of our domain. We acknowledge that this idealization will have effects on how turbulence is generated near the surface, but the idealization also allows us to isolate the effects of the terrain itself, which is the main focus of this study. Interesting follow-on work could be to study the effect of surface roughness value, to examine the use of canopy models versus surface roughness alone, and to investigate the application of heterogeneous roughness based on a land-use model.
Following Lee et al. [12], in order to consider only the turbulent contributions from the surface ruggedness, only neutral conditions are considered. While we recognize that neutral conditions are less common over land, convective and stable conditions will add additional sources and sinks of turbulent development that would convolute the analysis of how well rugged terrain is able to produce turbulence within the ABL.
A neutral profile is chosen above the capping inversion (rather than the typical stable lapse rate found in the free atmosphere) due to the excitement of gravity waves within the model when a stable lapse rate was applied [22]. By considering a neutral profile above the capping inversion, no gravity waves are excited and the simulation is able to reach a quasi-steady state. The initial wind speed is set to a constant 10 m s −1 throughout the column.

Inflow Turbulence Precursor Simulations
Because one of our inflow turbulence seeding methods is to perform a precursor LES to generate turbulent inflow, we explain the precursor LES here and its relation to the terrain simulation schematically in Figure 1. A flat precursor simulation with periodic boundary conditions is initially run to develop the turbulent neutral boundary layer over a period of four hours. At the end of the 4-h spin-up period, the precursors are then run for an additional 4.5 h while collecting vertical slices along the y-z plane every 0.5 s at the eastern boundary. These slices are used as the inflow in the turbulent cases. The precursor domain, indicated by a II in Figure 1, is 6 km by 4.16 km in the eastwest and north-south directions, respectively. It is periodic in the lateral directions. The numerical grid consists of four vertical layers with the grid near the surface spanning the lowest 1.5 km vertically with 20 m cubic grid spacing (∆ x,y,z ). Within the successive layers spanning from 1.5 km to 2.02 km, 2.02 km to 6.02 km, and 6.02 km to 9.86 km, ∆ x,y,z increases to 40 m, 80 m, and then 160 m, respectively. The surface sensible heat flux is set to zero to develop a neutral boundary layer. The time step considered within the precursors is 0.5 s. A prescribed surface roughness of 0.15 m, a value typical of long grasses, is used. This value is suitable for the grassland at the origin of the domain, however, it may not be indicative of the true land use over the extent of the domain.
The initial conditions in the precursor are important, as they greatly influence the forward evolution of the flow, especially for potential temperature, and they are indicated in Figure 1 as I. They consist of a potential temperature profile set to 300 K from the surface to the bottom of the capping inversion, followed by a 5 K increase over 200 m, and then a constant 305 K up to 10 km. The initial velocity is set to constant with height.
The turbulent wind and potential temperature field are used as the inflow for the "turbulent" cases outlined in Table 1. The spanwise mean of this precursor generated flow is used as the inflow profile for the "non-turbulent" and "perturbed" cases, the difference between the two being that the "perturbed" case also contains organized random perturbations meant to trip turbulence.

Perturbation Technique
For the "perturbed" inflow cases, some sort of perturbations must be added to the nonturbulent inflow that is derived from the mean of precursor flow. Because several strategies have been developed for adding perturbations to the inflow in order to decrease fetch, this subject merits some discussion. These methods include the cell perturbation strategy [8], adding stochastic turbulence from TurbSim [23], and utilizing a kinematic simulation using discrete Fourier-Gabor modes approach [24]. While all have been shown to be effective, a variation of the cell perturbation strategy is utilized in this study. This strategy has been shown to be computationally efficient and effective in generating turbulence in the neutral stability regime [8,9]. The cell perturbation as described in Muñoz-Esparza et al. [8] consists of adding random perturbations to blocks of cells at a time, specifically, an 8 × 8 block of grid points on the horizontal plane. The basis of this approach is due to the need to compensate for the numerical damping within the solver in which smaller scales of fluctuations would be damped out and unrealized. In Muñoz-Esparza et al. [8], it was shown that this method was very computationally inexpensive as well as the best performer in developing turbulence that (1) compares favorably with the control simulation on a periodic domain, and (2) is generated in the shortest distance (fetch) from the inflow boundary.
Inflow perturbations have been implemented in a fashion similar to that of Muñoz-Esparza et al. [8], but with important differences. Muñoz-Esparza et al. [8] apply rectangular "cells" of temperature perturbation that encompass multiple grid points horizontally. Multiple layers of these perturbation cells are applied extending from the inflow boundary inward into the domain. Muñoz-Esparza et al. [8] and Muñoz-Esparza et al. [9] have also optimized the size and magnitude of their perturbation cells based on local conditions. Our implementation differs first by the fact that we do not apply temperature perturbations on the field interior of the inflow boundaries. Instead, we apply temperature perturbations to the inflow temperature field directly on the boundary. We select rectangular patches of inflow boundary faces on which the inflow temperature is perturbed by the same value, and that perturbation value is held constant for a certain amount of time before being randomly updated. This gives an effect similar to "extruding" cells of uniform temperature perturbation into the interior of the domain. The width and height of these rectangular perturbation patches on the inflow boundary dictates the cross-stream width and height of the "extruded" perturbation cells that enter the domain. The inflow velocity and time between perturbation updates dictates the along-flow length of the "extruded" perturbation cells. We choose this approach because of its relative simplicity for implementation into an arbitrarily unstructured grid code, such as SOWFA. The method of Muñoz-Esparza et al. [8] was applied to WRF-a structured code in which locating perturbation cells only requires specifying ranges of i, j, k grid indices, logical ordering which is not present in unstructured codes.

Terrain Cases
The main focus of this study are the simulations of flow over complex terrain, indicated by III in Figure 1. To assess a more realistic terrain forcing, observed terrain data from over the Rocky Mountains near the Denver metropolitan area is utilized due to the complex features of the mountains to the west of Denver and the more subtle terrain variations to the east. The terrain data set utilized here is the 1 arc second data from the Shuttle Radar Topography Mission (SRTM; Farr et al. [25]) data set. This amounts to approximately a 30 m horizontal grid size for the topographic data. When considering the raw SRTM data, there is a mean decline from west to east across the Rocky Mountains east of the Continental Divide. To remove this mean declination such that the mountain peaks and river valleys are perturbations about zero, a high-pass uniform convolution filter is applied with a window of just over 4 km (135 × ∆ x,y,z ). The resultant SRTM data can be seen in Figure 2a. In order to generate intermediate levels of complexity from the same data set, the SRTM terrain data is filtered using a window size of 2010 m (67 × ∆ x,y,z ; see Figure 2b).
Two further modifications to the terrain are made. At the western edge of the domains, terrain height is set to zero along the north-south direction for the first 800 m and is then gradually transitioned to the value of the complex terrain through linear weighting with distance over two kilometers. This allows the western inflow boundary to have the same shape and size as the precursor simulation. Finally, the north and south edges of the terrain are smoothed to zero at 80 m from the boundary over 500 m such that a periodic boundary can be utilized, which avoids the use of a generalized inflow/outflow side boundary, that well handles the turbulence and lateral flow that forms along that boundary.  Due to the smoothing of the north and south edges of the domain, it would be unfair to expect the non-turbulent simulations to develop adequate turbulence within these areas over the same distance and time as the turbulent and perturbed inflow cases. Thus, in order to compare only the areas of the domain where all simulations are interacting with the complex topography, a subdomain is designated as seen by the dotted black lines in Figure 2a,b; this region will be referred to as the analysis domain for the remainder of this paper. Between the dotted black lines are colored circles representing the locations of virtual towers that store the wind and temperature fields at each model height within the lowest 1.5 km every 1 s within the simulation. The virtual tower locations are colored in groups referred to as Tower Lines which represent different distances from the inlet. The analysis domain spans the lateral extent of the Tower Lines from −1.11 km to 1.11 km in the north-south direction. The origin of this domain is defined as the location at which turbulence generation will be primarily analyzed, thus the locations upstream are all relative to this location that is past the complex terrain features.
The initial fields of the terrain cases are set to match the initial inflow non-turbulent profiles of velocity and temperature, but with local profile stretching to account for changes in column height caused by the terrain. The SGS turbulent kinetic energy is set to 0.1 m 2 /s 2 uniformly.

Precursor
After the initial four hours of spin-up for the three precursor cases, initial transients have passed, the mean velocity profile has good law-of-the-wall agreement, and turbulence extends to the capping inversion. At this point, velocity and potential temperature data are collected from the western domain boundary for an additional 4.5 h. The terrain simulations are then run for 4.5 h with the initial 1.5 h considered model spin-up where initial transients are able to pass through the terrain domain.
The planar-averaged vertical profiles shown in Figure 3 depict the average inflow conditions at the start of data collection. The average temperature profiles (left column) show temperatures of 300 K beneath the capping inversion and 305 K above. The Base (top) and Fast (middle) simulations both show a boundary layer height centered around 550 m with the Tall (bottom) case having the capping inversion centered around 1100 m. As can be seen in the center column, the averaged v-component wind speed is negligible due to the absence of the Coriolis force. However, the average u-component vertical profiles show the typical neutral boundary layer profile shape with a drag-induced shear layer near the surface and a deceleration beneath the capping inversion. Further, the Fast case produces stronger wind speeds throughout the column than the Base and Tall cases. Lastly, the average variance of the u-, v-, and w-component winds is shown in the right column with peaks in variance near the surface and then a decrease toward zero until a small increase by the capping inversion in the u-component variance. Above the capping inversion, the variances all decrease to near-zero, though not exactly zero due to small variations that develop within the free atmosphere. From this time, the boundary layers simulated within the precursor simulations continue to evolve and flow into the terrain simulations discussed below.   Figure 4 shows horizontal cross sections of the u-component wind field at 110 m for the turbulent (a), non-turbulent (b), and perturbed (c) inflow cases at 4.5 h into the complex terrain simulations (total time = 8.5 h). As can be seen near the inlet (west side), the turbulent case shows a heterogeneous turbulent field while the non-turbulent case shows a constant inflow which only becomes turbulent at around x = −11 km. The perturbed case shows a more gradual transition from an unresolved-to resolvedturbulent state. Although the temperature perturbations are applied at the inflow boundary, fluctuations in the horizontal wind field appear after some streamwise evolution of the flow. Once the flow reaches the large complex terrain features, the flow at 110 m is nearly indistinguishable between the three cases. The topography appears to have triggered turbulence development that is independent of upstream flow. However, it is necessary to inspect these simulations with time and height in order to determine exactly how similar they are.   Here, and for the analyses involving time series to follow, the first 1.5 h of the terrain simulation are discarded to allow for the passage of the initial transient. Comparing the time series, all three cases appear similar at each tower line location, except Tower Line 1. At Tower Line 1, only 2 km downwind of the inlet, the fluctuations are overall relatively small; they are largest for the turbulent inflow case and very small for both the cell perturbation and non-turbulent cases as the complex terrain has only just begun. Once the flow reaches Tower Line 5, the ranges are much larger and the three cases appear very similar. It is at this location that the flow has just passed the most complex of the terrain, which explains how the fluctuations are the largest and appear to be very similar due to the terrain determining the turbulence generated. By Tower Line 7, the terrain has begun to flatten out again and, thus, the fluctuations relax back to smaller amplitudes, but remain similar for each case. This is further shown in Figure 6, where, by Tower Line 5, the mean and variance (Figure 6a and b, respectively) have converged to a similar value.   For simplification, we removed the Coriolis acceleration from our simulations and provided a pressure gradient along the flow. This causes the flow to slowly accelerate over time. To compute the mean of this accelerating flow for the calculation variance, skewness, and kurtosis, we apply a double exponentially weighted moving average of the time series based on Croarkin et al. [26] as

Base Case Simulations over Complex Terrain
where with α and γ being constants with values between 0 and 1. The initial values are set to be U t=0 = U t=0 and b t=1 = U 2 − U 1 with the averaging starting after the first time step. Last, the constants α and γ are defined as with ts 1 = 400 s and ts 2 = 950 s. These parameter values were found by trial and error in order to produce a running mean that does not lag behind the signal but is smooth enough such that the computed fluctuations about the mean have the correct magnitude. This technique was justified due to the slow response to changes in the mean flow when using a running average. Using an averaging with too much lag behind the slowly accelerating mean wind speeds yielded abnormally high variances as the entire flow field was seen to be above the running mean. Notice that variance in the non-turbulent and perturbed inflow cases start near zero, then quickly surpass and finally converge to match the values of the turbulent inflow case ( Figure 6). This effect is often seen when turbulence is initialized in an LES before it comes to equilibrium.
Examining the skewness and kurtosis, no striking pattern in convergence is revealed although the difference between the values for the three cases remains small after Tower Line 4 or 5. In general, the values from the perturbed case are closer to those of the turbulent inflow case, specifically after the first few tower lines.
In studies examining methods to perturb non-turbulent inflow to LES, spectral analysis of the wind components at different downstream locations is commonly utilized in order to determine at which point the perturbed flow came close to resolving a similar amount of energy as flow that was fully turbulent at the inlet [2,7,8]. This technique is utilized here where the spectra at each tower along a tower line is calculated and then averaged so that an analysis of the resolved energy as the flow moves downstream can be analyzed. Figure 7 shows the power spectral density of the u-component wind at 110 m averaged along Tower Lines 1, 3, 5, and 7, (Figure 7a-d, respectively). Focusing first on Tower Line 1, it is clear that the non-turbulent inflow case (light blue) is deficient in energy across all scales when compared to the turbulent inflow (dark blue). Adding the cell perturbations improves this deficiency, however, the average spectra shows that the energy is less than that of the turbulent inflow. Somewhat surprisingly, however, by Tower Line 3, the spectral densities from both the perturbation and the non-turbulent inflow runs have improved to the point where they are nearly identical with the exception of an overabundance of energy at the lower frequencies, which agrees with the trends in variance observed in Figure 6. This abundance will be explained in more detail in coming discussions. By Tower Line 5, it is difficult to determine any significant differences between the three spectral densities, and as the complex terrain subsides before Tower Line 7, the energy returns to something closer to that shown in Tower Line 1 for the turbulent inflow.
Thus far, the simulation results have been analyzed strictly from data at 110 m above the surface. While it has been shown that the three simulations appear to converge statistically on the same value at this height, it remains to be seen how well turbulence has been generated throughout the boundary layer. Figure 8 shows vertical cross sections along the x-z plane at y = 10.0 m for the u-component wind speed at 4.5 h into the terrain simulations. The red line indicates the location of the boundary layer height calculated by taking advantage of the fact that this flow is within a neutral regime; the height at which the temperature profile first becomes larger or equal to 302.5 K (the middle of the capping inversion) is set as the boundary layer height, z i . The inflow regions, again, show the largest differences in the flow field where the turbulent inflow (top row) can be clearly distinguished from the non-turbulent (middle row) and cell-perturbation method (bottom row) inflows. Considering the non-turbulent inflow, we see turbulence generation near the surface as the flow begins to encounter topographic features. Once past the major peaks in the complex terrain, the flow fields all exhibit fluctuations in the wind field throughout the boundary layer. This shows us that the turbulence generated by the terrain along this slice does in fact reach the top of the boundary layer by the time the flow reaches the origin (x = 0 km).  To see this more clearly, Figure 9 shows the same vertical cross section but of the turbulent kinetic energy (TKE) field. Here, TKE is calculated at every time step based on variance from the mean described with Equations (1)-(4), and then a running average is performed over the last 3.5 h of simulation time. This is performed in order to produce a clear depiction of the TKE field in each case. The turbulent inflow case shows a clear amount of increased TKE within the boundary layer at the inlet while the non-turbulent and perturbation cases takes some time to develop measurable TKE. The non-turbulent case clearly develops most of its resolved TKE from the surface. Meanwhile, the cell perturbation case appears to develop a substantial amount of TKE throughout the vertical column after only 1 km or so along with the TKE generated at the surface due to interactions with the complex terrain. It was previously noted that within the precursor simulations there were non-zero velocity variances resolved above the capping inversion. This can be seen by the turbulent inflow case above the boundary layer where TKE is larger than that within the non-turbulent and perturbed inflow cases. By the time the flow reaches the origin, the TKE generated within the boundary layer looks very similar from case to case. It should be noted, however, that this is simply one slice through the domain and may not be indicative of the entire simulation domain. A thorough, full-domain analysis of TKE generation will be presented in Section 3.4.

Base Case Simulations over Smoothed Terrain
When flow passes over the smoothed terrain, turbulence development is very different. Figure 10 shows a horizontal cross section through 110 m for the turbulent (top row), nonturbulent (middle row), and perturbation (bottom row) inflow cases at 4.5 h into the terrain simulation. One of the most striking differences between flow over the complex terrain and the smoothed terrain can be seen in the non-turbulent inflow case in which turbulence does not develop fully at 110 m. Spotty patches of turbulent flow appear throughout the domain; however, there is no region of persistent turbulence at this level. The perturbed case, on the other hand, develops a turbulent field at 110 m that looks similar to that of the turbulent inflow case. Inspecting the flow field statistics (Figure 11), the mean of the perturbation case follows closely with the turbulent inflow case. On the other hand, the non-turbulent inflow case has a mean u-component wind speed greater than that of the other cases from Tower Line 3, onward. We suspect that this is due to the lack of turbulent transport between the surface drag and its momentum sink and the flow at 110 m thus preventing a deceleration of the wind field. The average variance of the u-component winds shows that the non-turbulent inflow case at 110 m remains non-turbulent until between Tower Lines 3 and 4. At Tower Line 4, the skewness and kurtosis both spike as the flow oscillates between non-turbulent and turbulent over the time period of interest. This is more clearly seen in Figure 12 where the probability distribution functions (PDFs) of the perturbation of u-component velocity for the turbulent inflow (red) and non-turbulent inflow (yellow) cases are shown. In the non-turbulent case, the PDF appears to be singular at Tower Line 3 (a) but then has a very long negative tail by Tower Line 4 (b). The perturbations were calculated as before by removing the exponentially weighted moving average from the time series for the last three hours of simulation. This tail reduces by Tower Line 5 (c) but the distribution does not recover to what is seen in the turbulent inflow case. From here, the variance of the non-turbulent inflow case ( Figure 11) fails to reach that of the other two cases. Meanwhile, the skewness and kurtosis appear to all converge on a similar value around Tower Line 7, although the scale is distorted due to large fluctuations around Tower Line 4.  Spectral analysis of the tower lines in the smoothed terrain cases show what one would expect from these flow fields ( Figure 13). Near the inflow (Tower Line 1; a), the nonturbulent inflow is largely energy-deficient across all frequencies. As the flow progresses, the perturbed inflow cases come close to matching the turbulent inflow rather quickly at this level. Meanwhile, the non-turbulent inflow case never fully develops and the energy deficiency remains significant across all frequencies.  The non-turbulent inflow case shows the lack of development of TKE in the vertical in Figure 14b where after the first of the larger terrain features (around x = −8 km) TKE near the surface is generated and very slowly propagates vertically. As with the flow over complex terrain, the perturbed inflow and turbulent inflow cases appear to compare well everywhere except near the inflow and above the capping inversion. However, as shown in Figure 10, the turbulence that develops in the non-turbulent inflow simulation is contained to merely the center of the domain. Thus, it is necessary to devise a way to compare these simulations at all locations in order to quantify how well turbulence develops in these domains.

Base Case Simulations
In order to objectively quantify how well turbulence is generated over complex terrain, a specific strategy was developed. First, the capping inversion was estimated using the same strategy as before where z i is assumed to be the height at which the temperature profile first becomes larger or equal to 302.5 K. Using z i , we are able to normalize height (z/z i ) at each location such that the surface is 0.0 and the boundary layer height is 1.0 (non-dimensional). This effectively allows us to compare each point on the turbulent inflow domain with each respective point on both the non-turbulent inflow and perturbed inflow domains. Next, the vertical profile of TKE is interpolated utilizing a cubic spline onto the common vertical coordinate between zero and one as defined above. The now normalized and interpolated fields for both the non-turbulent inflow and perturbation inflow cases are divided by the turbulent inflow case and multiplied by 100.0 to produce a field showing the percentage of TKE that is developed, with respect to the turbulent inflow case, for both the non-turbulent and perturbed inflow cases. Finally, the averaged TKE percentage field is averaged along the y-direction over the analysis domain (as seen in Figure 2) between y = −1110.0 m and y = 1110.0 m to produce a x-z plane showing how well TKE is generated with respect to the turbulent case. The average is only considered over this region along the y-direction, as there are no complex terrain features near the northern and southern edges of the domains, which will unfairly expect the non-turbulent cases to produce turbulence over a flat surface. This procedure produces the fields shown in Figure 15. For reference, Tower Lines 1, 3, 5, and 7 are shown in dark red, yellow, green, and violet, respectively. Additionally, lines that contour the 90th percentile (dark blue) and 110th percentile (red) are added for reference to where the TKE values are 100.0%±10%. Focusing first on Figure 15a, the percentage of TKE produced by the non-turbulent inflow is shown to be nearly zero by the inlet. As the flow starts to interact with the terrain features (as can be seen in Figure  15e) the TKE percentage begins to increase to, and past, 100% very quickly. The growth in TKE in the vertical reaches the boundary layer top (z/z i = 1), but has produced, by and large, more than double the amount of TKE as in the turbulent inflow case. From here, the TKE values begin to relax toward 100% throughout the column and reach values within ±10% of 100% below z/z i = 0.5. When perturbations are added to the inflow (Figure 15c), the TKE percentage remains near-zero near the inlet but quickly increases throughout the column. As flow reaches the complex terrain, TKE again increases to, and past, 100% and propagates vertically. In this case, however, the values do not overshoot 100% by nearly as much as in the non-turbulent inflow case. Further, they relax to within the 20% window much faster and over a considerably larger depth.
As for the cases over smoothed terrain (right column in Figure 15), the non-turbulent case is largely non-turbulent and only begins to reach around 100% in the last 1 km of the domain before overshooting. It can also be seen that TKE levels are far below the turbulent inflow case at all locations above z/z i = 0.3. Concurrently, when perturbations are added to the temperature field (Figure 15d), the TKE percentage improves. An overshoot is still produced, though the amplitude is much lower than the overshoot in the other cases. That said, near the capping inversion the percentage is below 90% throughout the extent of the domain. If the scope is limited to simply the lower half of the boundary layer, the TKE percentage matches closely at the origin and remains above 90%. In both instances, flow over complex and smooth terrain, adding perturbations at the inlet improve the TKE production over the non-turbulent inflow method for this particular boundary layer height and wind speed.

Fast Case Simulations
By increasing the driving pressure gradient magnitude, it is expected that the resultant higher wind speeds will trip turbulence more than with slower wind speeds. Figure 16 compares the instantaneous u-component velocity field at 4.5 h into the terrain simulation for the non-turbulent cases over complex terrain (a) and smoothed terrain (b). As before, turbulence appears to be generated at the center of the domain and through the boundary layer in the complex terrain case, but not fully developed within the non-turbulent inflow case.
To see this in the horizontal plane, Figure 17 shows time-averaged TKE at 4.5 h into the terrain simulations for the turbulent (a), non-turbulent (b), and perturbation (c) inflow cases at z = 100 m. The turbulent inflow case contains a TKE field that is relatively homogeneous at this level. When the inflow is non-turbulent, turbulence appears to develop around the center of the domain (this is just after the crest of the first terrain feature) and then spreads laterally as the flow propagates downstream. By the very end of the simulation domain, the TKE field appears to be developed throughout the lateral extent of the domain. The perturbed inflow case shows streaks near the inlet that are related to the size of the perturbation cell size. After some distance, these streaks develop into a more uniform field and it becomes difficult to distinguish between the perturbation and turbulent cases. These streaks are present in all perturbation simulations, although not shown here.
Using the same method as before, Figure 18 shows the percentage of TKE produced (relative to the turbulent inflow cases) for the non-turbulent (top row) and perturbed (middle row) inflow cases over both complex (left column) and smoothed (right column) terrain. The same general patterns as before are displayed where little to no TKE is produced near the inlet for all cases but turbulence develops from the surface and propagates vertically. The amount of TKE surpasses 110% shortly after becoming turbulent in the non-turbulent fast case (Figure 18a). From there, the overshoot remains much smaller than in the base cases; however, it fails to decrease back to between 90% and 110% within the simulation domain. Considering the perturbed inflow (Figure 18c) case, the overshoot is comparable to that of the non-turbulent case, though the TKE percentage decreases to between 90% and 110% nearly throughout the column by the time the flow reaches the x = 0, our location of interest.
As expected for faster mean boundary layer flow, the turbulence generated on average in the smooth terrain cases (Figure 18b,d) is improved (when compared to similar panels in Figure 15). Further, by adding perturbations to the inflow, the TKE percentage is closer to 100% and exists throughout a deeper layer within the boundary layer. Generally, the turbulence is higher and closer to 100%, within the bottom half of the boundary layer in this case. It remains to be seen if the turbulence produced over the smoothed terrain simulation can fill the vertical extent of the domain to reasonable levels. Thus, it is to be expected that when the capping inversion height is doubled, it will be even more difficult for turbulence to fill the vertical depth of the boundary layer.

Tall Case Simulations
In these cases, the original driving pressure gradient is restored but the height of the capping inversion is doubled to z i = 1100.0 m. Figure 19 shows vertical slices of TKE through y= 10.0 m for the turbulent (a), non-turbulent (b), and perturbed (c) inflow cases over smoothed terrain with the boundary layer height denoted by a red line for reference. By raising the boundary layer height, the vertical depth over which TKE develops in the non-turbulent inflow simulations is reduced, relative to z i . In this case, the flow appears to remain almost laminar near the top of the boundary layer as opposed to in Figure 14 where TKE reaches the boundary layer top, albeit at much lower levels than the respective turbulent inflow case. In the perturbed case, turbulence fills the depth of the boundary layer but at lower levels than in the turbulent inflow case.
These patterns can more readily be realized when examining Figure 20. Here, even in the complex terrain cases (left column), the average amount of TKE (relative to the turbulent inflow case) remains low and does not fill the vertical extent of the boundary layer. For both complex terrain cases, TKE reaches above 90% by the time the flow reaches the origin only in the lower half of the boundary layer. By adding cell perturbations, the TKE values are overall improved toward 100% throughout the domain. Overshooting 100% is minimal in these simulations, which suggests that by increasing z i , there is less energy buildup when the flow first becomes turbulent. When considering smooth terrain (right column), there appears to be improvement by adding cell perturbations; however, the overall amount of resolved TKE decreases when compared to simulations with a smaller initial z i .

Summary and Conclusions
We studied flow over complex terrain in which the inflow is either turbulent, nonturbulent, or contains perturbations to the temperature field throughout the boundary layer. The turbulent inflow conditions represent the control, where resolved, LES-generated turbulence is applied at the inflow boundary. Non-turbulent inflow is meant to represent inflow from a mesoscale model in which turbulence is accounted for within the average profile, but not explicitly resolved. Finally, the perturbed inflow cases are included in order to assess the impact of perturbing the non-turbulent inflow before the flow reaches the complex terrain. The complexity of the terrain is also varied between rugged and smooth in order to identify how well rugged topographic features initiate turbulence. Last, the height of the boundary layer and the strength of the driving pressure gradient force are sequentially doubled to assess the sensitivity of these factors to turbulence development.
It was shown in the Base case simulations that flow over complex topography generates and augments turbulence in a way that after a short distance, the average wind fields, variance, kurtosis, skewness, and spectra all match well at lower levels. However, once the full height and analysis domain are considered, the levels of TKE are found to spike initially throughout the column and then relax into values within 100 ± 10% of the turbulent inflow values by the origin. By smoothing the terrain, the non-turbulent inflow cases fail to become turbulent in the majority of the domain, under-resolve the u-component spectra, and produce negligible amounts of TKE. Skewness and kurtosis show large spikes around Tower Line 4 due to this region containing areas of turbulent flow and laminar flow, which skew the overall statistics. The statistics begin to recover but fail to ever match the turbulent and perturbation cases. Adding perturbations improves the results drastically throughout the boundary layer and produces results that are much more in line with the turbulent inflow case.
Increasing the mean wind speed was hypothesized to increase the production of turbulence by causing it to "trip" over features that slower flow would not. Additionally, by increasing the mean wind speeds, the vertical growth of turbulence was expected to shift further downstream as the faster wind speeds advect the flow further east. These simulations show that the vertical propagation of turbulence in the non-turbulent and perturbed cases is indeed advected further downstream; however, the normalized TKE profile comparisons show better agreement with the turbulent inflow simulation than in the slower-flow Base case. Without perturbations, the overshoot that is realized is just over 110% but broadly dispersed across the boundary layer. With perturbations, the overshoot is again reduced as in the Base case, and a large portion of the boundary layer is within the 100 ± 10% range. Further, the increase in wind speed over the smooth terrain does appear to improve the turbulence generation over the Base case, albeit not to the point that non-turbulent inflow comes close to matching that of the turbulent inflow. Once again, the addition of perturbations appears to improve the results, although, in the Fast case, the area within the 100 ± 10% range is located strictly within the bottom half of the boundary layer.
When the boundary layer height is increased, the relative height of the terrain to the boundary layer effectively reduces. Due to this, the turbulence that is generated by topography will have a longer vertical path to travel in order to fill the boundary layer. This holds true within the Tall simulations where TKE takes a longer distance to fill the vertical extent of the boundary layer in the complex terrain simulations. In both the non-turbulent and turbulent cases, the TKE generated falls within the 100 ± 10% range within the lower half of the boundary layer by the time the flow reaches the downstream point of interest. Improvement in decreased fetch and an overall deeper boundary layer is achieved through the inclusion of cell perturbations, though even this addition cannot aid in producing enough turbulence throughout the entire depth of the boundary layer. When considering the smoothed terrain, the layer of developed turbulence within the 100 ± 10% range remains very small and adding perturbations helps only slightly. Tuning the perturbation magnitude and size could remedy this problem, though.
Overall, it is shown here that relying on terrain alone to initiate turbulence does not produce desirable amounts of turbulence in any of the simulations considered. By including perturbations at the inlet, the amount of variance, resolved spectra, and magnitude and location of TKE are improved across the board. Considering real data cases in which wind direction changes with time and height, it is important to note that the resulting LES could have periods in which winds are coming from areas of complex terrain, and others in which the winds are coming from smoother terrain resulting in insufficient turbulence generation. Thus, the use of the cell perturbation method is important in that it improves performance for both cases and the resulting turbulence field can be more robust.
We also show that in addition to monitoring the downwind evolution of turbulence, it is also important to understand how turbulence fills in the boundary layer vertically. This study only includes simulations with the cell temperature perturbation strategy and sensitivity to the magnitude of the perturbation is left for future work. The focus of this work is not necessarily to find the perfect perturbation strategy and setup for this case, but rather to investigate if perturbation methods are necessary, and if so, effective in complex terrain.
Within this idealized study, a constant roughness length was used over the entirety of the domain in order to isolate impacts of terrain. It is expected that changing the roughness length would alter the resulting turbulent flow field, albeit at a smaller scale than that of the complex terrain features. Additional studies could focus on the effects of roughness length and the use of canopy models in generating turbulence. Based on the results shown here, it is suggested that further simulations considering non-turbulent inflow into a microscale domain should include a perturbation strategy at the inlet to improve turbulence generation, even in very complex terrain. Additionally, the perturbation strategy may have to be tuned in order to achieve the desired amount of turbulence based on the mean wind speed, stability, and boundary layer depth [9,27].
Author Contributions: P.H. performed the simulations conducted herein with the guidance of M.C. Analysis code was written by P.H. All authors have read and agreed to the published version of the manuscript.
Funding: This work was authored by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding provided by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes. The research was performed using computational resources sponsored by the Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory.

Data Availability Statement:
The analysis code and data that support the findings of this study will be made openly available through Zenodo.

Conflicts of Interest:
The authors declare no conflicts of interest.