Simulation of wind gust structure in the atmospheric boundary layer with Lattice Boltzmann Method

Cold fronts occur in northern East Asia during winter and spring. After cold frontal passage, airflow is downward and accompanying strong winds fluctuate significantly; this is termed wind gusts. Analysis of observation data shows that wind gust structure has coherent characteristics. This is important for entrainment of spring dust storms into the upper boundary layer, where they are transported great distances. The Lattice Boltzmann Method (LBM) is a computational fluid technique based on the Boltzmann transport equation. The LBM has been used to study complex motion such as turbulence, because it describes motion at the micro level. In this paper, Large eddy simulation is introduced in the LBM, enabling simulation of turbulent flow in the atmospheric boundary layer. The formation and development of wind gusts are simulated, and a coherent structure with a combination of wave and vortex is obtained. This explains the mechanism of soil erosion and sand entrainment by the coherent structure of wind gusts.

Air motion in the atmospheric boundary layer is very complicated because of effects such as earth rotation, temperature stratification, water-gas transport and complex underlying surfaces [1]. Stull classified flow into three categoriesaverage wind speeds, turbulence and waves. The wave fluctuations can be formed by shear flow or average flow over obstacles, and they can be transported very far from thunderstorms and other weather systems [2]. Atmospheric fluctuations sometimes strengthen local wind shear, forming turbulence. World Meteorological Organization (WMO) [3] defines fluctuations larger than turbulence but less than average wind speeds as wind gusts. Weather patterns that produce wind gusts include strong convection, fronts, radiation inversion type low-level jets, strong low-level wind shear generated by geography, and others. Instantaneous speeds of gusts can reach 9 on the Beaufort scale, or about 80 km per hour, which can be highly destructive. For example, large equipment in a port must withstand unexpected strong gusts, and an aircraft in flight will suffer additional load generated by gusts.
In addition to this destructive power, gusts are accompanied by severe weather. The outbreak of dust storms is closely related to development of a dry squall line (a violent gust) ahead of a strong cold front [1]. This is a serious phenomenon that is common in spring in Northeast Asia. Even more common is strong wind after cold frontal passage, which occurs very suddenly with gustiness [2][3][4][5][6][7][8]. The coherent structure of wind gusts is important in the mechanism of sand-dust erosion and entrainment [9][10][11]. However, there has been little detailed study of wind gust structure, especially of the coherent characteristics of horizontal and vertical velocities.
Wind gusts are usually forecast by weather forecast models [12][13][14], time series analysis [15,16] and other methods. In engineering, the destructive effects of gusts are studied through the gust spectrum [17], gust parameterization [18][19][20][21], wind tunnel simulation and other methods [22,23]. In these methods, wind gust fluctuation is assumed to obey a Gaussian distribution. However, recent studies show that this assumption is not suitable for wind gusts after cold frontal passage [24][25][26][27][28]. During this period of strong winds, there are often regular wind gust wave packets, with periods equal to about 3-6 min. These wave packets are superimposed on the non-stationary strong basic airflow, although the unsteady basic flow also has some gustiness (very low frequency disturbances). A wind gust possesses a clearly coherent structure, i.e. vertical velocity is downward when horizontal velocity is at its peak, but is upward when horizontal velocity is at a minimum. The downward horizontal momentum is very conducive to soil erosion and sand/dust emissions. However, the strong basic flow with descending motion suppresses entrainment of dust particles, by keeping them at low levels of the atmospheric boundary layer (about 200 m height). Owing to the coherent structure of wind gusts, dust particles can effectively overcome the systematic descending air motion and penetrate mid and upper levels of the boundary layer. These particles can then propagate further and diffuse into the troposphere, where ascending air motion prevails. The coherent structure of wind gusts is an effective mechanism of soil erosion and dust entrainment. However, the related theory is obtained from analysis of observation data at single stations. It needs verification by further experimental and numerical simulation.
In the 1980s, studies were done of vertical profiles of post-cold frontal gust structure, using data from a Beijing 325 m meteorological tower [29][30][31]. In recent years, some have used advanced numerical methods, such as large eddy simulation (LES) modeling, to simulate microbursts, hurricanes, and other phenomena [32,33]. The Lattice Boltzmann Method (LBM) describes motion at the micro level, and can directly simulate complex motion such as turbulence [34][35][36][37][38][39][40]. However, for the atmospheric boundary layer, the grids cannot be fine enough, and the LBM cannot numerically simulate airflow directly. Therefore, we introduce a LES model to simulate the fine structure of wind gusts.
1 Numerical method

Lattice Boltzmann equation
The Boltzmann transport equation can be discretized with a lattice gas method, attaining the lattice Boltzmann equation, where subscript i denotes the direction of particle motion. Because of enormous computation of the collision term  i , Bhatnager et al. [41] proposed in 1954 that detailed description of the interaction between vortices and the true, unknown collision effect is unimportant. They gave a simplified collision mode, so that the Boltzmann equation was simplified to In this form of the equation, the collision term, called the relaxation term, includes a characteristic time scale  and local equilibrium distribution f eq in time. This is called the BGK model. It describes the physical nature of molecular interaction. τ is molecular collision time, which is also called relaxation time. The BGK model seems only suitable for the local equilibrium state; however, it was realized that the approximation can expand the limit of the theory as far as the duration of the relaxation time covering related physical characteristics. In the turbulent BGK model,  is substituted for by a typical turbulent relaxation time  turb . Higuera et al. [42,43] revised the lattice Boltzmann equation according to the BGK model as LBGK: The macro-density of the fluid  (x,t) and macro-velocity u (x,t) can be obtained from the particle distribution function The pressure can be obtained by p= c s 2 ; c s is velocity of sound. From a formal viewpoint, the evolution of the LBGK model is a relaxation process. Through the relaxation acceleration of the micro-particle density distribution function f i to its equilibrium state f i eq , the system quickly reaches the objective state.

DdQq model
In the LBGK model, the DdQq series of Qian et al. [44] are most widely used, wherein d denotes the space dimension and q the numbers of discrete velocities. The local equilibrium distribution function of such a model can be written as where  is fluid density, u is fluid velocity, c s is the velocity of sound, c i =ce i is discrete velocity, and  i are weights in different discrete velocity directions. The viscosity coefficient of the DdQq model is where v is the viscosity coefficient. In two-dimensional space, the DdQq model includes D2Q7 and D2Q9. In threedimensional space, it includes D3Q15, D3Q19 and D3Q27. As mentioned in the introduction, analysis of observation data shows that coherence of the gust in the vertical direc-tion (viz. the coherent structure between horizontal and vertical velocity) is primary, and the wind gust is the propagation process of two-dimensional wave and vortex. Therefore, we simulate here the two-dimensional boundary layer wind field, using the D2Q9 model. This model has a two-dimensional, square grid structure, as shown in Figure 1. Its particle velocity has nine directions, given by the following equation: (0, 0), 0; The equilibrium distribution function of the model is wherein  i =4/9 when i=0,  i =1/9 when i=1,2,3,4, and  i = 1/36 when i=4, 6,7,8. Macro-density  and macro-velocity u can be obtained: In this model, the velocity of sound is c s 2 =c 2 /3, pressure is p= c 2 /3, and viscosity is v=(    0.5)c 2 t/3.

Large eddy model
For high Re number turbulent flow, the LBM is less stable. An effective solution is to combine it with a turbulent flow model. Here, the large-eddy model is introduced into the

Test example
To test the accuracy of the calculation, we first simulate Poiseuille flow, with constant pressure gradient-driven, twodimensional cavity flow. Then we compare the simulation results with analytical solutions or data from the literature.
(i) Poiseuille flow. Between two parallel plates, there is flow with viscosity v, and a given constant pressure gradient between inlet and outlet. The analytical solution of flow velocity is where L is the distance between plates, and G p x    is the pressure gradient.
The region of simulation is 0   x   2, 0   y   1. The grids are 128×64. The initial velocity is zero. The initial density is 1.0. Re=LU max /v=1500, and U max =L 2 G/8v is the flow velocity at the center of the plates. The pressure gradient is p=0.1. After several steps, the flow field becomes stable (Figure 2).
Because of the limit of Poiseuille flow, the Reynolds number is small. Although the simulation result with the D2Q9 model is good, the advantage of LES is not demonstrated. (ii) Two-dimensional cavity flow. Next, we simulated two-dimensional cavity flow with a larger Reynolds number. The flow field is a two-dimensional square cavity. The upper boundary of the cavity moves at a constant horizontal velocity, while the other three boundaries remain stationary. The cavity is characterized by the appearance of a first-class large vortex in the central cavity, and of two secondary vortices at the two bottom angles. The value of the stream function and the central positions of these vortices are functions of Reynolds number [50] Re=LU/ν, where L is the cavity height, U is the drag speed, and ν is the viscosity coefficient. In 1982, Ghia et al. [51] discussed in detail two-dimensional cavity flow. We simulated cavity flow for Re = 5000, and compared the result with that of Ghia et al.
For smaller Reynolds number, only three vortices appear, the primary vortex in the cavity center and a pair of secondary vortices in the lower left and lower right corners. As the Reynolds number increases, a third secondary vortex appears at the upper-left corner, and the center of the primary vortex moves to the cavity center. This is consistent with previous studies (Figure 3). Comparing the horizontal velocity on the centerline, the calculation agrees with the results of Ghia et al., as shown in Figure 4. The accuracy is satisfactory.

Simulation of wind gusts in atmospheric boundary layer
(i) Boundary condition. The wind profile in the boundary layer after cold frontal passage is obtained by analyzing data from the Beijing 325 m meteorological tower. The data includes fifteen levels of low-frequency wind speed (0.05 Hz), and three levels of high-frequency, ultrasonic anemometer wind speed (10 Hz) [52]. Figure 5(a) and (b) shows vertical profiles of horizontal wind speed ( u ) and vertical velocity ( w ). They were obtained from an average of 133 hours of  There is sinking air after cold frontal passage, as shown by the figures. We used the ensemble profiles as boundary conditions. The computational domain is the boundary layer with horizontal distance x = 2 km, and vertical height z is 1 km.
Appropriate boundary conditions are necessary when solving the flow problem. Usually these boundary conditions are macro-physical quantities. However, in LBM, the evolution is the distribution function. One must determine the distribution function on the boundaries according to macro boundary conditions. According to the macro profile, the inlet distribution function is , the outlet is determined by the finite-difference velocity gradient method [54], is the tensor, I is the unit tensor, S u   is the strain tensor, and the symbol : denotes contraction between two tensors.
The lower boundary is the wall and the upper boundary is the slip condition.
(ii) Scale transformation. As with other CFD methods, there are similarities between the simulation model and actual flow field during LBM simulation of practical problems. The initial and boundary conditions of the actual flow field are transformed into model initial and boundary conditions, model calculations are done, finally reverting to the actual flow, based on the similarity relationship. Since this calculation is for the turbulent boundary layer, the Reynolds num-ber should be kept constant during the conversion process.
Assuming the actual flow field parameters are characteristic length L, characteristic velocity U, time step T, and kinematic viscosity v p during the simulation, we take the characteristic length N (if N grids, each of length l=1), characteristic velocity u, and l L N Other parameters can be similarly determined, attaining , .
Through v, the relaxation time can be obtained as

Calculation result
Based on the above calculations, we obtain the formation of wind gusts. In windy conditions after cold frontal passage, sinking air produces a strong shear. This makes the vortex oscillate up and down, causing gusty wind fluctuations. The calculations with and without the LES model were compared. We found more small-scale eddies were generated in the LES mode, and the results were closer to the actual flow field. Figure 6 is the vorticity field at the 60-min point. The upper figure was produced with the LES model, and the lower figure without it. It is seen that there are vortices at 300 m height and these oscillate up and down, because the sinking air is about 300 m thick. Since the maximum observation height is 280 m, the observed downdraft thickness may be slightly less than actual. We determined the thickness of the simulated sinking air at 300 m through radar observations. In the figure, x is horizontal distance, z is height, with units in m. The legend shows vorticity values by a gray scale; light colors show a counterclockwise vortex, dark a clockwise vortex. Figure 7 shows wind speed vectors from actual observation (10 Hz), and rising and sinking airflow.
We also calculated the vertical profile of wind speed at horizontal distances (x) of 100, 500 and 1000 m. Since 100 m is too close to the entrance, the flow does not fully develop and is influenced by inlet conditions, so the result there is not given. Figure 8 shows the vertical profile of wind speed at 500 and 1000 m. The figure shows that horizontal wind speed decreases at the air separation height (z = 280 m), which is consistent with the observed results. In Figure 5(a), the horizontal wind speed measured by wind  cups (solid dots) at 240 m height is reduced. From Figure 6, we see that in windy conditions after cold frontal passage, gustiness is strongest at the level of air separation, leading to a decrease in horizontal wind speed. This results from the sinking air producing a strong shear, causing the upward and downward motion of the vortex, which in turn generates the gustiness. As to the vertical wind speed, there is a downdraft below 320 m, with speed first increasing then decreasing with height. Above 320 m, the airflow is upward. This finding is consistent with the observations ( Figure  5(b)).
As mentioned in the introduction, we found from the data that there were gusts of main period 3-6 min in post-frontal windy conditions, and that the flow could be decomposed into basic flow, wind gust and turbulence [52]. Using this method, we decomposed the wind speed time series at x = 1000 m, z = 280 m, as shown in Figure 9. This is similar to Figure 3 of [52]. We use the extracted time series of gusts and turbulence to calculate the gust friction velocity u g* and turbulence friction velocity u t* at various heights [52]. Results are shown in Figure 10. Both are the same order of magnitude, indicating that the energy contained in wind gusts is comparable to that in turbulence. In windy conditions, therefore, wind gusts should not be ignored. Gust friction velocity is large at the height of air separation. Turbulence friction velocity is large close to the wall and decreases with height; it increases again at the height of air separation. This is consistent with observations [52]. In addition, as the horizontal distance increases, the gust friction velocity gradually increases, while the turbulent friction velocity wanes.

Conclusion
We use the LBM with the large eddy model to simulate the generation and development of wind gusts. The simulation results show that in the windy conditions after cold frontal   passage, there is airflow separation upward and downward. This creates strong shear, forms upward and downward moving wavy vortices, and causes gusty wind fluctuations. The large-eddy simulation model produces results closer to experimental observations. Further analysis shows that at the level of air separation, wind speed decreases, the speed of sinking maximizes slightly below that level, and decreases with distance from that level. By extracting the gust structure, we see that the magnitude of gust friction velocity is comparable to turbulence friction velocity. Turbulence friction velocity is large near the ground, resulting in sand erosion. As height increases, the gust friction velocity gradually increases and equals turbulence friction velocity, so that both contribute to dust entrainment. Also with increasing height, gust friction velocity becomes larger than turbulent friction velocity, and thereby becomes the major factor in dust entrainment. This is more obvious in the downstream direction.