Numerical Simulation of Coherent Structures in the Turbulent Boundary Layer under Different Stability Conditions

Coherent structures in the turbulent boundary layer were investigated under different stability conditions. Qualitative analyses of the flow field, spatial correlation coefficient field and pre-multiplied wind velocity spectrum showed that the dominant turbulent eddy structure changed from small-scale motions to large- and very-large-scale motions and then to thermal plumes as the stability changed from strong stable to neutral and then to strong unstable. A quantitative analysis of the size characteristics of the three-dimensional turbulent eddy structure based on the spatial correlation coefficient field showed that under near-neutral stability, the streamwise, wall-normal and spanwise extents remained constant at approximately 0.3 δ , 0.1 δ and 0.2 δ ( δ , boundary layer height), respectively, while for other conditions, the extent in each direction varied in a log-linear manner with stability; only the spanwise extent under stable conditions was also independent of stability. The peak wavenumber of the pre-multiplied wind velocity spectrum moves towards small values from stable conditions to neutral condition and then to unstable conditions; thus, for the wind velocity spectrum, another form is needed that takes account the effects of the stability condition.


Introduction
Over flat and homogeneous terrain, turbulence is a very important factor affecting turbulent motions near the ground. Consequently, the structure of the turbulent boundary layer depends mainly on vertical turbulence transport which has a great influence on mass, momentum, heat, water vapor and energy exchanges [1] with turbulent coherent structures playing an important role [2,3]. Field measurements, laboratory experiments and numerical simulations of the turbulent boundary layer have revealed the existence of a wide variety of coherent structures such as low-speed streaks, hairpin vortices, large-scale motions and very large-scale motions in canonical wall-bounded flows including pipe flows, channel flows and turbulent boundary layer flows [2,[4][5][6][7][8][9][10]. Actually, coherent structures are characterized by a specific scale and lifetime, the recurrence of characteristic and flow-specific organization and regions of concentrated vorticity. Specifically, Theodorsen [4] proposed the existence of hairpin vortices, and such vortices, together with streamwise packets, have been suggested to be the primary elements forming coherent structures in the neutral turbulent boundary logarithmic layer under experimental conditions [6]. However, most studies have focused on neutral conditions. Under neutral conditions, turbulent structures observed in the field are similar to those found in experiments. Furthermore, Marusic et al. [3] indicated that turbulence and even the pressure field in the turbulent boundary layer are also similar to those in canonical wall-bounded flows, and they suggested that the reason for this is the presence of coherent structures.
The prevalent view regarding the formation of coherent structures in the turbulent surface layer is the "bottom-up" mechanism, in which active structures emerge mainly at the surface and then grow up into the outer layer [8,[11][12][13]. As is well known, the turbulent boundary layer is usually classified as being stable, neutral or unstable on the basis of the dominant mechanism of turbulence generation and stability [14,15]. Differences in stability will affect energy transport and fluxes between the near-surface layer and upper free layer, through the formation and development of coherent structures [16][17][18][19][20]. The stability depends on its ability to resist three-dimensional turbulent motions. In stable conditions, even slight three-dimensional turbulent motions are difficult or even impossible. By contrast, in unstable conditions, slight three-dimensional turbulent motions tend to become stronger, leading to turbulence and convective activity. Li and Bou-Zeid [18] concluded that coherent structures in the atmosphere must be different from those found in experiments, because in the turbulent boundary layer they can be affected by the surface buoyancy flux under stable (nocturnal) or unstable (daytime) conditions which are not taken into account in most experiments. Specifically, unstable turbulent boundary layers include three kinds of coherent structures: under slightly unstable conditions, large-scale horizontal roll vortices occur, owing to a combination of shear and buoyancy [21][22][23][24]; under moderately unstable conditions, individual hairpins arise from the ground, leading to an increase in the inclination angle of the vortex packets [25][26][27][28]; and under strongly unstable conditions, plume-shaped structures appear, together with streamwise elongated structures, and the magnitude of the vertical turbulent momentum flux increases [17,25,29,30]. Stable turbulent boundary layers can generally be divided into two kinds: weakly stable and strongly stable [31]. Under weakly stable conditions, because turbulence is nearly continuous and quasi-steady, Monin-Obukhov similarity theory applies [32]. Under strongly stable conditions, small eddy motions dominate in the boundary layer, and the layer is much shallower than those under neutral or unstable conditions; however, our understanding of this is still quite limited [31,33,34] partly because of the intrinsic complexity of fluid dynamical phenomena such as intermittency and Kelvin-Helmholtz instability [31,33,35].
In the present study, to obtain a better understanding of turbulent boundary layer flow under different stability conditions, a large eddy simulation (LES) approach was adopted. The LES technique was proposed for meteorological applications by Deardorff [36,37]; since then, it has become an indispensable tool to study the turbulent boundary layer [38][39][40][41]. In LES, large-scale structures are resolved explicitly, while the effect of small-scale structures to large-scale structures is not resolved but is parameterized using a subgrid-scale (SGS) model. Thus, because the unstable turbulent boundary layer is dominated by large-scale structures, LES is an appropriate approach [42][43][44]. Large eddy simulation has also been used successfully for neutral conditions [17,20,42,44,45]. However, simulation of the stable turbulent boundary layer faces a number of challenges. These arise because the way in which the radiative cooling of the ground leads to stable stratification modifying the flow turbulence is still not completely understood [46], and because the characteristic size of eddies decreases with increasing stability. The application of LES to stable conditions has, therefore, been mostly confined to weakly or moderately stable cases [47,48], and this approach has not been considered appropriate for strongly stable conditions [34,49]. However, Churchfield et al. [50] have suggested an approach that is capable of simulating complex terrain and stable conditions using wall model LES (WMLES), and it is this approach that is adopted in the present study.
Turbulent boundary layers have a complex structure and involve some phenomena that are still not completely understood, especially at high Reynolds numbers [9]. Therefore, the objective of this study was to investigate the characteristics of the turbulent boundary layer under different stability conditions and to quantify the variations in the shapes and sizes of eddy structures with stability. The LES model and simulation setup are described in Section 2, and the model is validated in Section 3. Simulation results are discussed in Section 4. Finally, conclusions are given in Section 5.

Simulation Details
The governing equations of the LES model are: For potential temperature transport, the governing equation is: where x i represent Cartesian coordinates, u i represents the corresponding velocity components, p represents the pressure which, divided by the fluid density and τ D ij , represents the SGS stress, and f T i represents the density-normalized forces. The Coriolis force was not considered, the same as in Fang and Porté-Agel [20]. A bar on a variable indicates the mean grid-filtered quantity. θ is the resolved potential temperature, and q j represents the flux of temperature due to the viscous and sub-filter-scale effects. The details of how the SGS momentum flux τ D ij was modelled is shown in Ren et al. [51]. The heat flux q j at the surface was modelled using a similar method [38]: where, in the simulation, the total average temperature flux is specified, and the fluctuating temperature flux q tot 3 is created by the wall model. For numerical discretization, the finite volume method was used; for the viscous and convective terms, a second-order central difference scheme was adopted; and for the unsteady term, a second-order backward scheme was adopted.
Schumann's model [52] was used as the wall model to model the near wall shear stress, and Moeng's model [38] was adopted for the total temperature flux at the surface. The subgrid stress (SGS) model is the one-equation eddy viscosity model [53]; this is because the one-equation eddy viscosity model does not require the assumption in the algebraic eddy viscosity models, that is, local balance between the production and dissipation of SGS energy. The SGS stress can be modelled as: where τ D ij,sgs is the SGS stress part of the total stress, τ D ij . S ij is the resolved-scale strain rate tensor, δ ij is the Kronecker delta function, k sgs is the SGS kinetic energy; and υ sgs is the SGS stress eddy viscosity.
Accounting for the historic effect of k sgs , because of diffusion, production and dissipation, a transportation equation was derived: where ∆ represents the filter scale, and C k and C ε are the model constants. In the present study, the size of the computational domain was 5 × 1 × 0.3 km for all cases. When the hexahedral mesh was applied, the mesh grid size was 5 m in each direction from the ground surface to 0.2 km, and above 0.2 km, the mesh grid size was 10 m. For the streamwise and spanwise directions, periodic boundary conditions were applied; this agrees well with Fang and Porté-Agel [20] of LES of turbulent boundary layer flow and with Lee et al. [54] of DNS of channel flow; for the upper boundary, the slip boundary condition was adopted; and for the solid wall, the wall model was used. For the initial conditions, the logarithmic mean wind speed profile was used in all cases (gradient wind height of 250 m and a gradient wind speed of 8 m/s).
The computation timestep was 0.5 s which ensures the stability and viscous stability by a given Courant number. First, to obtain a fully developed flow field and quasi-steady state, a spin-up simulation was performed. For the statistical analysis, to ensure statistical convergence, over 40,000 timesteps were simulated. For the perturbation component, the initial peak was set at the height of 0.015% of the computational domain, and for the horizontal direction, the maximum perturbation was 0.25 m/s; at the beginning of the simulation, for the streamwise and spanwise direction, the number of periods of the perturbations were 4.0 and 20.0, respectively.
The aim of this numerical simulation was to calculate the forms of coherent structures under specifically assumed stability conditions and to quantify the variations in the shapes and sizes of eddy structures with stability. The simulation was conducted within a pre-defined stability range. As is well known, stability refers to the air's tendency to either resist vertical movement (stability) or to rise and create storms (instability) which can be expressed in terms of a dimensionless stability parameter, −z/L, where z represents the height above the surface ground, and L represents the Obukhov length, given by L = −u 3 * T v /(kgQ v0 ), where k represents the von Kármán constant, u * represents the friction velocity, Q v0 is the kinematic virtual temperature flux at the surface, T v is a reference virtual temperature, and g is the gravitational acceleration. L expresses the relative roles of shear and buoyancy in the production/consumption of turbulence kinetic energy with positive and negative values indicating stable and unstable conditions, respectively, and it approaches zero in the limit of neutral stratification.  turbulent boundary layer is under stable conditions, the turbulent motions are weakly connected to the ground, and the height is no longer a similarity parameter of this motion, except in the region near the ground. As can be seen from Figure 2c, the normalized standard deviation of the vertical velocity component behaves as expected for stable conditions. Also, similar to the ratio of horizontal and vertical heat fluxes, the peak of * / w u σ shifts to smaller / z L − as stability increases, except for neutral conditions; this is related to the upper location of thermal plumes. . In (a), for unstable conditions, the blue circles, red crosses, green plus signs, cyan stars, magenta squares, yellow asterisks, black diamonds and blue rightward-pointing triangles represent cases 17, 16, 15, 14, 13, 12, 11 and 10, respectively, from very to slightly unstable. The curve corresponds to Equation (7). In (b), for stable conditions, the blue circles, red crosses, green plus signs, cyan stars, magenta squares, yellow asterisks, black diamonds and blue rightward-pointing triangles represents cases 8, 7, 6, 5, 4, 3, 2, and 1, respectively, going from slightly to very stable. Ratio of horizontal and vertical heat fluxes u θ /w θ versus −z/L. In (a), for unstable conditions, the blue circles, red crosses, green plus signs, cyan stars, magenta squares, yellow asterisks, black diamonds and blue rightward-pointing triangles represent cases 17, 16, 15, 14, 13, 12, 11 and 10, respectively, from very to slightly unstable. The curve corresponds to Equation (7). In (b), for stable conditions, the blue circles, red crosses, green plus signs, cyan stars, magenta squares, yellow asterisks, black diamonds and blue rightward-pointing triangles represents cases 8, 7, 6, 5, 4, 3, 2, and 1, respectively, going from slightly to very stable. neutral stability conditions; and (c) stable conditions. The symbols are the same as in Figure 1, except that the curve in (a) now corresponds to Equation (8).

Results
The above verification analysis has illustrated the feasibility of the WMLES approach and the accuracy of simulation of turbulent boundary layer flow under different stability conditions. In the following, we investigate the characteristics of turbulent boundary layers under different stability conditions and quantify the variations in the shape and size of eddy structures with stability.

Validation
The basic statistical properties of the velocity under neutral stability conditions were compared with those from previous studies, and good agreement was found when dimensioned both by inner and outer variables; the details are shown in previous studies by Ren et al. [51] and Ren et al. [55].
The ratio of the heat fluxes in the horizontal and vertical directions under unstable conditions is shown in Figure 2. Wyngaard and Cote [56] proposed the following formula for this ratio: where ϕ m = (kz/u * )(∂u/∂z) represents the non-dimensional wind shear function, u * is the surface friction velocity, ϕ h = (kz/u * )(∂θ/∂z) is the non-dimensional potential temperature gradient function, and a is a constant. In the present study, [57,58], and a = 5 [57][58][59]. In Figure 1, the ratio of the horizontal and vertical heat fluxes decreases with increasing −z/L away from the near-surface region, and u θ /w θ tends towards 0 which indicates that at higher elevations the shear stress is very small compared with the thermal plume force. The simulation results show good agreement with theory under strongly unstable conditions, but under moderately or slightly unstable conditions, they deviate from the theoretical values. This deviation occurs because under these conditions, the shear stress driving the horizontal motions near the surface is greater than under strongly unstable conditions, while the vertical transport is weaker than under strongly unstable conditions. Therefore, Equation (7) is appropriate for strongly unstable conditions only, and perhaps, after the forms of a, ϕ m and ϕ h are changed in the present study, it will also prove to be appropriate for moderately or slightly unstable conditions. Another phenomenon revealed in Figure 1 is a shift in the peak of w θ /w θ towards larger −z/L; this is because as conditions become more unstable, thermal plumes reach higher. Similarly, the peak of u θ /w θ shifts to smaller −z/L as the stability increases under stable conditions. This suggests that as conditions become more unstable, −z/L becomes larger which is related to the upper location of thermal plumes. Figure 2 shows the normalized standard deviation of the vertical velocity component, σ w /u * , which is important for describing stability. The curves in Figure 2 represent the formula proposed by Panofsky et al. [59] in the case of the unstable condition: where u * is the friction velocity, L represents the Obukhov length, and z is the height above the ground. It can be seen that σ w /u * remains constant under slightly or moderately unstable conditions but increases with increasing height under strongly unstable conditions. The simulation results show good agreement with Equation (8) which means that σ w /u * obeys Monin-Obukhov similarity theory under unstable conditions, fitting well with the 1/3 power law. Under neutral conditions, σ w /u * is constant, independent of height and equal to 1.32 (shown in Figure 2b); indeed, a large number of observations show that under neutral conditions, σ w /u * should be constant [58,[60][61][62][63][64]. In particular, Bowne and Ball [60] and Counihan [61] both concluded that this constant value is 1.3. Roth [64] summarized the values of different researchers, which ranged from 1.1 to 1.7. When the turbulent boundary layer is under stable conditions, the turbulent motions are weakly connected to the ground, and the height is no longer a similarity parameter of this motion, except in the region near the ground. As can be seen from Figure 2c, the normalized standard deviation of the vertical velocity component behaves as expected for stable conditions. Also, similar to the ratio of horizontal and vertical heat fluxes, the peak of σ w /u * shifts to smaller −z/L as stability increases, except for neutral conditions; this is related to the upper location of thermal plumes.

Results
The above verification analysis has illustrated the feasibility of the WMLES approach and the accuracy of simulation of turbulent boundary layer flow under different stability conditions. In the following, we investigate the characteristics of turbulent boundary layers under different stability conditions and quantify the variations in the shape and size of eddy structures with stability.

Flow Field
We now turn to an analysis of the flow field with the aim of providing a more intuitive view of its basic characteristics and those of the turbulent eddy scale. We consider a horizontal plane (whole domain) at the height of 0.2δ. Figure 3 shows the flow field for all 17 cases. On the whole, under stable conditions, the dominant turbulent eddy structure was small which is consistent with many previous results [31,33,34]. This is because, under stable conditions, the ground temperature is lower than that of the near-surface flow, and the vertical development of the turbulent eddy structure is limited which suppresses the development of three-dimensional turbulent eddy structure thus small-scale motion is prevalent. The dominant turbulent eddy structure scale increases as the stability decreases, and for the weakly stable condition in case 8, very-large-scale motion structures [8,10,65] even occur, just as under neutral conditions. It is obvious under neutral conditions that both large-scale and very-large-scale motion structures are visualized, based on the streamwise-elongated volumes associated with negative streamwise velocity, which is consistent with the results of Dennis and Nickels [66,67], Lee and Sung [65], Baltzer et al. [68] and Lee et al. [54]. Under slightly and moderately unstable conditions, obvious large-scale horizontal roll vortices form which are usually attributed to the combined effects of shear and buoyancy [21,22,24]. These can be seen as corresponding to the extended turbulent eddy structure that appears under neutral stability conditions. Under strongly unstable conditions, the effects of buoyancy dominate those of shear, so there is an obvious distinction between low-velocity and high-velocity regions which extends over the whole spanwise direction. This is consistent with the experimental results of Hommema and Adrian [25].
In Figure 3, going from the strongly stable condition to the neutral stable condition to the strongly unstable condition, some interesting phenomena can be seen. First, the form and shape of the dominant turbulent eddy structure changed from small-scale motions to large-scale and very-large-scale motions and then to regional motions as the turbulent boundary layer became unstable. Second, the dominant size of the turbulent eddy structure varied as the stability changed; a quantitative analysis of this is presented in Section 4.2 based on the three-dimensional spatial correlation coefficient field.

Spatial Correlation Coefficient Field
To quantify the size characteristics of the three-dimensional turbulent eddy structure under different stability conditions, a flow field correlation analysis was performed with reference location (x, y, z) = (2500 m, 500 m, 35 m). Figure 4 presents contours of the spatial correlation coefficient for all 17 cases. The plots in the left column are for the xz plane, and those in the right column are for the yz plane; thus, it is possible to visualize the three-dimensional turbulent eddy structure. The conclusion is similar to that from the flow field analysis, namely, that with increasing stability, the dominant size of the turbulent eddy decreases. This clearly shows that under weakly stable, neutral and slightly and moderately unstable conditions, very-large-scale motion structures appeared, while strongly unstable conditions led to thermal plumes. Meanwhile, an interesting phenomenon was that the spanwise size of the dominant vortices remained almost constant, except under strongly unstable conditions which may be due to the fact that the effect of the shear force was mainly in the streamwise direction and that of the buoyancy force was mainly in the wall-normal direction; however, under strongly unstable conditions, the temperature of the ground was higher than that of the adjacent near-surface air, and so the buoyancy force had a positive effect on the upward motion of air near the ground in the whole spanwise direction. The spanwise and wall-normal velocity correlations were different from the streamwise velocity correlation. For the spanwise velocity correlation under stable conditions, the spanwise size of the correlation structure was smaller than the streamwise velocity correlation; moreover, the degree of tilt was smaller than for the streamwise velocity correlation. Under unstable conditions, even light unstable conditions, the plume-shaped structure occurred, meaning that this structure form occurred earlier than the streamwise velocity correlation. For the wall-normal velocity correlation, from stable to unstable conditions, the spanwise and wall-normal sizes of the correlation structure became wide and deep which is very different to the other two direction velocity correlations. These results are not shown here. and slightly and moderately unstable conditions, very-large-scale motion structures appeared, while strongly unstable conditions led to thermal plumes. Meanwhile, an interesting phenomenon was that the spanwise size of the dominant vortices remained almost constant, except under strongly unstable conditions which may be due to the fact that the effect of the shear force was mainly in the streamwise direction and that of the buoyancy force was mainly in the wall-normal direction; however, under strongly unstable conditions, the temperature of the ground was higher than that of the adjacent nearsurface air, and so the buoyancy force had a positive effect on the upward motion of air near the ground in the whole spanwise direction. The spanwise and wall-normal velocity correlations were different from the streamwise velocity correlation. For the spanwise velocity correlation under stable conditions, the spanwise size of the correlation structure was smaller than the streamwise velocity correlation; moreover, the degree of tilt was smaller than for the streamwise velocity correlation. Under unstable conditions, even light unstable conditions, the plume-shaped structure occurred, meaning that this structure form occurred earlier than the streamwise velocity correlation. For the wall-normal velocity correlation, from stable to unstable conditions, the spanwise and wall-normal sizes of the correlation structure became wide and deep which is very different to the other two direction velocity correlations. These results are not shown here. (e) (f)             [70] as twice the distances from the self-correlation peak to the most distant locations on the 0.5 contour of the two-point correlation of the fluctuating velocity in the downstream, spanwise and wall-normal directions, respectively. We can see that under strongly and moderately stable conditions, the streamwise extent Lx uu increased in a log-linear manner as −z/L increased, and that it behaved similarly under strongly unstable conditions. However, under near-neutral stability conditions, Lx uu was independent of −z/L and remained constant at approximately 0.3δ. Similar behaviour was found for the wall-normal extent Lz uu which remained constant at approximately 0.1δ. From strongly stable through neutral to slightly unstable conditions, the spanwise extent Ly uu remained constant at approximately 0.2δ, while under strongly unstable conditions, Ly uu also increased in a log-linear manner with increasing −z/L. This corresponds to the behaviour of the spatial correlation coefficient field.

Pre-Multiplied Spectral Analysis
Pre-multiplied spectral analysis is a popular approach for determining the distribution characteristics of vortices of different scales in canonical wall-bounded flows [8,10,71,72]. Figure 6 shows the multipoint averaged [73] pre-multiplied wind velocity spectra for all 17 cases. As we know, the most energetic point in the pre-multiplied wind velocity spectrum represents the dominant turbulent eddy structure. Therefore, under stable conditions, the peak wavenumber of the wind velocity spectrum is located in the high-wavenumber region which indicates that the dominant vortices are small. With decreasing stability, the peak wavenumber moved to smaller values which means that the size of the dominant turbulent eddy structure increased and very-large-scale motion structures occurred under the slightly stable conditions of case 8. Here, it should be noted that verylarge-scale motion structures are those with turbulent eddy size larger than 3δ . Under near-neutral stability conditions (slightly stable, neutral and slightly unstable conditions), the stationary stage of the pre-multiplied wind velocity spectrum has a large width, because the stationary stage corresponds to the very-large-scale motion structure [6,8,10,12]. Under strongly unstable conditions, the peak wavenumber is located in the small wavenumber region which means that the dominant turbulent eddy structure is very large, and the shape of the wind velocity spectrum is similar to that of a typhoon wind velocity spectrum near the eye-wall region. This may be because the eye-wall region of a typhoon is also characterized by strongly convective conditions, accompanied by a strong updraft. For spanwise and wall-normal pre-multiplied wind velocity spectra, the maximum energy points are consistent with the streamwise pre-multiplied wind velocity spectrum, and as the instability increases, it moves to a small wavenumber. The difference is that the steady stage is more obvious in the streamwise pre-multiplied wind velocity spectrum than the other two directions, and that the spanwise and wall-normal pre-multiplied wind velocity spectra have obvious cross-over phenomenon at different heights, which does not occur for streamwise pre-multiplied wind velocity spectrum. These results are not shown here.

Pre-Multiplied Spectral Analysis
Pre-multiplied spectral analysis is a popular approach for determining the distribution characteristics of vortices of different scales in canonical wall-bounded flows [8,10,71,72]. Figure 6 shows the multipoint averaged [73] pre-multiplied wind velocity spectra for all 17 cases. As we know, the most energetic point in the pre-multiplied wind velocity spectrum represents the dominant turbulent eddy structure. Therefore, under stable conditions, the peak wavenumber of the wind velocity spectrum is located in the high-wavenumber region which indicates that the dominant vortices are small. With decreasing stability, the peak wavenumber moved to smaller values which means that the size of the dominant turbulent eddy structure increased and very-large-scale motion structures occurred under the slightly stable conditions of case 8. Here, it should be noted that very-large-scale motion structures are those with turbulent eddy size larger than 3δ. Under near-neutral stability conditions (slightly stable, neutral and slightly unstable conditions), the stationary stage of the pre-multiplied wind velocity spectrum has a large width, because the stationary stage corresponds to the very-large-scale motion structure [6,8,10,12]. Under strongly unstable conditions, the peak wavenumber is located in the small wavenumber region which means that the dominant turbulent eddy structure is very large, and the shape of the wind velocity spectrum is similar to that of a typhoon wind velocity spectrum near the eye-wall region. This may be because the eye-wall region of a typhoon is also characterized by strongly convective conditions, accompanied by a strong updraft. For spanwise and wall-normal pre-multiplied wind velocity spectra, the maximum energy points are consistent with the streamwise pre-multiplied wind velocity spectrum, and as the instability increases, it moves to a small wavenumber. The difference is that the steady stage is more obvious in the streamwise pre-multiplied wind velocity spectrum than the other two directions, and that the spanwise and wall-normal pre-multiplied wind velocity spectra have obvious cross-over phenomenon at different heights, which does not occur for streamwise pre-multiplied wind velocity spectrum. These results are not shown here.
Energies 2020, 13, x FOR PEER REVIEW 14 of 19 differences in the wind velocity spectra and allows us to obtain a universal wind velocity spectrum that is suitable for all stability conditions.   Figure 6. Pre-multiplied wind velocity spectra of the streamwise velocity fluctuation for each case at different heights: (a) to (q) are for cases 1 to 17, respectively. The dashed and solid black lines represent 0 and -2/3 slope, respectively.

Conclusions
We investigated the coherent structures in the turbulent boundary layer under different stability conditions and performed qualitative and quantitative analyses of the shape and size of eddy structures.
First, the WMLES approach and the basic variables of the simulation results were validated by comparing these results with theoretical predictions, observed data, experimental data and other simulation results.
Then, based on a flow field analysis, it was found that the form and shape of the dominant turbulent eddy structure changed from small-scale motions to large-scale and very-large-scale motions and then to thermal plumes as the stability condition changed from strong stable conditions to neutral stability conditions and then, to strongly unstable conditions. A spatial correlation analysis was performed, and the results led to similar conclusions as those obtained from the flow field It can be seen from Figure 6 that as the stability conditions changed from stable to neutral to unstable, the peak wavenumber of the pre-multiplied wind velocity spectrum moved towards smaller values which is consistent with the above analyses of the flow field and the spatial correlation coefficient field, and the stationary stage of the wind velocity spectrum existed only under near-neutral conditions. Under strongly stable or unstable conditions, the small-scale and large-scale turbulent eddy structures were the dominant structures, respectively, and only one dominant peak existed in the pre-multiplied wind velocity spectrum.
The canonical wind velocity spectral styles are based on the neutral stability condition, but in this paper, we found that the different stability conditions gave different wind velocity spectral forms. Therefore, in future studies, we should try to find a stability parameter that represents the differences in the wind velocity spectra and allows us to obtain a universal wind velocity spectrum that is suitable for all stability conditions.

Conclusions
We investigated the coherent structures in the turbulent boundary layer under different stability conditions and performed qualitative and quantitative analyses of the shape and size of eddy structures.
First, the WMLES approach and the basic variables of the simulation results were validated by comparing these results with theoretical predictions, observed data, experimental data and other simulation results.
Then, based on a flow field analysis, it was found that the form and shape of the dominant turbulent eddy structure changed from small-scale motions to large-scale and very-large-scale motions and then to thermal plumes as the stability condition changed from strong stable conditions to neutral stability conditions and then, to strongly unstable conditions. A spatial correlation analysis was performed, and the results led to similar conclusions as those obtained from the flow field analysis. An interesting phenomenon was also found, namely, that the spanwise size of the dominant vortices remained almost constant, except under strongly unstable conditions which may be because the main effect of shear is in the streamwise direction while that of buoyancy is in the wall-normal direction. From an analysis of the pre-multiplied wind velocity spectrum, it was found that as the stability condition goes from stable to neutral to unstable, the peak wavenumber of the pre-multiplied wind velocity spectrum moves towards smaller values, and the stationary stage of wind velocity spectrum only exists under near-neutral stability, whereas under strongly stable or unstable conditions, the smalland large-scale turbulent eddy structures have an overwhelming advantage, so only one dominant peak exists in the pre-multiplied wind velocity spectrum. In future work, we need to consider a new style for the wind velocity spectrum that takes account of the effects of stability.
Finally, a quantitative analysis of the size characteristics of the three-dimensional turbulent eddy structure under different stability conditions based on the spatial correlation coefficient field showed that for near-neutral stability, the streamwise, wall-normal and spanwise extents remain constant at approximately 0.3δ, 0.1δ and 0.2δ, respectively. For other conditions, it was found that the streamwise, wall-normal and spanwise extents increase in a log-linear manner with increasing stability, except for the spanwise extent under stable conditions which is independent of stability.