Simulation of stratified flows over a ridge using a lattice Boltzmann model

A three-dimensional thermal lattice Boltzmann model (TLBM) using multi-relaxation time method was used to simulate stratified atmospheric flows over a ridge. The main objective was to study the efficacy of this method for turbulent flows in the atmospheric boundary layer, complex terrain flows in particular. The simulation results were compared with results obtained using a traditional finite difference method based on the Navier–Stokes equations and with previous laboratory results on stably stratified flows over an isolated ridge. The initial density profile is neutral stratification in the boundary layer, topped with a stable cap and stable stratification aloft. The TLBM simulations produced waves, rotors, and hydraulic jumps in the lee side of the ridge for stably stratified flows, depending on the governing stability parameters. The Smagorinsky turbulence parameterization produced typical turbulence spectra for the velocity components at the lee side of the ridge, and the turbulent flow characteristics of varied stratifications were also analyzed. The comparison of TLBM simulations with other numerical simulations and laboratory studies indicated that TLBM is a viable method for numerical modeling of stratified atmospheric flows. To our knowledge, this is the first TLBM simulation of stratified atmospheric flow over a ridge. The details of the TLBM, its implementation of complex boundaries and the subgrid turbulence parameterizations used in this study are also described in this article.


Introduction
The Lattice Boltzmann method (LBM) has been developed in recent years as an effective computational fluid dynamics (CFD) tool for fluid flow modeling. Unlike traditional CFD methods, such as finite-volume or finite-element methods, which discretize the Navier-Stokes (NS) equations in a macroscopic continuum, the LBM is a mesoscopic kinetic method based on statistical mechanics [1][2][3][4]. The LBM solves the Boltzmann equation for a particle at each grid point by performing collision and propagation calculations of the particle's probability distribution function (PDF) over a discrete lattice mesh with certain fixed directions. Using the principles of statistical mechanics, the macroscopic flow properties such as mass, velocity, and momentum fluxes are computed from the statistical moments of the particle's PDF. There are several advantages of the LBM over traditional CFD (finite-volume and finite-difference) methods that solve NS equations, including the simplicities in code development, in implementation of complex boundary conditions, and intrinsic parallelism [1][2][3][4]. The simplest LBM was originated from lattice gas automata (LGA) [5] to remedy shortcomings of the LGA. The usage of the Bhatnagar-Cross-Krook (BGK) [6] relaxation in the collision term makes the transition from LGA to LBM possible. Further development of the LBM using BGK (LBGK) single relaxation [7][8][9] eliminated statistical noise, and created a viable method for simulating fluid flow. Soon after, the LBM was directly derived from the Boltzmann equation by discretization in both time and phase space [10]. There are many publications that detail further developments of LBM such as the reduction of compressibility effects [11][12][13][14], the implementation of boundary conditions [15][16][17][18], turbulent flow simulations [19][20][21][22], energy and mass transport, multiphase flow [23][24][25][26], and in parallel computation using graphical processing unit (GPU) [27,28]. The potential for drastic speed up using GPU is a very attractive point for the large-eddy simulation of a large domain [29] with a huge number of grid points. A detailed review of the literature, including historical and recent developments in LBM, can be found in [1][2][3][4], with abundant citations.
In spite of many successful simulations of laminar flow problems using LBGK, there are two inherent defects in the LBGK method: the numerical instability in high Reynold number flows and inaccurate boundary conditions [30]. The problem is rooted in the single relaxation time of LBGK. To remedy this shortcoming, d'Humières developed the multiple-relaxationtime (MRT) LBM [31] which has been demonstrated [30,[32][33][34] to be superior over LBGK in numerical stability and accuracy. The MRT-LBM is more stable than the LBGK because different relaxation times can be individually tuned to achieve 'optimal' numerical stability [30,32]. There are many successful simulations using the MRT-LBM, such as turbulent flows [35][36][37][38] and thermal driven convective flows [37][38][39][40]. The apparent advantages of the MRT-LBM provided an incentive for the exploratory simulation study of this paper.
In this study, we investigate the suitability of thermal LBM (TLBM) using the MRT method for modeling stratified flow in complex terrain, using a simple flow configuration. Specifically, we focus on stratified flow past a long ridge, with approaching flow perpendicular to it. To our knowledge, MRT-LBM numerical modeling of atmospheric lee-waves has not been conducted. In this study several aspects of the implementation of MRT-LBM for the problem in hand require careful attention, such as the coupling of momentum with temperature stratification, buoyancy wave generation, treatments of arbitrary curved 3D boundaries, and a turbulent flow parametrization.
Atmospheric lee-waves and rotors over a long mountain ridge have been active research topics for many decades. The first analytical study using linear solutions was derived for a bell-shaped ridge by Queney [41]. Numerous studies on this problem have been carried out in analytical solutions [42][43][44][45][46], both based on linear small-amplitude [41] and nonlinear finiteamplitude formulations of NS equations [43]. The numerical models based on NS equations, including non-linear effects were carried out in computer simulations [47][48][49][50]. Laboratory experiments are also abound [51][52][53]. The mountain waves and associated phenomena, such as rotors, are challenging research problems because of the complexity of its dynamics and problem configurations. Buoyancy waves, rotors, and hydraulic jumps are not only dependent on large scale weather, but also on the characteristics of the topography. Stratified flow over a realistic mountain ridge was selected for this preliminary study hoping that the TLBM technique, a novel tool in atmospheric flow research, will garner increased popularity within the community. In this study, we simulated stratified flow over a bell-shaped hill where the stratification is neutral in the boundary layer, topped with a stable cap and stably stratification aloft. This case has been the subject of numerical simulations [48] and laboratory modeling [53].
Organization of this paper is as follows: The second section describes the thermal MRT-LBM in detail, various boundary conditions of MRT-LBM, and the turbulent parameterization used in this study. The third section presents results of the MRT-LBM simulations of stratified flow over a simple idealized ridge/hill. The last section includes a summary and conclusions.

Lattice Boltzmann model-the MRT method
Before description of the LBM modeling, it is worthwhile to discuss the scaling non-dimensional form NS equations because it is related to our flow variable transformation, and thermal force coupling in the LBM framework. For laboratory and numerical simulations of fluid flows, a set of similarity parameters must be satisfied to represent the real flow conditions. The following normalization can be applied [54] to the thermally coupled incompressible NS equations in Cartesian coordinate with Boussinesq approximation for atmospheric flow, without radiative and latent heat flux terms. Taken the original Boussinesq approximation of the NS equations, the velocity (U i ), length (X i ), time ( ), pressure (P), and temperature (θ) can be normalized by a set of corresponding reference scales u 0 , l 0 , t 0 = l 0 ∕u 0 , 1 − 0 as The resulting non-dimensional NS system is written in tensor form where i3 is the Kronecher delta, St, R e , R a , and P r are the dimensionless parameters which are defined as the Strouhal, Reynolds, Rayleigh, and Prandtl numbers respectively: where 0 is the reference air density; p and ′ are the pressure and potential temperature deviations from the reference state; and, g, , , are the gravitational acceleration, kinematic molecular viscosity and thermal diffusivity, and the atmospheric expansion coefficient. In our study of the atmospheric boundary layer, the St number is not important because it is a dimensionless number for the oscillating of the flow. The other similarity numbers R e , R a , P r must be the same in the simulations as in the real atmospheric boundary layer flows. In this simulation study, the reference scale parameters for velocity (u 0 ), time (t 0 ), and length (l 0 ) are used to transform the model simulated flow variables to the atmospheric wind variables. We focus on 3D atmospheric flows in which the following descriptions of MRT-LBM are 3D in a Cartesian coordinate frame. It should be noted that original LBM is weakly compressible [9,11,12]. However, the compressibility is reduced when using the MRT-LBM [33], therefore the Boussinesq approximation is valid.

Coupled MRT-LBM for stratified flow
In this simulation, we used a 3D thermal MRT model with nineteen lattice directions (D3Q19). The discrete velocity is that of D3Q19 model as follows [31] (Fig. 1): . Fig. 1 The geometry of the D3Q19 lattice with lattice vector c α defined in Eqs. (6) and (7). The lattice is 3D and has 19 discrete collision directions 1 3 The MRT-TLBM equation with α discrete velocities in 19 directions that is equivalent to NS Eq. (2) and (3) can be written as follows with the following equations with MRT equations [31,32]: The symbol in bold-face is defined as a column vector with a number of components in lattice direction α. f ⃗ x, t is the particle velocity PDF vector at a spatial location ⃗ x at a time t; and m and m eq are the velocity moment vector and its equilibrium function vector, respectively. F is an external force vector in which the buoyancy force is coupled. All column vectors in Eq. (7) are expressed as follows, in which superscript T represents the transpose of a vector: M is the transformation matrix that maps the PDF vector f in the velocity space to the vector m in the moment space (see "Appendix" for the matrix). i.e. m = Mf , f = M −1 m. S is a diagonal relaxation matrix, and I is an identity matrix. It is noted here that coupling of the buoyancy force follows Guo et al. [4,55] scheme which satisfies the mass and momentum conservation constraints. The moment vector components can be interpreted as the hydrodynamic quantities and their fluxes [31,32]. The moment vector is arranged in the following order and each component is related to a physical property of the fluid: where ′ is the density fluctuation, e is related to energy, is related to energy squared, j x , j y , and j z are the components of flow momentum, j x , j y , j z = 0 ⃗ u, q x , q y , q z are the third order moments corresponding to energy fluxes in x, y and z directions respectively, p xy , p yz and p zx are related to the components of the symmetric and traceless strain rate tensor, xx and ww are fourth order moments, m x , m y , m z are the third order moment [32]. The orthogonal transformation matrix M is derived using the Gram-Schmidt orthogonalization procedure [31,32] and listed in the Appendix. Since M is an orthogonal matrix, its inverse M −1 can be easily determined by multiplying by its transpose.
There are four conserved quantities in this MRT model-the density ′ , and three components of momentum j x , j y and j z . By using the following decomposition of density, the compressibility effect of the model [33,38] and round-off error is reduced, The equilibrium functions for the other non-conserved moments are given as follows [31,32]: The parameters in the Eq. (11) are chosen to optimize the linear stability of the MRT model [30]: = xx = 0 and j = − 475 63 . With the equilibrium moments above, the speed of sound c s , in the unit of lattice speed, c ≡ dx dt = 1 , is The diagonal relaxation matrix S for the collision process in the moment space is where the relaxation rates of the non-conserved moments are in the range of (0, 2) and are specifically related to the shear viscosity ( ) and bulk viscosity (ζ) [30,32] as, In this simulation, we use the following relaxation parameters for optimized stability [30,32]: As before in (3), the thermal and momentum equation coupling in the NS equations is carried out using a buoyancy forcing term. In the MRT-LBM system, the coupling of the thermal and momentum equation uses an external buoyancy force F as in Eq. (7). The buoyancy force F b has to be computed as a non-dimensional force term in Eq. (7). The similarity numbers, R a , Re, and Pr, computed from the real atmospheric conditions, are used for the buoyancy force F b . F b is then transformed to F . using the following relationship [4,56,57] to satisfy the hydrodynamic mass and momentum conservation: Consequently, this term is added when computing the wind velocity [4,56,57]: 3.

3
The temperature deviation ′ was comped with a MRT model using D3Q7 lattice (see Fig. 1 for α = 0,…,6 lattice setup). The evolution of the probability distribution function h for temperature in the MRT method is [58,59] where c has only seven discrete velocities. The transformation matrix N T is defined by [58] and listed in the "Appendix".
The n and n eq are the thermal moment vector and its equilibrium vector, respectively. The equilibrium moments for temperature are [58,59]: where a is a constant related to thermal diffusivity = √ 3(4 + a)∕60 . To avoid numerical instability of the D3Q7 model, a < 1 must be maintained [37,39], where u, v, and w are the flow velocity components in Eq. (19). The diagonal relaxation matrix R is expressed as [58,60] where relaxation rate k is determined by the heat diffusivity [37,39], , and the Once the shear viscosity is determined, the heat diffusivity is determined by the Prandtl number. In this study, Pr was set to 0.7. The above values for relaxation rates were adopted from previous studies [37,39]. The temperature ′ is computed by [37,39,60]

Boundary conditions
Since the model's accuracy is dependent on the boundary conditions (BC), they are a critical component of numerical modeling. The BCs of the PDF in TLBM are specified to correspond to the macroscopic BC of the fluid, i.e.⃗ u, p ′ , ′ . We applied four types of boundary conditions, (i) the macroscopic no-slip BC, which is equivalent to bounce back BC in the LBM; (ii) the velocity inflow BC; (iii) outflow BC; and (iv) free-slip BC at the top of the domain. The implementation of the boundary conditions is shown in Fig. 2. The prescribed inflow boundary condition uses observed values of wind profiles and in this case is a constant wind velocity. Details on the open boundary condition and free-slip boundary conditions are discussed in the literature [1,4,[13][14][15][16][17]25]. Here we list the free-slip, periodic, and open boundary conditions used in this study for the velocity MRT model, the boundary condition are in similar form for the temperature component MRT (see Fig. 1 for reference of lattice directions).
u are the known wind velocity at west and east boundaries. The arbitrary curved boundary condition will be described in detail since it is critical for this study due to the complexity of the ground surface in atmospheric flows. We have extended Guo's [17] curved boundary condition to arbitrary 3D terrain. First the three dimensional terrain surface is triangulated into many small triangles by drawing southwest to northeast diagonal lines for each ground surface square. Each triangle is an approximation of the ground surface. The terrain surface is preprocessed and all the spatial coordinates for all the triangles are stored, and the intersection point coordinates for each lattice lines are also stored. At this point, all the lattice line intersections with the ground surface is treated as the same in Guo's curved boundary algorithm [17]. (25) Top free-slip boundary conditions, The viscous no-slip surface BC is realized by a simple bouncing back of the PDF in the lattice directions. On a flat surface, the second order accuracy of the no-slip BC is achieved by arranging the physical boundary at the half height of the first computation grid. The specification of an arbitrary curved no-slip boundary condition is much more complex, which, for solid surfaces, has been developed by several researchers [17,61,62,18]. In this research, we use the Guo et al. [17] BC for a curved boundary due to its more robust features that can be applied in both velocity and temperature. Yan and Zu [25] applied this curved BC for a coupled flow and heat transfer problem that obtained good results. The Guo et al. [17] curved BC is based on an extrapolation method and is of second order accuracy. In this method, the PDF for velocity or temperature is divided into equilibrium and non-equilibrium parts. The equilibrium PDF at the solid node (Fig. 3, w node) right below the corresponding boundary node is approximated by extrapolation of the corresponding fluid nodes (f, and ff nodes) with a known boundary point (Fig. 3, b node) wind velocity ( ⃗ u b = 0 non-moving boundary). The nonequilibrium PDFs are extrapolated using the non-equilibrium PDF from the corresponding neighbor nodes in the fluid (Fig. 3, f, and ff nodes). The final assembled PDF for the boundary wall node (w) uses the LBM Eqs. (5) and (15) for wind and temperature respectively.
The above description is put into mathematical equations using Fig. 3 as a reference. The fraction of intersected distance is defined by The extrapolation of the wind speed and temperature at the boundary node (w) are determined by and The equilibrium PDFs for the boundary node (w) can then be computed using the following equations, for velocity (29) and temperature (30), respectively, where c is the lattice vector defined in (4) and is the weight factors in all lattice directions defined in (16).
The non-equilibrium part of the PDFs for velocity and temperature at the boundary nodes are extrapolated by It should be noted that the formulations for different are used to avoid numerical instability. Finally the equilibrium and non-equilibrium part of PDFs at the boundary node are used for BCs for the PDF, f + ⃗ x w , t for the post-collision step, i.e.,

Subgrid turbulence parameterization
In this research, the Lilly-Smagorinsky [63,65] turbulence model is used. It uses the eddy viscosity, which is related to the local strain rate tensor, in analogy to the molecular momentum diffusion determined by molecular viscosity. The equation for total momentum diffusion at subgrid scales is determined by a total effective viscosity ( total ) that includes molecular viscosity ( ) and turbulent viscosity ( t ): where C is the Smagorinsky constant, which is taken as 0.1, Δ the filter size that is taken as grid size, and S is the magnitude of the filtered strain rate tensor The strain rate tensor components can be computed directly from the non-equilibrium moments [28]: Note that in implementing the Smagorinsky subgrid model in LBM, the strain rate tensor is computed from the moments directly, instead of from the flow velocities which saves computational time. For the subgrid heat fluxes, we use the turbulent Prandtl number (Pr = 0.7) to compute the turbulent heat diffusivity (κ t ) The relaxation time for the thermal LBM Eq. (11) is adjusted to account for the subgrid turbulence.

Results
As in any other atmospheric numerical model, the TLBM simulation results and the actual flow have to go through a scale transformation. In addition to the domain of the real flow, which has to be scaled to the model, the initial and boundary conditions also have to be scaled accordingly. After the model simulation, the results are converted back to the real atmospheric flow domains using the scaling relationships in Eq. (1). It should be noted that we have tested different physical grid sizes in the range of 1-50 m for various flows such as 2 m resolution in urban and 50 m in a larger complex terrain domains. The results (not shown here) of varied grid resolutions from 1 to 50 m (much smaller than molecular speed of atmosphere) does not affect the results.

Lee flow regimes with an inversion cap and neutral boundary layer
The simple ridge set-up is shown in Fig. 2 (300 m height) with a horizontal half-width of 1.5 km. The boundary conditions are also shown in the Fig. 2. A neutrally stratified layer is setup between the ground and inversion layer, and an inversion cap is located between the boundary layer and the stably stratified upper layer. This is a more realistic condition for smaller scale lee-wave development that often happens in nature, where a neutral boundary layer develops due to surface layer mixing.
Using the laboratory experiments of Knigge et al. [53] as guidance, the following parameters were used for the simulation. The first non-dimensional parameter is the obstacle height, defined by H/Z i , where H is the ridge height and the Z i is the neutral layer height. The second non-dimensional parameter is the inversion Froude number, defined by where U is mean wind speed, Δ is the potential temperature jump at the inversion cap at the top of the neutral layer of height Z i , and 0 is a background temperature taken as zero. (37) S xx = − 1 38 0 dt s 1 m 1 + 19s 9 m 9 , S yy,zz = − 1 76 0 dt 2s 1 m 1 − 19s 9 (m 9 ∓ 3m 11 , S xy,yz,zx = − 3s 9 2 0 dt m 13,14,15 .
The reference temperature is taken as 290 K. The third non-dimensional parameter is the ridge Froude number defined by where H is the ridge height, and N is the Brunt Väisälä frequency defined as follows: The domain setup for this simulation is as follows: 501, 21, and 101 grid points in x, y, and z, respectively, and dx = dy = dz = 20 m, and H = 300 m. In this setup, the scaling parameters are l 0 = 20, t 0 = 0.2 and u 0 = 100. The inflow wind speed was scaled with the same ratio of 100, providing a model inflow velocity of u = 0.1c, where c is the lattice velocity. The Mach numbers are 0.03 for the physical flow and 0.173 for the numerical model. This value of Mach number has been tested by Wang et al. [39] in their study and produced correct results. In order to compare with previous numerical modeling [48] and laboratory results [53], we use a similar idealized ridge height h(x) of cosine shape: where K = 2π/L, and L = 3 km is the width of the ridge (Fig. 2).
To demonstrate the capability of the model in simulating flow regimes under different atmospheric stability parameters and inversion heights, we first show three representative simulations with varied physical parameters that resulted in three different flow regimes. The physical similarity parameters were defined for the three cases of simulations as listed in Table 1 and were selected based on the numerical simulations of Vosper [48] and laboratory water channel experiments of Knigge et al. [53]. These three cases are for different flow phenomena, including the appearance of lee-waves in which no reversal u component is observed in the wake area, rotors in which reversal u component is observed in the wake area, and hydraulic jumps in which there is a steep gradient of u velocity (10 ms −1 per 100 m increments) at the temperature cap and large turbulent intensity (> 1) in   Fig. 4. The top panel x-z cross section (Fig. 4a) showed that transient, low amplitude lee-waves developed in the simulation. However, there was no reversal of flow under the wave train. The rotor case is shown in the middle panel (Fig. 4b), in which there are greater wave amplitudes than the top panel. Also, a rotor and near-ground recirculation is developed at the lee side of the ridge. The hydraulic jump case is shown in the bottom panel (Fig. 4c). In this case there is much greater wave amplitude than in the top and middle panels, and upper level flows of high velocity penetrated the wake region near the ground surface. The vortices elongated in the vertical direction, and a vortex was present near the ridge peak. In all cases, the nonlinear lee-flow phenomena were frequently transient. For the x-y cross sections, the flow fields are consistent to their flow regimes with very small turbulence variations in y direction because the terrain are a 2D ridge without variation in y direction. However, the 3D dimension simulation is critically necessary for producing turbulent flow field simulation since the turbulent vorticity transform is a 3D process. The flow regimes were also established in a range of physical similarity parameters as other researchers have done [48,53]. In this analysis, we fixed Fr, U/NL, and H/L, to values similar to those used by Vosper [53]. By doing simulations with varying Fi and H/z i , an array of flow patterns were derived (Fig. 5). Our TLBM simulations produced flow patterns of lee-waves, rotors, and hydraulic jumps, similar to those observed by Vosper [48] and Knigge et al. [53] over an equivalent range of F i and H/Z i . The results indicated that a high inversion height Z i with moderate F i promote the generation of lee-waves, while a medium to low inversion height z i with medium to high F i promote rotors, and a moderate to low inversion height coupled with low F i were derived (Fig. 5). Some differences in the flow patterns were observed in the region of F i > 0.8 and H/z i > 0.5 compared to the results of [48], in that our simulation always produced some lee-waves vis-à-vis no lee-waves in [48]. The TLBM simulations are closer to the results observed in the water channel by Knigge et al. [53], wherein the laboratory results also show lee-waves in the same physical parameter range.  [53]. The red colored symbols represent the results from Vosper's numerical simulation [48] In this simulation study, we kept the Reynolds number in the range of 15,000 to 24,000 in concurrence with the laboratory study of Knigge et al. [53]. It is important to point out the differences in the laboratory setup [53] and in this numerical simulation. The laboratory test [53] used a towing ridge model, in which only the model surface is a no-slip boundary condition, and other the channel flow is at rest before the towing model reaches the location. On the contrary, in this numerical simulation the entire ground surface in the domain is a no-slip boundary. The consequence is that the laboratory flows at Re = 20,000 were mostly laminar, except for sporadic turbulence at the immediate lee side of the ridge. The literature indicated [66], for this type of the towing model laboratory setup, the critical Re for the transition from laminar to turbulence is 50,000. Observations [67] and numerical simulations [68] indicated that boundary layer separation and turbulence play very important roles in the formation of lee-wave and rotors. The 2D version of our model simulations (not shown here) did not produce much turbulence except in the wake region because the 2D model suppressed turbulence production. Obviously, a detailed 3D simulation study is needed for investigating complexities such as the interaction of boundary layer turbulence and lee-waves. The waves and rotors produced by buoyancy acted as a transport mechanism which increased the dimensions of the turbulent wake region in which it was taller and longer at the lee side. The wave and turbulence interaction also probably enhanced the turbulent region in the lee side boundary layer.
To compare our simulations to the laboratory results of [53], we have ran our 3D model simulation in the laminar regime by setting the Re = 2000 while preserving keeping the other parameters constant. The results are shown in Fig. 6. For the large lee side flow structures, the laminar simulation results are in better agreement with the laboratory test which is not surprising as the numerical simulation is in a laminar flow regime. In this laminar regime, only shallow laminar boundary layers are developed at the frontal and lee sides of the ridge. The features of the lee waves, rotors and hydraulic jump structures are stronger than those for the turbulent regime. This can be explained by the vigorous turbulent diffusion process in the turbulent boundary layer developed in the large Reynold number.
One point should be noted here that a damping layer near the top boundary was not used. The damping layer is often used in lee wave numerical modeling using NS equations to avoid wave reflection by the top boundary which does not appear to be a problem in this MRT based LBM modeling simulation.

Turbulent characteristic of the simulation results
Turbulent flows in the atmosphere are always an important subject to discuss since the turbulence dominates the momentum, heat, and scalar transport in the atmospheric boundary layer. In this simulation, we used the Lily-Smagorinsky turbulence model for the subgrid scale. To evaluate the performance of the turbulence model, we plot the power spectra of a specific location in the simulation domain. The point coordinates are [x, y, z] = [4.5, 0.2, 0.3 km], where the turbulence produced by waves, rotors, or hydraulic jumpss are vigorous. In this analysis, a time series of winds is transformed into a spectral representation using the Fourier transform: where R a (t) is the auto-correlation function of the wind component a, f the cyclic frequency, t the time. Figure 7 shows an un-smoothed power spectra for the u, v, w wind components for different flow simulations as depicted in Fig. 4. The power spectra for all the wind components show an inertial subrange as has been observed in most atmospheric  Fig. 4. The xy planes are not shown because the terrain is same in y direction Fig. 7 The power spectra of u, v, and w for simulated time series at coordinate [4.5, 0.2, 0.3 km], where black color for wave only case, green for rotor case, and red for hydraulic jump case corresponding to Fig. 4. The xy planes are not shown because the terrain is same in y direction turbulent flows [59], indicating a TKE cascading process from large scale eddies to small scale eddies. The turbulence power spectra for the u component is greater than those for the v and w component which indicates greater turbulent kinetic energy (TKE). The lower frequency range length of the simulation.
In order to compare the turbulent flows at the lee-side of the ridge, we have computed averaged turbulent intensities for entire flow fields. First the TKE for each computation grid was computed using the fluctuation part of the time series. i.e.: Then the turbulent intensity (TI) is computed as follows: The TI is plotted in Fig. 8 for the x-z cross section at y = 0.2 km. Since the TI in the y direction had very little changes, the x-y cross section is not plotted. Comparing with three different flow regimes, the turbulent intensity is much higher for the wave and rotor cases. This indicated that the hydraulic jump case is the most effective in terms of kinetic energy dissipation for the stratified flow cases. It is also interesting to point out the TI above the ridges also displayed different values. This indicated that the stratification also affects the entire computation domain, not only at the lee side of the ridge.
It is should be noted that the fully three-dimensional model is important for turbulent boundary layer simulation. In our previous 2D MRTLBM simulations (not shown here), the model captured the basic structure of lee-side flow, but the turbulence was very week for the 2D simulation. The turbulent boundary layer flow is a three-dimensional process, twodimensional will restrict the turbulent vortex stretching in one y direction (north-south) so that the turbulent is reduced unrealistically.

Summary and Conclusions
A novel large-eddy simulation model is being developed at the Army Research Laboratory for simulating atmospheric flows based on the lattice Boltzmann method. To our knowledge this is the first application of a LBM model for simulating stratified atmospheric flow over a hill. The waves, rotors, and hydraulic jumps at the lee side were successfully simulated with the model. The results compared well the results of previous laboratory experiments as well as a traditional atmospheric numerical model based on finite-difference formulations, terrain following vertical coordinate systems and incompressible NS equations. Various flow regimes at the lee side of the hill including, lee-waves, rotors, and hydraulic jump, were reproduced in multiple runs of the LBM model in the physical parameter ranges used for previous laboratory and numerical experiments.
The LBM model was described in detail in the paper, including the basic formulation based on statistical principles and associated boundary conditions for simulating atmospheric stratified flow problems. The subgrid turbulence parameterization was also described. The LBM model offers simplicity in development, simplified treatment of complex boundaries, and intrinsic massive parallelism. We believe this method will be useful for many atmospheric applications, especially when very complex boundary (44) . Fig. 8 Average turbulent intensities corresponding to simulations delayed in Fig. 4. a For wave only case, b for rotor case, and c for hydraulic jump case. The xy planes are not shown because the terrain is same in y direction effects have to be resolved. We are in the process of extending LBM to three-dimensional, more realistic atmospheric flows in complex terrain, and urban environments. It should be cautioned that this model currently does not include variation of ground surface temperature, and a radiative transfer model. For a longer period of simulation, the radiative transfer process and associated land surface model is definitely required. We are in the process of developing those model components for more realistic atmospheric boundary layer simulations.
1. The transformation matrix for wind velocity between moment space and velocity space [32]: The transformation matrix for temperature [58]: