Long-range transport of terrain-induced turbulence from high-resolution numerical simulations

Over complex terrain, an important question is how various topographic features may generate or alter wind turbulence and how far the influence can be extended downstream. Current measurement technology limits the capability in providing a long-range snapshot of turbulence as atmospheric eddies travel over terrain, interact with each other, change their productive and dissipative properties, and are then observed tens of kilometers downstream of their source. In this study, we investigate through high-resolution numerical simulations the atmospheric transport of terraingenerated turbulence in an atmosphere that is neutrally stratified. The simulations are two-dimensional with an isotropic spatial resolution of 15 m and run to a quasi-steady state. They are designed in such a way to allow an examination of the effects of a bell-shaped experimental hill with varying height and aspect ratio on turbulence properties generated by another hill 20 km upstream. Averaged fields of the turbulent kinetic energy (TKE) imply that terrain could have a large influence on velocity perturbations at least 30 H (H is the terrain height) upstream and downstream of the terrain, with the largest effect happening in the area of the largest pressure perturbations. The results also show that downstream of the terrain the TKE fields are sensitive to the terrain’s aspect ratio with larger enhancement in turbulence by higher aspect ratio, while upstream there is a suppression of turbulence that does not appear to be sensitive to the terrain aspect ratio. Instantaneous vorticity fields shows very detailed flow structures that resemble a multitude of eddy scales dynamically interacting Correspondence to: M. Katurji (marwan.katurji@gmail.com) while shearing oppositely paired vortices. The knowledge of the turbulence production and modifications by topography from these high-resolution simulations can be helpful in understanding long-range terrain-induced turbulence and improving turbulence parameterizations used in lower resolution weather prediction models.


Introduction
Through 2-D high-resolution computational fluid dynamic simulations, we examine the transport of turbulence generated by topographic features in the atmospheric boundary layer.Atmospheric flows have high Reynolds number (Re) on the order of 10 9 and the turbulent eddies generated typically range in scale from a few metres (surface eddies) to boundary layer height (a couple of kilometres).
Atmospheric turbulence is very complex in nature and is often a combination of many spatial and time scale phenomena, (Lin, 2007;Pielke, 2002;Stull, 1988;Orlanski, 1975), ranging from daily cycles of tens of kilometers to sub-second events at the sub-millimeter scale.The interest in atmospheric turbulence comes from its importance in mixing phenomena and high frequency energy transport.Applications can include air pollution, aviation safety, wind energy, which have been greatly aided by the development of supercomputing capabilities (Anderson et al., 2009).In the boundary layer, turbulent kinetic energy (TKE) per unit mass (m 2 s −2 ) as shown below statistically characterizes turbulence, which is the average sum of the squared perturbation velocities.
Topography produces TKE by mechanical shear, which represents the transformation of the mean kinetic energy of the flow to the turbulent kinetic energy through turbulent momentum flux interaction with the mean wind shear.
Resolving the smallest-scale atmospheric boundary layer flow poses a huge computational challenge, requiring simplifications of the physical system to reduce the solution so as not to cover the full turbulence scale.An example is the ongoing development of large eddy simulations (LES).LES was first pioneered by Lilly (1967), in which large-scale turbulence is fully resolved and the sub-grid scale effect is parameterized, based on the energy cascading principle.Details of the numerical process can be found in Mason (1994) and Sagaut (2001).Current application trends in applied LES can be found in cross-disciplinary studies with various success rates.Bouffanais (2010) points out the importance of the mechanisms and difficulties of the cross-disciplinary application of LES, one of which is the choice of the sub-grid parameterization aimed at better representing flow features at smaller scales, in particular near the surface.Theoretically, running LES code with very fine resolution makes the solution invariant to the sub-grid scale parameterization, but improvements can be made by using a hybrid model of more than one sub-grid parameterization scheme or running an ensemble of LES codes (Iizuka and Kondo, 2004;Fedorovich et al., 2004;Beare et al., 2006).
A practical use of LES is its application in the wind energy industry.Turbulence information is crucial for characterizing the cyclic structural loads that ultimately lead to fatigue and failure of the turbine structure.LES provides more accurate dynamic turbulence information as opposed to the traditional methods used in the European Wind Atlas system (Troen and Petersen, 1989) that were ultimately designed for flat terrain, attached flows, and relied on linear models from Jackson and Hunt (1975) and Mason and Sykes (1979).Many studies reported an underestimation of flow separation, turbulence intensities and gust measurements in the linear model solver, thus stressing the importance of using nonlinear models in the wind resource prediction (Palma et al., 2008;Castro et al., 2003).The importance of LES was highlighted 10 years ago in the review of Wood (2000).In more recent times, the application of LES in the wind energy industry has taken a big stride forward, and its economic application in a full-scale resource assessment has been examined by Tamura (2008), in which the author pointed out that computational costs are not the hindrance but it is the parallelization algorithm that will make a good practical use of LES.Another application of using LES output in the wind energy industry was reported by Kelly et al. (2004) in which the authors used LES output for atmospheric waves as input for the wind turbine structural design.The above are all examples of the practical use of LES in wind energy applications.
Improvements of LES capabilities as a result of rapid advancements in hardware technology and parallel computing algorithms permit a detailed examination of atmospheric tur-bulence phenomena that are difficult to observe, but are very important to understand for practical applications.Acquiring real measurements of turbulent eddies, or coherent turbulent structures in the atmospheric boundary layer can be a very difficult task, especially if the dynamic information is required.The limitation is the absence of an instrument that can take a long range snapshot of turbulence as the eddies travel over terrain, interact with each other, change in their productive and dissipative properties, and are then observed tens of kilometres downstream of their source.
To our knowledge, little is known about long-range transports of turbulence, as indicated by turbulent kinetic energy (TKE), from high-resolution model simulations, as is presented in this paper.Previous work has considered the details of rotor dynamics in relation with stability and upstream wind profiles on the immediate lee side of mountains, but did not focus on the long-range transport of TKE.Hertenstein and Kuettner (2005) investigated rotor types produced on the lee side of an idealized mountain using two-dimensional (2-D) LES simulations with horizontal grid spacing of 150 m.Doyle and Durran (2007) extended rotor dynamic analysis to include 2-D and 3-D simulations; their horizontal resolution was 60 m and also in a stable atmosphere.Moreover, Harman and Finnigan (2010) in their numerical investigation of turbulence within canopies concluded that to properly assess local turbulence it is important to have knowledge of the terrain in a larger area, which can have influence on the local turbulence measurements.
This work focuses on the evolution of TKE field downstream of the turbulent rotor formed by a hill.The turbulent intensity of the rotors as they travel downstream over terrain with varying steepness or aspect ratio (AR) is also examined.The diagnosis of the TKE at this high-resolution is crucial in quantifying the velocity intermittency for practical applications such as wind energy.The TKE derived from the high-resolution simulation can also provide a good framework for the development of high-order turbulence closure schemes, which are required for lower-resolution numerical weather prediction models.The simulations carried out in this work are 2-D turbulence resolving simulations.2-D and 3-D turbulence resolving simulations are very different in nature (physically and numerically).Practically, there might be no 2-D turbulence in the physical world, and it might only exist when geometric length scales permit the dominance of one dimension over the other.One fundamental difference between 3-D and 2-D turbulence is the vortex stretching in the 3rd dimension, which is not present in 2-D.The second fundamental difference is that 3-D turbulence is characterized by downward energy cascade into smaller scale, while in 2-D there exists an inverse cascade from smaller to larger energetic scales (Davidson, 2007).The above mentioned differences are to be taken into consideration when interpreting the results of the simulations.

Methodology
The research approach in this work is purely numerical, and this is due to the technical and economic difficulty in measuring TKE at such a large scale.The model used was the LES version of the Atmospheric Regional Prediction System (ARPS) (Xue et al., 2000(Xue et al., , 2001) ) that was specifically designed to simulate small-scale atmospheric phenomena like thunderstorms and tornadoes.ARPS is a non-hydrostatic, compressible atmospheric flow model based on the Navier-Stokes equations and uses a terrain following coordinate system.ARPS was validated successfully in numerous studies (e.g.Chow et al., 2006;Weigel et al., 2006), and also cross-compared with other LES codes for different atmospheric stability conditions (Fedorovich et al., 2004;Beare et al., 2006).ARPS was also used in some idealized highresolution experiments to investigate turbulent structures over forested hills (Dupont et al., 2008;Fesquet et al., 2009).
The ARPS simulations are 2-D with four grid points in the cross wind direction (north-south).The solution is carried out at one central point and the rest are for boundary conditions.The 2-D simulations resolve turbulent eddies in the wake of hills at 15 m spatial resolution with north-south periodic boundary conditions.The choice of periodic boundary conditions allow the solution to represent eddies as infinitely long circulation rolls in the north-south direction.This type of setup is very typical of a 2-D simulation.An eddy resolving model, or LES model, to have any practical significance or application should be truly 3-D.Only then it could be compared to DNS or observations.By definition, only 3-D simualtions can be classified as LES (Moeng et al., 2004).This work's simulations are classified as 2-D turbulence resolving simulations, with subgrid scale (SGS) parameterizations for the unresolved turbulence.Moreover, Moeng et al. (2004) show that, in their 2-D turbulence resolving simulations, that turbulent statistics of 2-D simulations can compare well with 3-D-LES when proper tuning and selection of the SGS model is performed.
The basic idea behind the current design of numerical simulations is to be able to simulate high wind speed scenario that is able to generate consistent turbulent eddies or rotors in the wake of a terrain feature.The perturbations produced by the terrain feature travel downstream over a relatively long distance of approximately 20 km before they encounter an experimental series of terrain configurations with varying aspect ratios.The turbulent perturbations, assessed by analyzing the resolved TKE, are compared for zones 15 km upstream and downstream of the experimental terrain to study the long-range transport and modification of the TKE field.

Model simulation setup
The simulations are 2-D with an isotropic spatial resolution of 15 m and run to a quasi-steady state in a neutral atmospheric background, and with horizontal and vertical number of grid points of 4099 and 353, respectively.The simulations are designed in such a way to allow an examination of the effects of a bell-shaped hill with varying steepness or aspect ratio on turbulence properties generated by another hill 20 km upstream.The inflow (westerly) velocity is assumed to be uniform at 15 ms −1 and constant in time.The model domain is 61.5 km in horizontal dimension and 5.3 km in vertical dimension.Bottom boundary condition is non-slip rigid wall, while the top is a zero gradient condition with a Rayleigh damping sponge layer in the upper third of the vertical domain.The lateral inflow (west) and outflow (east) boundary conditions are set to zero gradient while periodic boundary conditions are used at the north and south boundaries, The surface heat exchange processes were turned off so as to consider a pure mechanical turbulent production at this stage of the current research.For the sub-grid turbulence parameterizations the 1.5-order TKE turbulence closure scheme is selected.In this scheme an additional prognostic equation of TKE is solved, and the horizontal and vertical turbulent mixing coefficients are represented as a function of TKE and length scales as outlined in Deardorff (1980).The computational time step is 0.1 s and all simulations are run to a steady state that was usually reached after 5 h of simulation time.The analysis results correspond to half hour averages after steady state is reached, which typically takes about 5 days real time on a 32-processor computer cluster.

Terrain setup
The model domain contains two bell-shaped hills separated (peak to peak) by 20 km.The topography of both hills can be described by the following function, which is commonly called in the mathematical community by the "Witch of Agnesi": where H max is the maximum height of the hill, X c is the center position or the coordinate corresponding to H max , and, X hw is the hill half width.The first hill (referred to as turbulence generator) is responsible for producing the background turbulence that is then modified by the downstream or the second hill (referred to as experimental hill).Two sets of experiments are carried out and constitute in total 7 simulations.Set (I) consists of 5 simulations with an eddy generator of fixed height (H = 500 m) and aspect ratio (AR = 0.25) and 4 experimental hills of the same height (500 m) but with different aspect ratios.The aspect ratio is calculated based on the ratio of the height of the bell-shaped hill to the hill half width (H max /X hw ).So technically the slope of the experimental hill changes from one experiment to another, and is based on the steepest face facing the incoming flow.One of the experiments in set (I) is the base simulation, and this only contains the eddy generator without     any downwind experimental terrain.The base simulation is necessary to provide the background flow for estimating the modifications by the experimental hill of the turbulence produced by the eddy generator and advected downstream to the experimental hill.A summary of the terrain setup for set (I) is shown in Fig. 1a; the steepness increases gradually for experiments t02E01, t02E02, t02E03 and t02E04 from 14, 18, 26, to 45 degrees.The second set (II) of experiments, shown in Fig. 1b, consists of two simulations in which the sensitivity of turbulence fields to the height of the terrain generator and to the height of the experimental hills is examined.Also shown in Fig. 2 are selected zones from which the half-hour averaged fields are extracted.Zone A, located in between the eddy generator and the experimental hill, highlights the upstream effects of the experimental hill on the flow field, while zone B is downstream of the experimental hill.Also shown, in Fig. 2, are locations of several vertical probes from which vertical profiles are extracted.

Results and analysis
In this section we first present results for the average fields produced by the simulations to show agreement with previous published work, and subsequently present the TKE fields.The results are presented in separate sections for set (I) and set (II) of the experimental runs, and all correspond to 30-min averages.The resolved TKE is calculated from horizontal and vertical velocity perturbations sampled at a frequency of 1 Hz.The velocity perturbations are derived by subtracting the one-minute mean velocity from the total resolved velocity, which are then used in deriving the oneminute TKE.The reason behind choosing a one-minute averaging period for the TKE calculation is based on the turnover cycle of the largest eddies present in the domain.This is deduced from time series spectral analysis of the vertical velocity that indicated a frequency peak at around 7 to 11 min.So choosing a sampling interval of one-minute should be sufficient to statistically represent the average.An inspection of the parameterized sub-grid scale near surface (15 m above ground level) TKE, deduced from the turbulence closure scheme, reveals that the model is resolving 97 % of the TKE and only 3 % is parameterized at the current spatial resolution.For this reason the TKE used in the analysis will be considered to be the resolved TKE.

Vorticity and rotor dynamics
Figure 3 presents three instantaneous fields at three consecutive times (2 min apart) of the vorticity extracted from one of the experiments of set (I).The extraction location is between distances 40H to 90H (H is the maximum height of experimental hill (Fig. 2) covering a distance of 60 km downstream of the experimental hill.Immediately downstream of the hills, flow separation introduces mechanical shear leading to the formation of eddies that rotate clockwise (negative vorticity).These eddies typically scale to the height of the experimental hill.Between 35 and 38 km small rotor begins to form as the incoming flow over the hill shears and separates from the lee side zone characterizing the circulation zone in the wake of the hill.Smaller rotors of similar vorticity signs join together as shown in the time evolution of rotors (r1), ( r2) and (r3) of Fig. 3a.
As eddies propagate downstream, they become larger because they gain energy from the background flow.The centre of the large rotors is around the peak of the experimental hill and their sizes extend to around 1000 m in the vertical and about 1500 m in the horizontal.The large eddies leave behind a trail of smaller eddies that rotate in the same direction.The energetic large scale rotors rotate in a clockwise direction with a clear shear perimeter on the outside of the eddy characterized by the red color; this is evident in the trailing edge of rotor (R3) in Fig. 3a.These large rotors, as in (R3) and (R4), join together where rotor (R3) bottom half section tilts forward as shown in Fig. 3b to join with the slower moving rotor (R4) as shown in Fig. 3c.
Of interest is the shedding of smaller counter-rotating eddies from larger eddies.These interesting features are probably produced by the shear instability at the outer edge of the larger eddies; smaller rotors detach on the upper edge to produce rotors (v2) and (v1) in Fig. 3a, which rotate in the opposite direction and are around 300 m in the vertical and 500 m in the horizontal.They tend to shear off and get embedded in the background flow, moving faster downstream than the larger eddies that are attached to the ground leaving a trail of negative vorticity in their wake that appears to interact with the new rotors coming from upstream.It is worth noting that the formation of the above mentioned rotors represent a time snapshot in the simulation, and not a spatial or temporal average.Vorticity was derived from sequenced time snapshots during the simulation and not averaged.
Another interesting feature in the vorticity diagrams is the evidence of very weak background waves (please refer to the area enclosed by the square in Fig. 3) dispersed throughout the domain trailing behind eddies.These numerical imprints resemble wind generated sand cusps on the beach.The rotors forming on the upper edge of the circulation zone in the lee side of the experimental hill, rotors (r1), (r2) and (r3), are similar to those previously documented in numerical simulations of lee waves and rotors for stable atmospheres (Doyle et al., 2009).The rotors in the previous studies had similar dynamical properties as the ones in our simulations, except that they were forming further downstream of the hill as an interaction of strong down slope winds with the lifted flow due to the mountain wave in the stable atmosphere, this interaction causes mechanical shearing of subrotors that travel along the crest of the wave.In our case the atmosphere is neutral and hence the absence of the lee side mountain wave causes an earlier formation of the rotors that shed off the relatively slower winds in the lee side of the hill.Also, lidar observations during the T-REX field campaign (De Wekker et al., 2009) showed evidence of sub-rotor activity on the leading edge of the main rotor that is embedded in the large-scale wave and rotates counter clockwise with the www.atmos-chem-phys.net/11/11793/2011/flow.In our simulations the sub-rotors forming downstream of the experimental hill, are rotating counter clockwise, but additional smaller rotors are formed on the upper edge of the main rotor and then drawn downwards to the surface by the pressure field created by the large rotor advection.

Mean wind speeds
The introduction of topography produces areas of positive and negative pressure perturbation, Fig. 4. The perturbations are responsible for accelerating or decelerating the flow field depending on direction of the pressure gradient.The highest negative pressure perturbation is found in the lee side of both the eddy generator and the experimental hill where horizontal mass convergence is observed.Downstream of the eddy generator, Fig. 4a, c, the pressure perturbation field is quite complex.Below 200 m, there are surface-based patches of positive and negative pressure perturbations.These roughly correspond to the flow feature of the boundary layer eddies as they roll over the surface (i.e.causing a surface return flow).As will be discussed later, these effects manifest as reduced surface winds and increased wind shear.In Fig. 4c, the introduced experimental hill causes an increase in the positive pressure perturbation downstream of the eddy generator with the strongest effect in the 5.5 km mark in zone A, Fig. 4c.Increase in the upstream penetration of positive pressure perturbations at higher altitudes in zone A also occurs.
The mean horizontal U-velocity for the base simulation (t02) and the ensemble average of the 4 experiments (t02E01, t02E02, t02E03 and t02E04) is shown in Fig. 5a, b and Fig. 5c, d, respectively, and the difference between the base simulation and ensemble average results can be found in Fig. 5e, f.Vertical wind shear of U-velocity extends to a height of 565 m for the base simulation in both zones A and B in Fig. 5a, b.The introduction of the experimental hills thickens the layer affected by horizontal wind shear to 752 m in zone A and B in Fig. 5c, d.This thickening increases the gradient of the horizontal wind shear immediately above the return flow in Fig. 5a.Horizontal wind shear of U-velocity is depicted in the wind speed acceleration zone above the return flow (blue color) in Fig. 5a.
Another important feature is the separation zone in the lee side of the terrain hills.As expected from the high wind speed flow regime over a sloped terrain, the flow detaches from its streamlined regime and creates a recirculation zone in the lee side characterized by negative pressure perturbations and negative pressure gradient force as explained earlier.This zone extends, for the base simulation in Fig. 5a, to 6.5 km or 13H downstream of the eddy generator.The presence of the hill downstream of the eddy generator causes the separation zone to shrink towards the hill to 4.4 km or 8.8H as shown in Fig. 5e.Breuer et al. (2009) showed in a series of direct numerical simulations (DNS) and LES of flows over periodic hills that the recirculation zone was around 5H downstream of the hill and pushed back as Reynolds number increased.Their simulations were carried out for hill separation distance of 9H with a maximum Re of 10595, but in our simulations the hill separation distance was 40H with Re on the order of 108.The return flow intensifies in this case by 4 ms −1 (Fig. 5e) at the 2.5 km mark, which equates to a relative increase of 57 %.This reveals a direct relation between the intensity of the surface return flow and the pressure gradient flow in the lee side of the hill, as was previously pointed out in the 2-D simulations by Doyle and Durran (2002).A general decrease in mean velocity above the surface layer was induced throughout zone A, Fig. 5e, from 8.5 to 17.5 km between the surface to 1000 m height level.
In zone B, Fig. 5f, the experimental hill caused an overall decrease in wind speed, but the interesting features were the surface patches of increased wind speed that were about 3 km apart and correlate to the surface pressure perturbations discussed earlier in this section.
The influence of the hill (experimental or eddy generator) height on the mean flow is revealed in Fig. 6 that show the results from set (II) experiments.In the case of the 1000 m experimental hill, zone A exhibits a larger decrease in the mean wind speed upstream of the experimental hill (Fig. 6a), as compared to the simulation in set (I) shown in Fig. 5c.The major modification that a 1000 m high terrain had on the flow field is depicted in the lee side of the hill in Fig. 6b, c. Between 5.5 and 8.5 km (Fig. 6c), and between 6.5 and 10.5 km for the steeper hill (Fig. 6b), there is an area of high wind speeds that extended from the surface up to about 1000 m.This region was also characterized by very low vertical shear of U-velocity.This phenomenon is produced for the case of www.atmos-chem-phys.net/11/11793/2011/      the 1000 m hill and indicates a feature that is missing in set (I).The pressure perturbations in Fig. 7 indicate a high negative pressure perturbation zone between 2.5 and 4 km.This causes subsidence of high momentum flow from aloft suggested by the mean vertical velocity in Fig. 8, thus enhancing the wind speed in the area between 5.5 and 8.5 km.These simulations highlight the importance of the pressure gradient produced between terrain features and suggest that this forcing can create large horizontal wind speed shear that should be taken into consideration for many practical reasons.

TKE and velocity perturbations
Turbulence generated by the eddy generator extended horizontally a distance of 8.5 km or 17H and vertically a distance of 1000 m or 2H .For zone A (Fig. 9a) and zone B (Fig. 9b) of the baseline simulation, there is a clear interference pattern that is depicted by the periodic high and low values of TKE, which suggests that as the eddies move from left to right and fluctuate during the travel, they leave behind some of their trailing influences that either suppresses or enhances the TKE of the approaching eddy.An examination of the interference patterns reveals that the spatial period (or horizontal wavelength) of the oscillation patterns is similar  to the height of the hill, between 500 m and 600 m.There also appears to be a generation of TKE at the surface every 3 km, which correlates with the surface pressure perturbation wave discussed earlier.The waves that are evident in the horizontal velocity and TKE diagrams (Fig. 9) are not atmospheric waves resulting from the restoring force of gravity due to atmospheric thermal stratification or stability.The simulations were carried out under an isentropic potential temperature background without surface energy balance parameterizations, therefore no moisture or heat exchange between the surface and the atmosphere, thus maintaining the neutral stratification throughout the simulation.The temporal averaging of the TKE field reveals a wavelike character, which only represents mean spatial distributions in the domain with increased or decreased magnitude of TKE.
The addition of the experimental terrain, 20 km downstream of the eddy generator, induces considerable modifications to the TKE field in both zones A and B, as shown  To quantify any cumulative change in the TKE field within the entire zones of A and B, a turbulence intensification factor (TIF) is calculated for each of the two zones.The TIF is the zonal or spatial sum of the TKE difference between the ensemble average of the 4 simulations in experimental set (I) and the base simulation (t02) normalized with the spatial mean of the wind speed.So it represents a non-dimensional parameter that indicates the changes in turbulence intensity averaged spatially and temporally.TIF is plotted separately for zones A and B versus the experimental hill aspect ratio as shown in Fig. 10.The TIF for zone A has a negative sign, indicating a reduction in the turbulence level upstream of the experimental hills, while TIF is positive for zone B and values are higher in zone B than zone A, indicating a large increase in turbulence downstream of the experimental hills.The results also reveal that the influence on turbulence in zone A is not sensitive to the aspect ratios, but downstream in zone B, the results vary with aspect ratios and the influence appears to be non-linear.
Vertical and horizontal velocity perturbations are extracted from 6 probes in zone A and B. The extraction locations are given in Fig. 2, and the vertical profiles of the velocity perturbations are given in Fig. 11.Over flat terrain with neutral atmosphere, the vertical velocity perturbation usually decreases as function of increasing distances from the terrain, since there is no thermal stratification to create vertical mo-tion other than the topographic lifting force.The vertical velocity perturbations of the base simulation (t02), as shown in the gray dotted lines in Fig. 11, reduce in intensity with distance downstream of the eddy generator.This is not similar for vertical velocity perturbation of the ensemble average of the four subsequent experiments with the presence of experimental hill (black dotted lines in Fig. 11), where there is a consistent peak in intensity near the height of the hill (500 m) throughout zones A and B. Below half of the hill height (200 m), TKE is primarily generated by the horizontal velocity perturbations, while above that vertical velocity perturbations dominate TKE production with a peak between 400 and 500 m.These results show the importance of vertical velocity perturbations induced by terrain even in a neutral atmosphere.
The 1000 m hill produces more intense turbulence (almost twice of the 500 m eddy generator), and the strong turbulence is seen in the entire zone A of experiment t03E01 in Fig. 12c.The TKE field is also characterized by periodic zones of high and low intensities, but this time they scale to the new hill height and are spatially separated by around 1000 m.In zone A (Fig. 12c) of experiment t03E01, a gap in turbulence intensity is created in the area of high wind speeds and low wind shear as indicated in the corresponding U-velocity (Fig. 6c).This feature produced a low frequency wave (wave length of ∼15 km that is equal to the separation distance between the eddy generator and experimental hill) in the TKE field in which a high frequency signal (wave length of 1000 m) is embedded.In zone B (Fig. 12d), the influence of the 500 m high experimental hill is to damp out the turbulence generated by the eddy generator, and break the long wavelength into smaller periods that now scale to the 500 m high terrain.
As for the experimental run t02E07, shown in Fig. 12a, b, the influence of the 1000 m high experimental hill appears to significantly reduce the upstream TKE, while in zone B, as shown in Fig. 12b, the TKE appears to be the superposition of the low frequency (high wave length, 15 km) coming from the 1000 m high terrain with the high frequency (low wave length, 500 m) TKE produced by the 500 m high eddy generator.The embodiment of high frequency TKE structures within a larger scale (low frequency) wave implies a longrange transport mechanism of turbulence via low frequency waves.

Conclusions
In this paper we have presented 2-D high-resolution numericalsimulations for atmospheric boundary layer flow in an attempt to highlight the importance of long-range transport and modification of the average and perturbation velocity fields by topography.The domain consists of two hills: an eddy generator and an experimental hill 15 km downstream and the atmosphere is assumed to be neutral with a constant wind speed of 15 ms −1 .
Only slight modification to the mean wind speed was recorded as a result of different hill aspect ratios of the experimental hill with the 500 m fixed height, while the largest influence was recorded for the intensification of the return flow 15 km upstream of the experimental hill.More prominent modifications in the mean wind speed were produced when the height of the experimental hill was doubled from 500 to 1000 m, in which the average wind speed in the region upstream of the experimental hill was significantly reduced and relatively strong horizontal wind shear was recorded downstream.The induced changes in the mean flow were directly correlated to the pressure gradient force introduced by the experimental hills, which had a long-range (15 km) effect on the velocities.There was also evidence of a shallow surface pressure perturbation wave that created distinct and consistent zones of accelerated wind speeds; this wave had a wavelength of around 3 km or 6H and was depicted for the set of simulations with 500 m high hills.Instantaneous vorticity fields, in the neutral atmospheric setting considered in the simulations, revealed the formation of small rotors mechanically shedding from the upper edges of the large boundary layer rotors with rotating direction opposite to that of the main rotor, thus showing a different formation concept to that of the rotors in stable atmospheres.
The averaged fields of the turbulent kinetic energy (TKE) implied that terrain can have a large influence on velocity perturbations at least 15 km upstream and downstream of the terrain, with the largest effect happening in the area of the largest pressure perturbations.Spatially and temporally integrated TKE showed that the results downstream of the terrain are sensitive to its aspect ratio with an enhancement in turbulence by higher aspect ratio (taller hill), while upstream there was a suppression of turbulence with little sensitivity to the terrain aspect ratio.Further experiments are required to develop the relationship between the terrain shape and the long-range turbulence field.The stronger pressure perturbations induced by the 1000 m high hill produced a TKE wave on the lee side of the hill with a wave length of ∼15 km.As a result, there was a large suppression of turbulence below its crest and high surface turbulence as its trough reached the surface.Eddies appeared to follow the path of the strongest wind speed avoiding the stagnation region, thus creating an arch resemblance in their trajectory (Fig. 13).Within the 15 km TKE wave there were higher frequency turbulence structures that scaled with the terrain.In one of the presented simulations, a 500 m hill acted as a turbulence suppressor for the TKE background field produced by a 1000 m hill upstream.The TKE field produced by a 1000 m hill and interacting with the TKE field from a 500 m hill shows a solution that suggests a superposition of two wave solutions, but this idea will require more simulations and analysis before can be used to develop a model for long-range TKE transport.
The knowledge of the turbulence production and modifications by topography from these high-resolution simulations can be helpful in improving turbulence parameterizations used in lower resolution weather prediction models.Moreover, the understanding of the TKE long-range transport can prove to be very helpful in practical problems like in the initial phases of wind farm prospecting in complex terrain or atmospheric pollution dispersion problems.But more importantly, the work can have implications towards analysis of wave modeling or multiple-scale perturbations to mathematically represent the transport of the high-frequency vortical structures within the low-frequency waves, which can be the key to understanding long-range terrain-induced turbulence.

Figure 1 .
Figure 1.Experimental terrain setup for set (I) and set (II) of simulations.

Figure 2 .
Figure 2. Dimensions and locations of the analysis zones A and B (inside the shaded boxes), and the probe locations for vertical profiles.

Figure 2 .
Figure 2. Dimensions and locations of the analysis zones A and B (inside the shaded boxes), and the probe locations for vertical profiles.

Fig. 2 .
Fig. 2. Dimensions and locations of the analysis zones A and B (inside the shaded boxes), and the probe locations for vertical profiles.

3 Fig. 3 .
Figure 3. Y-component vorticity (s -1 ) extracted from locations 40H to 90H (see Fig. 2) extending 25 km downstream from the center of the experimental hill.Red color indicates counter-clockwise eddy rotation, while blue color indicates clockwise eddy rotation.Note that the vertical and horizontal axes scales are different.The black square points to evidence of very weak background waves dispersed throughout the domain trailing behind eddies.

3 Fig. 4 .
Figure 4. 30-minute average pressure perturbations (Pa) for the base simulation (t02) in (a) and (b), and for the ensemble average of the 4 set (I) experiments in (c) and (d).Zone A corresponds to the left column while zone B corresponds to the right column (refer to Fig. 2).

3 Fig. 5 .
Figure 5. 30-minute average u-velocity (ms -1 ) for the base simulation (t02) in (a) and (b), and for the ensemble average of the 4 set (I) experiments in (c) and (d).The relative change of the ensemble average and the base simulation is presented in (e) and (f).Zone A corresponds to the left column while zone B corresponds to the right column (refer to Fig. 2).
Figure 6.30-minute average u-velocity (ms -1 ) for set (II) of experimental runs t02E07 (a, b), and t03E01 (c, d).Zone A corresponds to the left column while zone B corresponds to the right column (refer to Fig. 2).
Figure 6.30-minute average u-velocity (ms -1 ) for set (II) of experimental runs t02E07 (a, b), and t03E01 (c, d).Zone A corresponds to the left column while zone B corresponds to the right column (refer to Fig. 2).
30-minute average u-velocity (ms -1 ) for set (II) of experimental runs t02E07 b), and t03E01 (c, d).Zone A corresponds to the left column while zone B rresponds to the right column (refer to Fig. 2).ure 7. 30-minute average pressure rturbation (Pa) for set (II) experimental t03E01 zone A.

Fig. 9 .
Figure 9. 30-minute average TKE (m 2 s -2 ) for the base simulation (t02) in (a) and (b), and for the ensemble average of the 4 set (I) experiments in (c) and (d).The relative change of the ensemble average and the base simulation is presented in (e) and (f).Zone A corresponds to the left column while zone B corresponds to the right column (refer to Fig.2).

Fig. 10 .
Fig. 10.Turbulence intensification factor (TIF) plotted against experimental aspect ratio for the 4 runs of set (I) for zones A and B.

Figure 11 .Fig. 11 .
Figure 11.Horizontal (u) and vertical (w) velocity perturbations for the base simulation (black) and the ensemble average of the 4 experiments (grey).(a, b, c) correspond to probe locations 1, 2 and 3, and (d, e, f) for probe locations 5, 6 and 7. See Fig. 2 for probe locations.

Figure 12 Fig. 13 .
Figure 12. 30-minute average TKE (m 2 s -2 ) for set (II) of experimental runs t02E07 (a, b), and t03E01 (c, d).Zone A corresponds to the left column while zone B corresponds to the right column (refer to Fig. 2).