Numerical Investigation of Very-Large-Scale Motions in a Turbulent Boundary Layer for Different Roughness

Wall-model large eddy simulations (WMLES) are conducted to investigate the spatial features of large-scale and very-large-scale motions (LSMs and VLSMs) in turbulent boundary flow in different surface roughnesses at a very high Reynolds number, O (106–107). The results of the simulation of nearly smooth cases display good agreement with field observations and experimental data, both dimensioned using inner and outer variables. Using pre-multiplied spectral analysis, the size of VLSMs can be reduced or even disappear with increasing roughness, which indirectly supports the concept that the bottom-up mechanism is one of the origins of VLSMs. With increases in height, the power of pre-multiplied spectra at both high and low wavenumber regions decreases, which is consistent with most observational and experimental results. Furthermore, we find that the change in the spectrum scaling law from −1 to −5/3 is a gradual process. Due to the limitations of the computational domain and coarse grid that were adopted, some VLSMs and small-scale turbulence are truncated. However, the size of LSMs is fully accounted for. From the perspective of the spatial correlation of the flow field, the structural characteristics of VLSMs under various surface roughnesses, including three-dimensional length scales and inclination angles, are obtained intuitively, and the conclusions are found to be in good agreement with the velocity spectra. Finally, the generation, development and extinction of three-dimensional VLSMs are analyzed by instantaneous flow and vorticity field, and it shows that the instantaneous flow field gives evidence of low-speed streamwise-elongated flow structures with negative streamwise velocity fluctuation component, and which are flanked on each side by similarly high-speed streamwise-elongated flow structures. Moreover, each of the low-speed streamwise-elongated flow structure lies beneath many vortices.


Wall Model LES (WMLES)
For incompressible flow, the time-dependent filtered N-S equations are written as follows: where x i are Cartesian coordinates, u i are corresponding velocity components; p is pressure divided by fluid density; τ D ij is the total stress; and f T i is the other density-normalized forces. The bar on top of the variables represents a filtering operator.
In Equation (1), on the left hand side, the first term represents the time change rate, the second term represents convection, on the right hand side, the first term gives the density normalized pressure gradient, the second term represents the stresses, and the third term is the other density-normalized forces.
SGS model used to represent the effect of unresolved motions on the resolved motions is critical for accurate simulations, especially for the specific purpose of capturing the dynamics of large-scale structures. In the present study, the one-equation of the eddy viscosity model [31] was adopted. The reason for this being that the crucial feature of the Smagorinsky model lies in the assumption of the balance between SGS energy production and dissipation. On using the ensemble-averaging procedure in place of the filtering in LES, the balance of energy production and dissipation holds only in some Energies 2020, 13, 659 4 of 26 restricted flow situations. In fact, this balance breaks down in the vicinity of walls in channel and boundary-layer flows, as in jets and wakes. Moreover, for the Smagorinsky model, the eddy-viscosity representation for the Reynolds stress based on the ensemble mean cannot reproduce the anisotropy of turbulence observed in some typical flows like channel flows. However, the one-equation eddy viscosity model is not based on the assumption of SGS energy balance and uses a more accurate representation for SGS Reynolds stress, the deficiency of local balance assumption adopted in algebraic eddy viscosity models can be overcome. At high Reynolds number flows and for coarse grid resolution such a phenomenon may occur, which fits perfectly with the simulation considered in this study.
The SGS stress can be modeled as: where τ D ij,sgs SGS stress part of the total stress τ D ij . S ij resolved-scale strain rate tensor, δ ij Kronecker delta function, k sgs SGS kinetic energy; υ sgs SGS stress eddy viscosity.
Accounting for the historical effect of k sgs because of diffusion, production, and dissipation, a transportation equation is derived: where ∆ filter scale, C k and C ε model constants. C k = 0.0673 and C ε = 0.93 were adopted in this study. The computational requirement of high Reynolds number wall-bounded flow of fully-resolved LES dependents strongly on Reynolds number. Piomelli and Balaras [29] pointed out that at high Reynolds number o (10 6 ), for smooth boundary layers flow, 99% of the computational mesh grids are used to resolve the inner layer whose thickness is only 10% of the boundary layer height; furthermore, more than 99% of the computational mesh grids is required for the rough boundary layers flow. Consequently, the only economical way to perform LES of high Reynolds number turbulent boundary layer flow is just to compute the outer layer. The outer-flow eddies can determine the grid size in this case, while for the near wall eddies, the grid is too coarse, so must model the wall layer. In particular, near the wall, discrete differentiation cannot evaluate the momentum flux, since the coarse grid cannot resolve the sharp velocity gradient and the quasi-streamwise in this region. Moreover, Piomelli and Balaras concluded that the WMLES tend to be more accurate at high Reynolds number than at moderate or low Reynolds number.
Moreover, for the fully-resolved LES, the SGS model and numerical errors very large in the near-wall region usually; thus, it is conceivable that comparison with fully-resolved LES, using "perfect" wall model can get a more accurate predictions.
Another important point we need to consider involves the case of an atmospheric surface covered with high roughness elements (e.g., sands, rocks, vegetation) that would greatly increase the computation cost when used in conjunction with the fully-resolved LES.
Thus, in the present study, the surface stress model is the following matrix: and, the Schumann's model [30] was adopted: 1/2 (7) where 1/2 is at first level cell centers near the surface; the angle brackets represent horizontal average at a height. The friction velocity: u 2 * = ( τ tot 13 2 + τ tot 23 2 1/2 (8) which rough wall log law: where f (L) atmospheric stability-related function, L Obuhkov length, z 0 aerodynamics roughness height.
In this study, f (L) is zero for neutral stability, z 0 is related to roughness elements shape, height, and distribution.

Simulation Setup
1/2 * 13 23 =( + ) tot tot u   (8) which rough wall log law: where ( ) f L atmospheric stability-related function, L Obuhkov length, 0 z aerodynamics roughness height. In this study, ( ) f L is zero for neutral stability, 0 z is related to roughness elements shape, height, and distribution. In this paper, the five different cases of surface roughness are investigated: 0 z = 0.0002 m (case 1, Re   3.538  10 6 ), 0 z = 0.002 m (case 2, Re   4.894  10 6 ), 0 z = 0.02 m (case 3, Re   6.559  10 6 ), 0 z = 0.2 m (case 4, Re   9.927  10 6 ), and 0 z = 2.0 m (case 5, Re   1.199  10 7 ). Figure 1 gives the schematic diagram of the computational domain (5000  1000  300 m for case 1, 5000  1000  350 m for case 2, 5000  1000  400 m for case 3, 5000  1000  500 m for case 4, 5000  1000  600 m for case 5) adopted in this study. The streamwise, spanwise, and wall-normal coordinates are x , y , and z . At the solid wall, wall model was adopted, and at the upper surface, slip boundary conditions were adopted. For the streamwise and spanwise directions, periodic boundary conditions were applied, same with the LES TBL flow simulation by Fang and Porté-Agel [27], and DNS channel flow simulation conducted by Lee et al. [6]. The logarithmic mean wind-speed profile was used for the inflow, specifically for five cases:  The hexahedral meshes type was adopted in the simulation. The computational domain at the height below 200 m, the mesh grid size is 5 m in each direction. Except for this region, the mesh grid size is 10 m in each direction.

Simulation Setup
The time step is set to 0.5 s in order to satisfy the Courant number for the stability condition and viscous stability. To guarantee a fully developed flow field and quasi-steady conditions, a spin-up simulation is first made for enough time. To ensure statistical convergence, statistical calculations are then performed over 40,000 timesteps. At 0.015% of domain's height, the initial peak in perturbation magnitude occurs, with both the maximum perturbations of streamwise/spanwise flow near the surface running at 0.25 m/s, a 4.0/20.0 ratio was adopted at the beginning of the simulation for total periods of streamwise/spanwise perturbations.
The simulations were performed using the Simulator for Wind Farm Application (SOWFA) [32], which is developed using OpenFOAM C++ library. The numerical discretization was performed by a finite volume method. A second-order central difference scheme was used for the convection and The hexahedral meshes type was adopted in the simulation. The computational domain at the height below 200 m, the mesh grid size is 5 m in each direction. Except for this region, the mesh grid size is 10 m in each direction.
The time step is set to 0.5 s in order to satisfy the Courant number for the stability condition and viscous stability. To guarantee a fully developed flow field and quasi-steady conditions, a spin-up simulation is first made for enough time. To ensure statistical convergence, statistical calculations are then performed over 40,000 timesteps. At 0.015% of domain's height, the initial peak in perturbation magnitude occurs, with both the maximum perturbations of streamwise/spanwise flow near the surface Energies 2020, 13, 659 6 of 26 running at 0.25 m/s, a 4.0/20.0 ratio was adopted at the beginning of the simulation for total periods of streamwise/spanwise perturbations.
The simulations were performed using the Simulator for Wind Farm Application (SOWFA) [32], which is developed using OpenFOAM C++ library. The numerical discretization was performed by a finite volume method. A second-order central difference scheme was used for the convection and viscosity terms, and a second-order backward scheme was employed for the unsteady term.

Validations
Because the wall layer is modeled in WMLES, the accuracy of the simulation results must be validated before the study can proceed. Here, the simulation results based on non-dimensional inner variables z + = zu τ /v and u τ , and outer variables δ and U, are compared with the results of DNS, experiments, and field measurements.

Validation by Non-Dimensional Inner Variables
To confirm that the WMLES approach is suitable for TBL flows at high Reynolds numbers, variations in mean velocity, Reynolds stress, and turbulent variation with height and Reynolds number are compared with theoretical results and experimental, field-measured data for the TBL flow.
The non-dimensional mean streamwise velocity profiles, shown in Figure 2. The simulated non-dimensional mean velocities U + = U/u τ are displayed at different heights from case 1. The mean velocity profile agrees well with the logarithmic law and Qingtu Lake Observation Array (QLOA) field-measured data, indicating that the WMLES is highly accurate for TBL flow, despite the down shift in ∆U + due to ground roughness.

Validations
Because the wall layer is modeled in WMLES, the accuracy of the simulation results must be validated before the study can proceed. Here, the simulation results based on non-dimensional inner variables / z zu v    and u  , and outer variables  and U , are compared with the results of DNS, experiments, and field measurements.

Validation by Non-Dimensional Inner Variables
To confirm that the WMLES approach is suitable for TBL flows at high Reynolds numbers, variations in mean velocity, Reynolds stress, and turbulent variation with height and Reynolds number are compared with theoretical results and experimental, field-measured data for the TBL flow.
The non-dimensional mean streamwise velocity profiles, shown in Figure 2. The simulated nondimensional mean velocities are displayed at different heights from case 1. The mean velocity profile agrees well with the logarithmic law and Qingtu Lake Observation Array (QLOA) field-measured data, indicating that the WMLES is highly accurate for TBL flow, despite the down shift in due to ground roughness.
Comparison of mean velocity profiles. The solid colored symbols indicate current simulation data; the black dashed line indicates U + = ln(z + )/0.41 + 5.0 for a smooth-wall wind profile; the black open upward-facing triangles show Surface layer Turbulence and Environmental Science Test (SLTEST) field measured data, Re τ = 6.28 × 10 5 ; the solid black circles indicate the Melbourne wind tunnel experiment data, Re τ = 1.801 × 10 4 ; the solid black squares represent the William B. Morgan Large Cavitation Channel experiment data acquired by the US Navy, Re τ = 6.867 × 10 4 ; the solid black leftward-facing triangles show the Princeton Superpipe experiment data, Re τ = 9.819 × 10 4 , the observational and experimental data obtained by Marusic et al. [33]; the open symbols present the QLOA field-measured data from Wang and Zheng [9]; the black open squares show Re τ = 1.63 × 10 6 ; the black open circles represent Re τ = 2.19 × 10 6 ; the black open leftward-facing triangles give Re τ = 2.95 × 10 6 ; the black open rightward-facing triangles provide Re τ = 3.14 × 10 6 ; and the black open stars show Re τ = 4.19 × 10 6 . Figure 3 shows the non-dimensional Reynolds shear-stress profiles for case 1. These agree well with theoretical predictions [34,35] for canonical TBL flow, the field-observation data of Hutchins et al. [36] Energies 2020, 13, 659 7 of 26 at Re τ = 7.7 × 10 5 , and the QLOA field measured data with higher Reynolds number of Wang and Zheng [9].
Energies 2020, 13 In addition, 2 / uu u  is plotted against the inner-scaled height, as shown in Figure 4, including the field measurement data from the SLTEST site [36,38], laboratory data [39,40], and the similarity formulation presented in Marusic and Kunkel [41]. This reasonable agreement suggests that the simulated data agree well with the theoretical formulation [34] and previous measurements at the SLTEST. The magnitude of the normalized streamwise component of normal stress is strongly dependent on the Reynolds number, except for the sublayer. Specifically, the normalized streamwise component of normal stress, simulated at /~0.036 z   , was compared to the results of QLOA [9], SLTEST [36], and wind tunnel experiments [40,[42][43][44] in Figure 5. The current simulation result falls within the QLOA results because the Reynolds number for case 1 is located within the five cases of Wang and Zheng [9], indicating an approximate log-linear relationship between 2 / uu u  and the Reynolds number.    The solid colored symbols are current simulation data, the dashed line is Chauhan [34] for Re τ = 2.33 × 10 6 , the solid line represents Re τ = 3.76 × 10 6 , the solid black upward-facing triangles exhibit the wind tunnel results of Graaff and Eaton [37] at Re τ = 1.35 × 10 3 , the solid black downward-facing triangles give the results of the field observations of Hutchins et al. [36] at Re τ = 7.7 × 10 5 , the open symbols show the field-measured QLOA data of Wang and Zheng [9], the black open squares represent Re τ = 1.63 × In addition, uu/u 2 τ is plotted against the inner-scaled height, as shown in Figure 4, including the field measurement data from the SLTEST site [36,38], laboratory data [39,40], and the similarity formulation presented in Marusic and Kunkel [41]. This reasonable agreement suggests that the simulated data agree well with the theoretical formulation [34] and previous measurements at the SLTEST. The magnitude of the normalized streamwise component of normal stress is strongly dependent on the Reynolds number, except for the sublayer. Specifically, the normalized streamwise component of normal stress, simulated at z/δ ≈ ∼ 0.036, was compared to the results of QLOA [9], SLTEST [36], and wind tunnel experiments [40,[42][43][44] in Figure 5. The current simulation result falls within the QLOA results because the Reynolds number for case 1 is located within the five cases of Wang and Zheng [9], indicating an approximate log-linear relationship between uu/u 2 τ and the Reynolds number. , was compared to the results of QLOA [9], SLTEST [36], and wind tunnel experiments [40,[42][43][44] in Figure 5. The current simulation result falls within the QLOA results because the Reynolds number for case 1 is located within the five cases of Wang and Zheng [9], indicating an approximate log-linear relationship between 2 / uu u  and the Reynolds number.    [38] at Re τ ≈ 3.1 × 10 6 and 3.8 × 10 6 , respectively, the black pluses exhibit the result of Hutchins et al. [36] at Re τ = 7.7 × 10 5 , the times and open stars represent Metzger et al. [39] at Re τ ≈ 8.3 × 10 5 , the open upper triangles show the result of Hutchins et al. [40], Re τ = 2.8 × 10 3 , 7.3 × 10 3 , and 1.903 × 10 4 , the lines represent the similarity formulations proposed by Marusic and Kunkel [41] calculated at Re τ = 7.7 × 10 5 (solid black line), Re τ = 3.2 × 10 6 (black dash line), Re τ = 3.8 × 10 6 (black dotted line).    [42], Knobloch and Fernholz [43], Hutchins et al. [40], and Kulandaivelu and Marusic [44]; the solid black five-pointed stars represent the field observations of Hutchins et al. [36]; and the solid black squares are the field observations of Wang and Zheng [9].

Validation by Non-Dimensional Outer Variables
In the outer layer of TBL, the flow is assumed to be independent of viscosity and to depend on the global characteristics of the flow, as represented by  , U .
The Reynolds shear stress is shown in Figure 6 against the outer-scaled coordinates / z  . It is clear that the simulated result of case 1 is in a good agreement with the experimental results of Graaff and Eaton [37]; most importantly, around / 0.02 z   , a trend in the shift of Reynolds shear stress is also captured by this simulation, which indicates the accuracy and high precision of the WMLES adopted in this study.
Normalized streamwise component of normal stress is shown in Figure 7 against outer-scaled coordinates / z  . The simulated result of case 1 agrees well with the experimental results of Graaff and Eaton [37]; like the Reynolds shear stress, a normalized streamwise component of normal stress also shows a shift in the trend with increasing height.   [40], and Kulandaivelu and Marusic [44]; the solid black five-pointed stars represent the field observations of Hutchins et al. [36]; and the solid black squares are the field observations of Wang and Zheng [9].

Validation by Non-Dimensional Outer Variables
In the outer layer of TBL, the flow is assumed to be independent of viscosity and to depend on the global characteristics of the flow, as represented by δ, U.
The Reynolds shear stress is shown in Figure 6 against the outer-scaled coordinates z/δ. It is clear that the simulated result of case 1 is in a good agreement with the experimental results of Graaff and Eaton [37]; most importantly, around z/δ = 0.02, a trend in the shift of Reynolds shear stress is also captured by this simulation, which indicates the accuracy and high precision of the WMLES adopted in this study. Normalized streamwise component of normal stress is shown in Figure 7 against outer-scaled coordinates z/δ. The simulated result of case 1 agrees well with the experimental results of Graaff and Eaton [37]; like the Reynolds shear stress, a normalized streamwise component of normal stress also shows a shift in the trend with increasing height. Figure 8 shows the mean velocity defect profiles according to the outer coordinates. The profile of case 1 is in good agreement with the experiment results of Graaff and Eaton [37] and the DNS results of Lee and Sung [8] throughout the outer region of the boundary layer, albeit at different Reynolds numbers.       Figure 9 shows the wavenumber-normalized spectra of streamwise velocity 2 / u S u   for all five cases at different heights. Taylor's hypothesis of spatial-temporal transformations was adopted to calculate the wavenumber spectrum, to smooth the spectra, a multi-point averaging method was used, following Balasubramaniam [45]. As Figure 9a-c shows, for 2 k  (wavelengths larger than ), each spectral curve has an individual inflection point where the slope of the curves conspicuously changes. The inflection point represents a spectral peak in low-wavenumber region of the pre-multiplied spectra (as shown in Figure 10) associated with the VLSMs [3,9,46]. Thus, it can also be demonstrated that the VLSMs exist in the first three cases for the current TBL study. While, for cases 4 and 5, individual inflection points exist for 2 k  , which means that VLSMs does not exist in these two cases. Above all, with where roughness increases, the scale of the VLSMs (the width of −1 region) decreases or even disappears. We predict that the reason for this may be associated with roughness elements of similar or larger sizes, which undoubtedly disrupt near-wall quasi-streamwise streaks. This phenomenon is identical to what is shown in Figure 1 from Volino et al. [17], although it is not mentioned as such in the report of that study.   Figure 9 shows the wavenumber-normalized spectra of streamwise velocity S u /u 2 τ δ for all five cases at different heights. Taylor's hypothesis of spatial-temporal transformations was adopted to calculate the wavenumber spectrum, to smooth the spectra, a multi-point averaging method was used, following Balasubramaniam [45]. As Figure 9a-c shows, for kδ < 2 (wavelengths larger than 3δ), each spectral curve has an individual inflection point where the slope of the curves conspicuously changes. The inflection point represents a spectral peak in low-wavenumber region of the pre-multiplied spectra (as shown in Figure 10) associated with the VLSMs [3,9,46]. Thus, it can also be demonstrated that the VLSMs exist in the first three cases for the current TBL study. While, for cases 4 and 5, individual inflection points exist for kδ > 2, which means that VLSMs does not exist in these two cases. Above all, with where roughness increases, the scale of the VLSMs (the width of −1 region) decreases or even disappears. We predict that the reason for this may be associated with roughness elements of similar or larger sizes, which undoubtedly disrupt near-wall quasi-streamwise streaks. This phenomenon is identical to what is shown in Figure 1 from Volino et al. [17], although it is not mentioned as such in the report of that study.

Velocity Spectra Characteristics under Different Surface Roughnesses
Streamwise Velocity Spectra Figure 9 shows the wavenumber-normalized spectra of streamwise velocity 2 / u S u   for all five cases at different heights. Taylor's hypothesis of spatial-temporal transformations was adopted to calculate the wavenumber spectrum, to smooth the spectra, a multi-point averaging method was used, following Balasubramaniam [45]. As Figure 9a-c shows, for 2 k  (wavelengths larger than ), each spectral curve has an individual inflection point where the slope of the curves conspicuously changes. The inflection point represents a spectral peak in low-wavenumber region of the pre-multiplied spectra (as shown in Figure 10) associated with the VLSMs [3,9,46]. Thus, it can also be demonstrated that the VLSMs exist in the first three cases for the current TBL study. While, for cases 4 and 5, individual inflection points exist for 2 k  , which means that VLSMs does not exist in these two cases. Above all, with where roughness increases, the scale of the VLSMs (the width of −1 region) decreases or even disappears. We predict that the reason for this may be associated with roughness elements of similar or larger sizes, which undoubtedly disrupt near-wall quasi-streamwise streaks. This phenomenon is identical to what is shown in Figure 1 from Volino et al. [17], although it is not mentioned as such in the report of that study.   Figure 10 shows the pre-multiplied spectra for all five roughness cases. There is a distinct peak in the spectra for all the heights in Figure 10a-c in the low-wavenumber region that corresponds to VLSMs; another peak appears in high-wavenumber region that corresponds to the LSMs in the spectra for all heights in Figure 10a-e.
In addition, the product of wavenumber and power spectra u kS decreases with height throughout all of the spectra under various ground roughness, which is consistent with the result of the turbulent-boundary layer [2,46], and the pipe-flow result [4]. This phenomenon can be explained by "bottom-up" mechanism, and VLSMs weaken in the wall-normal direction, due to the TKE produced in near-wall region gradually dissipates.
The slope between the first and second peaks of the streamwise velocity spectra is also analyzed. This is governed by the transformation scaling law of VLSMs and LSMs (cases 1, 2, and 3) or the transformation scaling law of LSMs and the small-scale turbulence vortex (cases 4 and 5); instead of observing jumps from −1 to −5/3, we observe a gradual change as height increases. With increasing roughness, the scaling law grows, and the transformation rate from −1 to −5/3 increases, with the roughness increasing as well, while the VLSMs go extinct. Another interesting phenomenon observed here is the scaling law in cases 4 and 5, showing height increases of −5/3. For other cases, the scaling law is located between the values −1 and −5/3, with the smoother terrain yielding the larger scaling law.
Based on this finding, a new concept could be proposed, according to which the presence or absence of VLSMs cannot simply be defined by the value −1 or −5/3; instead, the scale factor /   (the length of VLSMs and LSMs to the boundary-layer height) and the scaling law should be combined to illustrate the location of the vortex structure in VLSMs, LSMs, or small-scale turbulence vortex regions.   Figure 10 shows the pre-multiplied spectra for all five roughness cases. There is a distinct peak in the spectra for all the heights in Figure 10a-c in the low-wavenumber region that corresponds to VLSMs; another peak appears in high-wavenumber region that corresponds to the LSMs in the spectra for all heights in Figure 10a-e.
In addition, the product of wavenumber and power spectra kS u decreases with height throughout all of the spectra under various ground roughness, which is consistent with the result of the turbulent-boundary layer [2,46], and the pipe-flow result [4]. This phenomenon can be explained by "bottom-up" mechanism, and VLSMs weaken in the wall-normal direction, due to the TKE produced in near-wall region gradually dissipates.
The slope between the first and second peaks of the streamwise velocity spectra is also analyzed. This is governed by the transformation scaling law of VLSMs and LSMs (cases 1, 2, and 3) or the transformation scaling law of LSMs and the small-scale turbulence vortex (cases 4 and 5); instead of observing jumps from −1 to −5/3, we observe a gradual change as height increases. With increasing roughness, the scaling law grows, and the transformation rate from −1 to −5/3 increases, with the roughness increasing as well, while the VLSMs go extinct. Another interesting phenomenon observed here is the scaling law in cases 4 and 5, showing height increases of −5/3. For other cases, the scaling law is located between the values −1 and −5/3, with the smoother terrain yielding the larger scaling law.
Based on this finding, a new concept could be proposed, according to which the presence or absence of VLSMs cannot simply be defined by the value −1 or −5/3; instead, the scale factor λ/δ (the length of VLSMs and LSMs to the boundary-layer height) and the scaling law should be combined to illustrate the location of the vortex structure in VLSMs, LSMs, or small-scale turbulence vortex regions. The k decreases at the first peak (meaning 0 ( ) k ) with increasing height; when the height reaches to a certain extent, 0 ( ) k begins to attain a constant value or even to increase, as determined by the characteristics of the VLSMs and is consistent with the literature [18,46]. These two phenomena can be seen more clearly in Figure 11, which plots the relationship between 0 ( ) k and the dimensionless height / z  . In Figure 11, the result of case 1 agrees with the field-measured data, as a result of the similarity of terrain between this and the QLOA site, as shown in Figure 2. With increases in roughness, cases 2 and 3 reach a certain degree of agreement with field data. While, for cases 4 and 5, the results diverge from the field data, implying that the surface roughness has a great influence on the structure of VLSMs. As we know, there is also a certain deviation between the experimental value and the field-measured value, which may be due to the difference in Reynolds number between them. The kδ decreases at the first peak (meaning (kδ) 0 ) with increasing height; when the height reaches to a certain extent, (kδ) 0 begins to attain a constant value or even to increase, as determined by the characteristics of the VLSMs and is consistent with the literature [18,46]. These two phenomena can be seen more clearly in Figure 11, which plots the relationship between (kδ) 0 and the dimensionless height z/δ.
In Figure 11, the result of case 1 agrees with the field-measured data, as a result of the similarity of terrain between this and the QLOA site, as shown in Figure 2. With increases in roughness, cases 2 and 3 reach a certain degree of agreement with field data. While, for cases 4 and 5, the results diverge from the field data, implying that the surface roughness has a great influence on the structure of VLSMs. As we know, there is also a certain deviation between the experimental value and the field-measured value, which may be due to the difference in Reynolds number between them.

Wall-Normal Velocity Spectra Characteristics
Wall-normal velocity fluctuation spectra are shown in Figure 12. The trend for wall-normal location is similar to the streamwise velocity fluctuations, with the exception of the crossover phenomenon, which is consistent with the finding of Guala et al. [4]. For the left part of the crossover point, energy density increases from the wall to a wall-normal location, and then it decreases again; for the right part of the crossover point, the energy density always decreases. In addition, it is clear that as the roughness increases, the crossover point moves to a large wavenumber, which indicates that the rougher the terrain, the less the vortex interacts at different heights. Comparison with streamwise velocity spectra, the wall-normal velocity spectra has lower values, which reflects the difference between the two components in the total energy, and the ultimate difference is close to one order of magnitude.

Wall-Normal Velocity Spectra Characteristics
Wall-normal velocity fluctuation spectra are shown in Figure 12. The trend for wall-normal location is similar to the streamwise velocity fluctuations, with the exception of the crossover phenomenon, which is consistent with the finding of Guala et al. [4]. For the left part of the crossover point, energy density increases from the wall to a wall-normal location, and then it decreases again; for the right part of the crossover point, the energy density always decreases. In addition, it is clear that as the roughness increases, the crossover point moves to a large wavenumber, which indicates that the rougher the terrain, the less the vortex interacts at different heights. Comparison with streamwise velocity spectra, the wall-normal velocity spectra has lower values, which reflects the difference between the two components in the total energy, and the ultimate difference is close to one order of magnitude.
Energies 2020, 13, x FOR PEER REVIEW 13 of 26 Figure 11. Variation in the peak wavenumber of VLSMs with height. The solid colored circle is the current simulation data, the open symbols show the field-measured data of Wang and Zheng [9], and the solid black symbols represent the experimental data of Valliviki et al. [46].

Wall-Normal Velocity Spectra Characteristics
Wall-normal velocity fluctuation spectra are shown in Figure 12. The trend for wall-normal location is similar to the streamwise velocity fluctuations, with the exception of the crossover phenomenon, which is consistent with the finding of Guala et al. [4]. For the left part of the crossover point, energy density increases from the wall to a wall-normal location, and then it decreases again; for the right part of the crossover point, the energy density always decreases. In addition, it is clear that as the roughness increases, the crossover point moves to a large wavenumber, which indicates that the rougher the terrain, the less the vortex interacts at different heights. Comparison with streamwise velocity spectra, the wall-normal velocity spectra has lower values, which reflects the difference between the two components in the total energy, and the ultimate difference is close to one order of magnitude.

Effects of VLSMs and LSMs on Turbulent Kinetic Energy
One of the principal questions addressed in this study is how far VLSMs and LSMs contribute to TKE. We will examine the distribution of TKE as a function of height.
The TKE of VLSMs and the corresponding energy fraction are plotted against height using a low-pass filter with a cutoff length of , as shown in Figure 13a,b), respectively. Similarly, bandpass filtering with cutoff lengths of 0.3 3

  
and 3    are used to demonstrate the energy fraction of the LSMs, as shown in Figure 13(c and d), respectively. The total kinetic energy of VLSMs of both case 1 in the present study and the experimental data of Balakumar and Adrian [2] slightly decreases with increases in height. However, the total kinetic energy of VLSMs increases with height increases as in Wang and Zheng [9], which may relate to fact that the case investigated by Wang and Zheng [9] found energy increases with height in the low-wavenumber region of the wind spectra, unlike what was found in other studies [2,46]. Case 1 shows a close agreement with the experimental data of Balakumar and Adrian [2], given the increasing ground roughness and nondimensional decreases in TKE. The reason for this is found in the fact that VLSMs decrease or even disappear as roughness increases. It is worth noticing that Figure 13 also shows a result found in cases 4 and 5 that does not exist in VLSMs. It simply is the mathematical expression of the region of VLSMs and LSMs, as shown in Figure 10.
In addition, the energy fraction of the VLSMs in the present study increases with height, which tends to close agreement with the wind-tunnel measurements of Balakumar and Adrian [2], as shown in Figure 13(b); however, little deviation from the field-measured data of Wang and Zheng [9] appears, due to the spectrum characteristics mentioned above. The phenomenon, shown in Figure 13 indicates, that contrary to most other studies [4,8,47], which have generally found that VLSMs contribute 40-50% of the total TKE, the present study shows that VLSMs contribute up to 40% only at the highest height, shown in Figure 13(b). In Figure 13(c), it is seen that in LSMs ( 0.3 3      , Guala et al. [4]), the energy fraction is larger than the field-measured data; but in Figure 13(d), the energy fraction of the LSMs ( 3      , Lee et al. [48]) is similar to the field measured data. It is clear that the energy fraction of the LSMs is independent of ground roughness, possibly owing to the fact that the size of the element roughness was adopted in the present study only to disrupt the VLSMs, yielding no effect on LSMs, at least from the viewpoint of energy fraction.
The reason why the energy fraction of LSMs is found to be large in our study can be explained as follows. Because the WMLES approach is used, the smallest scale falls in the inertial subrange, while some or even most of the inertial subrange remains unresolved. It can be seen from the wind spectra shown in Figure 10 that it quickly falls into the dissipative region. Therefore, the total energy in this study is small. Due to the relative limitations of the computational domain, the structure of VLSMs is not completely contained, so that the simulation of the energy of VLSMs is less than that of the real atmosphere or the sufficient computational-domain case. This might also explain why the experimental data of Balakumar and Adrian [2] are lower than the field-measured of from Wang and Zheng [9], as shown in Figure 13b, the experimental domain being also limited to some extent. The

Effects of VLSMs and LSMs on Turbulent Kinetic Energy
One of the principal questions addressed in this study is how far VLSMs and LSMs contribute to TKE. We will examine the distribution of TKE as a function of height.
The TKE of VLSMs and the corresponding energy fraction are plotted against height using a low-pass filter with a cutoff length of 3δ, as shown in Figure 13a,b, respectively. Similarly, band-pass filtering with cutoff lengths of 0.3δ − 3δ and δ − 3δ are used to demonstrate the energy fraction of the LSMs, as shown in Figure 13c,d, respectively. The total kinetic energy of VLSMs of both case 1 in the present study and the experimental data of Balakumar and Adrian [2] slightly decreases with increases in height. However, the total kinetic energy of VLSMs increases with height increases as in Wang and Zheng [9], which may relate to fact that the case investigated by Wang and Zheng [9] found energy increases with height in the low-wavenumber region of the wind spectra, unlike what was found in other studies [2,46]. Case 1 shows a close agreement with the experimental data of Balakumar and Adrian [2], given the increasing ground roughness and nondimensional decreases in TKE. The reason for this is found in the fact that VLSMs decrease or even disappear as roughness increases. It is worth noticing that Figure 13 also shows a result found in cases 4 and 5 that does not exist in VLSMs. It simply is the mathematical expression of the region of VLSMs and LSMs, as shown in Figure 10.
In addition, the energy fraction of the VLSMs in the present study increases with height, which tends to close agreement with the wind-tunnel measurements of Balakumar and Adrian [2], as shown in Figure 13b; however, little deviation from the field-measured data of Wang and Zheng [9] appears, due to the spectrum characteristics mentioned above. The phenomenon, shown in Figure 13 indicates, that contrary to most other studies [4,8,47], which have generally found that VLSMs contribute 40-50% of the total TKE, the present study shows that VLSMs contribute up to 40% only at the highest height, shown in Figure 13b. In Figure 13c, it is seen that in LSMs (0.3δ < λ < 3δ, Guala et al. [4]), the energy fraction is larger than the field-measured data; but in Figure 13d, the energy fraction of the LSMs (δ < λ < 3δ, Lee et al. [48]) is similar to the field measured data. It is clear that the energy fraction of the LSMs is independent of ground roughness, possibly owing to the fact that the size of the element roughness was adopted in the present study only to disrupt the VLSMs, yielding no effect on LSMs, at least from the viewpoint of energy fraction.
The reason why the energy fraction of LSMs is found to be large in our study can be explained as follows. Because the WMLES approach is used, the smallest scale falls in the inertial subrange, while some or even most of the inertial subrange remains unresolved. It can be seen from the wind spectra shown in Figure 10 that it quickly falls into the dissipative region. Therefore, the total energy in this study is small. Due to the relative limitations of the computational domain, the structure of VLSMs is not completely contained, so that the simulation of the energy of VLSMs is less than that of the real atmosphere or the sufficient computational-domain case. This might also explain why the experimental data of Balakumar and Adrian [2] are lower than the field-measured of from Wang and Zheng [9], as shown in Figure 13b, the experimental domain being also limited to some extent.
The scale of LSMs, ranging between 0.3δ − 3δ, can be fully resolved. In summary, the energy of VLSMs is relatively smaller than the field-measured data, while the energy of LSMs is larger.

 
 , can be fully resolved. In summary, the energy of VLSMs is relatively smaller than the field-measured data, while the energy of LSMs is larger. Figure 13. Variations in the total energy (a) and energy fraction of VLSMs (b) and LSMs (c, Solid colored circles are current simulation data, the open symbols represent the field-measured data of Wang and Zheng [9], and the solid black symbols give the experimental data of Balakumar and Adrian [2].

Flow Field Correlation Analysis
To intuitively visualize the structure of VLSMs, the coefficients of spatial correlation of the streamwise fluctuating velocity between the spatial and reference locations ( ( , , ) (2500 ,500 , 35 ) x y z m m m  ) were calculated. Figure 14 shows the contours of the spatial correlation coefficient for five types of roughness cases. The contour level of the correlation coefficient indicating the coherent structure edge was set to be 0.05, following Hutchins and Marusic [7]. As the figure shows, the length scale of the streamwise coherent structure satisfies / Along the streamwise direction, the correlation contour lines are tilted, which indicates the inclination of turbulent structures that observed by Marusic and Heuer [49] and Wang and Zheng [9] is also present in the current TBL simulation at higher Reynolds numbers. The inclination angle of the two-point correlation of the fluctuating velocity is related to the average extent of the inclination of the structures of the VLSMs or LSMs. For the present simulation, the inclination angles in Figure  14 199  10 7 ). The results of case 1 agree closely with those found by Wang and Zheng [9], who indicated that the inclination angles of field-measured data were 12.14° for 3.14  10 6 , while for 1.63  10 6 , the angle is 20.35°. The difference among all five cases is comparable to the range reported in the literature for smooth-and rough-wall boundary layers. Christensen and Wu [50] found an inclination angle of 11° for their smooth-wall channel flow. Head and Bandyopadhyay [51] observed inclination angles between 15-20°. Christensen and Adrian [52] reported 12-13°. Adrian et al. [53] found an inclination angle of 12°. Tomkins and Adrian [54] found inclination angles between 10-20°. Nakagawa and  Figure 13. Variations in the total energy (a) and energy fraction of VLSMs (b) and LSMs (c, 0.3δ < λ < 3δ; d, δ < λ < 3δ ) with height. Solid colored circles are current simulation data, the open symbols represent the field-measured data of Wang and Zheng [9], and the solid black symbols give the experimental data of Balakumar and Adrian [2].

Flow Field Correlation Analysis
To intuitively visualize the structure of VLSMs, the coefficients of spatial correlation of the streamwise fluctuating velocity between the spatial and reference locations ((x, y, z) = (2500m, 500m, 35m)) were calculated. Figure 14 shows the contours of the spatial correlation coefficient for five types of roughness cases. The contour level of the correlation coefficient indicating the coherent structure edge was set to be 0.05, following Hutchins and Marusic [7]. As the figure shows, the length scale of the streamwise coherent structure satisfies λ/δ > 3 for the first three cases, indicating the existence of the VLSMs. For cases 4 and 5, the 0.3 < λ/δ < 3 range only bears on the existence of the LSMs. For increasing roughness, the scale of VLSMs and LSMs decreases, indicating that the generation and development of VLSMs and LSMs are affected by the bottom terrain. It also indirectly illustrates that the bottom-up mechanism controls the generation and development of VLSMs.
Along the streamwise direction, the correlation contour lines are tilted, which indicates the inclination of turbulent structures that observed by Marusic and Heuer [49] and Wang and Zheng [9] is also present in the current TBL simulation at higher Reynolds numbers. The inclination angle of the two-point correlation of the fluctuating velocity is related to the average extent of the inclination of the structures of the VLSMs or LSMs. For the present simulation, the inclination angles in Figure 14 are 11.31 • for case 1 (Re τ = 3.538 × 10 6 ), 14.04 • for case 2 (Re τ = 4.894 × 10 6 ), 16.70 • for case 3 (Re τ = 6.559 × 10 6 ), 14.04 • for case 4 (Re τ = 9.927 × 10 6 ), and 19.29 • for case 5 (Re τ = 1.199 × 10 7 ). The results of case 1 agree closely with those found by Wang and Zheng [9], who indicated that the inclination angles of field-measured data were 12.14 • for Re τ = 3.14 × 10 6 , while for Re τ = 1.63 × 10 6 , the angle is 20.35 • . The difference among all five cases is comparable to the range reported in the literature for smooth-and rough-wall boundary layers. Christensen and Wu [50] found an inclination angle of 11 • for their smooth-wall channel flow. Head and Bandyopadhyay [51] observed inclination angles between 15-20 • . Christensen and Adrian [52] reported 12-13 • . Adrian et al. [53] found an inclination angle of 12 • . Tomkins and Adrian [54] found inclination angles between 10-20 • . Nakagawa and Hanratty [55] found an angle of 9 • on the rough wall in a channel flow, compared to smooth-wall values of 6-8 • . In general, the inclination angle of large-scale coherent structures increases with roughness. Hanratty [55] found an angle of 9° on the rough wall in a channel flow, compared to smooth-wall values of 6-8°. In general, the inclination angle of large-scale coherent structures increases with roughness. Figure 14. Two-point correlation contours of streamwise-fluctuating velocity.
In addition, the streamwise, wall normal, and spanwise extent of fluctuating-velocity two-point correlations are shown in Figure 15. The distance, uu Lx , twice the distance from the most downstream location to the self-correlation peak on two-point correlations of the fluctuating streamwise velocity equal to 0.5 contour, as defined in Christensen and Wu [50]. Volino et al. [17] found that the results for rough and smooth walls agree well with a value of about In addition, the streamwise, wall normal, and spanwise extent of fluctuating-velocity two-point correlations are shown in Figure 15. The distance, Lx uu , twice the distance from the most downstream location to the self-correlation peak on two-point correlations of the fluctuating streamwise velocity equal to 0.5 contour, as defined in Christensen and Wu [50]. Volino et al. [17] found that the results for rough and smooth walls agree well with a value of about Lx uu /δ = 0.65 between z/δ = 0.1 and 0.6 with their experimental data. Close to the wall, some differences are visible. In Figure 15a, the result of case 1 is much smaller than that found by Volino et al. [17]; however, the field observation data of Wang and Zheng [9] with solid symbols are shown in Figure 15a, for Re τ = 3.14 × 10 6 measured data, almost same as in case 1 (Re τ = 3.553 × 10 6 ). This indicates that the results of simulation are correct. This can be explained by the fact that field-measured and simulation data have a similarity ratio comparable to the experimental data. In Figure 15a, as roughness increases, Lx uu /δ decreases, which agrees with the results of Krogstad and Antonia [56], who found Lx uu /δ was about 50% lower on their rough wall. The wall-normal extent of the two-point correlations of the fluctuating velocity, Lz uu , is determined based on the wall-normal distance between the points farthest and closest from the surface at a 0.5 contour. As the contours merge with the wall, reliable estimates for Lz uu cannot be obtained for z/δ < 0.2. The conclusion is similar to the Lx uu , as shown in Figure 15b. The spanwise extent Ly uu shown in Figure 15c, with increasing roughness, the Ly uu /δ range from 0.2 to 0.1, which agrees with the DNS result of Lee and Sung [8], who concluded that the streamwise negative motions meander in the downstream direction with a characteristic width of approximately 0.1δ − 0.2δ in the spanwise direction. Comparing Figure 15a,b, the ratio of Lx uu /Lz uu obtained in Figure 15d, the value is roughly 2.5 for both the experimental data of Volino et al. [17] and the present simulation cases, also consistent with Nakagawa and Hanratty [55], and the smooth-wall results of Krogstad and Antonia [56].
Energies 2020, 13, x FOR PEER REVIEW 17 of 26 result of case 1 is much smaller than that found by Volino et al. [17]; however, the field observation data of Wang and Zheng [9] with solid symbols are shown in Figure 15(a), for 3.14  10 6 measured data, almost same as in case 1 ( 3.553  10 6 ). This indicates that the results of simulation are correct. This can be explained by the fact that field-measured and simulation data have a similarity ratio comparable to the experimental data. In Figure 15(a), as roughness increases, / uu Lx  decreases, which agrees with the results of Krogstad and Antonia [56], who found was about 50% lower on their rough wall. The wall-normal extent of the two-point correlations of the fluctuating velocity, uu Lz , is determined based on the wall-normal distance between the points farthest and closest from the surface at a 0.5 contour. As the contours merge with the wall, reliable estimates for cannot be obtained for / 0.2 z   . The conclusion is similar to the uu Lx , as shown in Figure 15 Figure 15(d), the value is roughly 2.5 for both the experimental data of Volino et al. [17] and the present simulation cases, also consistent with Nakagawa and Hanratty [55], and the smooth-wall results of Krogstad and Antonia [56].  [17], the circle symbols represent smooth wall condition, and squares represent roughness wall condition, the field measured data of Wang and Zheng [9], the anise star is R e   3.14  10 6 , and five star is R e   1.63  10 6 which show in (a). Figure 16 shows visualizations of the VLSMs or LSMs structures in the instantaneous flow fields for z = 30 m in the xy plane, both the VLSMs and LSMs were visualized based on the streamwiseelongated negative velocity, which is consistent with the results of previous studies [6,8,[57][58][59]. The direction of flow is from left to right. Figure 16a [17], the circle symbols represent smooth wall condition, and squares represent roughness wall condition, the field measured data of Wang and Zheng [9], the anise star is Re τ = 3.14 × 10 6 , and five star is Re τ = 1.63 × 10 6 which show in (a). Figure 16 shows visualizations of the VLSMs or LSMs structures in the instantaneous flow fields for z = 30 m in the xy plane, both the VLSMs and LSMs were visualized based on the streamwise-elongated negative velocity, which is consistent with the results of previous studies [6,8,[57][58][59]. The direction of flow is from left to right. Figure 16a-c show several apparent very long negative (dark color) motions meandering in the streamwise direction, which in the spanwise direction are flanked by positive streamwise fluctuation (white color), and these motions often extend more than 10δ in streamwise direction, and the widths is around 0.1δ − 0.2δ. This phenomenon is in agreement with TBL and turbulent channel DNS studies [7,8]. For Figure 16d Figure 16c in the xz plane, and the dashed line in Figure  17 corresponds to the horizontal plane height of Figure 16. In Figure 16c, line A locates at the very long streamwise structure, line B locates at the intermediate long streamwise structure in order to see the inclination angle of VLSMs, and line C locates at the highvelocity region. Comparing these three lines streamwise fluctuation in Figure 17, there are inclined at angles of 18.26° of line B, very close to the spatial correlation result, 16.70°, which is calculated by the reference point of streamwise direction. Note that the inclination angle is 10°-15° in an almost linear type [8], it also can be seen in the present study, the line A has nearly same inclination angle. Moreover, the flow patterns of these aligned packets are consistent with the observations [53]. For line C, which locates at the high-velocity region, the inclination angles at A and B do not exist at line C, this indicates that LSMs and VLSMs structures located at low momentum region, and comes from the near wall aligned packets, as described by "bottom-up" mechanism.   Figure 16c in the xz plane, and the dashed line in Figure 17 corresponds to the horizontal plane height of Figure 16. In Figure 16c, line A locates at the very long streamwise structure, line B locates at the intermediate long streamwise structure in order to see the inclination angle of VLSMs, and line C locates at the high-velocity region. Comparing these three lines streamwise fluctuation in Figure 17, there are inclined at angles of 18.26 • of line B, very close to the spatial correlation result, 16.70 • , which is calculated by the reference point of streamwise direction. Note that the inclination angle is 10 • -15 • in an almost linear type [8], it also can be seen in the present study, the line A has nearly same inclination angle. Moreover, the flow patterns of these aligned packets are consistent with the observations [53]. For line C, which locates at the high-velocity region, the inclination angles at A and B do not exist at line C, this indicates that LSMs and VLSMs structures located at low momentum region, and comes from the near wall aligned packets, as described by "bottom-up" mechanism.

The Generation and Development of VLSMs
The results of Section 4.4 illustrate the presence of VLSMs and their associate with lowmomentum regions. In this section, we address the creation and extinction procession of VLSMs in the flow field. As we know, the hairpin packet model could give a reasonable explanation for logarithmic layer features, for example, quasi-streamwise vortices and bursting processes. An idealized conceptual model Based on the hairpin packet paradigm and in the outer layer, most features of the coherent structures, an idealized conceptual model is built [53]. In contrast, Kim and Adrian [3] suggested that VLSMs is not a new type of turbulent motion. They conjectured that the alignments of hairpin packets could form VLSMs; however, the created and extincted of VLSMs in TBL is little to know. Thus, to obtain spatial information about the formation and disappearance of VLSMs is our object.
The mechanism of Kim and Adrian [3] means that in the streamwise direction continually spawning new hairpins to create packets. After growing, the packets become larger and longer, and then probably merge with adjacent packets or broken up into small structures. In this study, a representative case of merging and extincting of vortex packets is discussed in detail, which VLSMs showed in the dashed red rectangle of Figure 16c Figure 18a, two adjacent low-momentum regions of P1 and P2 are travelling downstream, then they merge together from Figure 18f, the similar phenomenon is observed at dashed line D2. For dashed line E, the P2 and P3 are merge together, as they move downstream, and the low-momentum regions separate in the streamwise direction, as shown in Figure 18j. For dashed line F1, two adjacent low-momentum regions of P3 and P4 are travelling downstream, and the first stage is that they merge together from Figure 18a-e, then separate from each other again in the Figure  18j, the similar phenomenon is observed at dashed line F2. For dashed line G1, as part P5 and P6 move downstream, and the first stage is that they separate from Figure 18c and then merge together again in the Figure 18f, the similar phenomenon is observed at dashed line G2. The four cases of D1, D2, E, F1, F2 and G, represent three typical development process of adjacent hairpin vortices. Due to the low-momentum regions merge in streamwise direction, resulting in a longer structure at the streamwise scale of the flow field in Figure 18f, g. Another phenomenon is that although the streamwise growth of the VLSMs occurs; however, the spanwise length of VLSMs becomes thinner slightly.

The Generation and Development of VLSMs
The results of Section 4.4 illustrate the presence of VLSMs and their associate with low-momentum regions. In this section, we address the creation and extinction procession of VLSMs in the flow field. As we know, the hairpin packet model could give a reasonable explanation for logarithmic layer features, for example, quasi-streamwise vortices and bursting processes. An idealized conceptual model Based on the hairpin packet paradigm and in the outer layer, most features of the coherent structures, an idealized conceptual model is built [53]. In contrast, Kim and Adrian [3] suggested that VLSMs is not a new type of turbulent motion. They conjectured that the alignments of hairpin packets could form VLSMs; however, the created and extincted of VLSMs in TBL is little to know. Thus, to obtain spatial information about the formation and disappearance of VLSMs is our object.
The mechanism of Kim and Adrian [3] means that in the streamwise direction continually spawning new hairpins to create packets. After growing, the packets become larger and longer, and then probably merge with adjacent packets or broken up into small structures. In this study, a representative case of merging and extincting of vortex packets is discussed in detail, which VLSMs showed in the dashed red rectangle of Figure 16c Figure 18a, two adjacent low-momentum regions of P1 and P2 are travelling downstream, then they merge together from Figure 18f, the similar phenomenon is observed at dashed line D2. For dashed line E, the P2 and P3 are merge together, as they move downstream, and the low-momentum regions separate in the streamwise direction, as shown in Figure 18j. For dashed line F1, two adjacent low-momentum regions of P3 and P4 are travelling downstream, and the first stage is that they merge together from Figure 18a-e, then separate from each other again in the Figure 18j, the similar phenomenon is observed at dashed line F2. For dashed line G1, as part P5 and P6 move downstream, and the first stage is that they separate from Figure 18c and then merge together again in the Figure 18f, the similar phenomenon is observed at dashed line G2. The four cases of D1, D2, E, F1, F2 and G, represent three typical development process of adjacent hairpin vortices. Due to the low-momentum regions merge in streamwise direction, resulting in a longer structure at the streamwise scale of the flow field in Figure 18f,g. Another phenomenon is that although the streamwise growth of the VLSMs occurs; however, the spanwise length of VLSMs becomes thinner slightly.   Figure 18. Combining the Figures 18-20, which indicates that the development and extinction of VLSMs structure is a three-dimensional process. As can be seen in these figures, in the downstream direction, these packets are growing continuously and further away from the wall, part of the upstream packets go through the higher mean flow velocity. The lower parts slower than these part flows, which leads lift vortex away from the wall by a higher velocity, and in the higher part of the upstream packet causes in greater stretching. Consequently, in downstream the adjacent packets are merged to form VLSMs with a shallow angle. These processes are same with hairpin structure dynamics [5] and DNS results [8].
For the formation of VLSMs, present simulations describe one possible mechanism, however, in the future, the analysis of the general mechanism of VLSMs formation is needed, because the trigger of the formation of new packets still no definitive evidence.  Figure 18. Combining the Figures 18-20, which indicates that the development and extinction of VLSMs structure is a three-dimensional process. As can be seen in these figures, in the downstream direction, these packets are growing continuously and further away from the wall, part of the upstream packets go through the higher mean flow velocity. The lower parts slower than these part flows, which leads lift vortex away from the wall by a higher velocity, and in the higher part of the upstream packet causes in greater stretching. Consequently, in downstream the adjacent packets are merged to form VLSMs with a shallow angle. These processes are same with hairpin structure dynamics [5] and DNS results [8].
For the formation of VLSMs, present simulations describe one possible mechanism, however, in the future, the analysis of the general mechanism of VLSMs formation is needed, because the trigger of the formation of new packets still no definitive evidence.

Vortices Field Analysis
The packet merging process is presented schematically in Figure 17 of Tomkins and Adrian [54]. The similar conclusion was obtained in Figure 25 of Adrian, Meinhart and Tomkins [53]. This is the reason for the vorticity field analysis, which could show the VLSMs structure is formed by hairpin vortex or hairpin vortex package development. Figure 21 shows three-dimensional vortical structures that contoured by 2  [60] in the horizontal plane. It is clear that several highly elongated low-speed regions are shown in the flow field, and many hairpin-type vortices lie above each of them. Corresponding to the Figure 17, here Figure 22 is shown as the vertical plane to illustrate the negative velocity region is covered by vortices, while positive velocity is not. By comparing A, B, and C, it is clearly shown that vortices are more density for A than B than C, which locates in high-speed region-this is consistent with the preceding results. These flow patterns, together with three-dimensional vortical structures to form longer structures that erupt from the wall and grow towards the outer region, are consist with Lee and Sung [8].

Vortices Field Analysis
The packet merging process is presented schematically in Figure 17 of Tomkins and Adrian [54]. The similar conclusion was obtained in figure 25 of Adrian, Meinhart and Tomkins [53]. This is the reason for the vorticity field analysis, which could show the VLSMs structure is formed by hairpin vortex or hairpin vortex package development. Figure 21 shows three-dimensional vortical structures that contoured by λ 2 [60] in the horizontal plane. It is clear that several highly elongated low-speed regions are shown in the flow field, and many hairpin-type vortices lie above each of them. Corresponding to the Figure 17, here Figure 22 is shown as the vertical plane to illustrate the negative velocity region is covered by vortices, while positive velocity is not. By comparing A, B, and C, it is clearly shown that vortices are more density for A than B than C, which locates in high-speed region-this is consistent with the preceding results. These flow patterns, together with three-dimensional vortical structures to form longer structures that erupt from the wall and grow towards the outer region, are consist with Lee and Sung [8].

Vortices Field Analysis
The packet merging process is presented schematically in Figure 17 of Tomkins and Adrian [54]. The similar conclusion was obtained in Figure 25 of Adrian, Meinhart and Tomkins [53]. This is the reason for the vorticity field analysis, which could show the VLSMs structure is formed by hairpin vortex or hairpin vortex package development. Figure 21 shows three-dimensional vortical structures that contoured by 2  [60] in the horizontal plane. It is clear that several highly elongated low-speed regions are shown in the flow field, and many hairpin-type vortices lie above each of them. Corresponding to the Figure 17, here Figure 22 is shown as the vertical plane to illustrate the negative velocity region is covered by vortices, while positive velocity is not. By comparing A, B, and C, it is clearly shown that vortices are more density for A than B than C, which locates in high-speed region-this is consistent with the preceding results. These flow patterns, together with three-dimensional vortical structures to form longer structures that erupt from the wall and grow towards the outer region, are consist with Lee and Sung [8].

Conclusions
In this paper, the spatial features for LSMs and VLSMs in the TBL flow under different surface roughness at very high Reynolds number, O (10 6 -10 7 ), were investigated using the WMLES approach. The VLSMs in TBL were studied progressively from a single point, in correlation with a field analysis. The statistical properties of the mean and fluctuating velocity for the nearly smooth case were compared with those found in the literature, including field observations, laboratory studies, and DNSs, and a close agreement was found using the inner and outer variables. These results, thus, provide a useful database for the turbulence statistics of TBL.
Our analysis of the streamwise velocity spectra demonstrates the tendency is that with increasing roughness, the scale of VLSMs decreases or even disappears. The For wall-normal velocity spectra, the crossover phenomenon occurs, and as roughness increases, the crossover point moves to larger wavenumber.
The VLSMs and LSMs carry a considerable portion of the TKE. While, the energy fraction of VLSMs in the present study is a little smaller than the field-measured data, but it is similar to experimental data, which can be explained by spectral characteristics.
The spatial correlation coefficient pertaining to five different roughness cases were calculated. The results have reached a close agreement with the wind spectra. For the ratio of the streamwise extent and wall-normal extent of the two-point correlations of the fluctuating velocity, the nearconstant is 2.5. Moreover, the inclination angle of nearly smooth cases is consistent with field observation data at almost the same Reynolds number. The tendency of the inclination angle of different cases of roughness is similar to what has been shown in other studies.
Then, we conducted the instantaneous flow field analysis to illustrate the generation and development of VLSMs. It is clearly shown that both the VLSMs and LSMs were visualized based on the streamwise-elongated volumes associated with the negative streamwise velocity, and which are flanked on either side by similarly elongated positive streamwise velocity. In addition, the generation and extinction process of VLSMs is seen by three-dimensional flow field analysis.
In the end, Vortices field analysis is conducted to explain further the VLSMs locate in negative streamwise velocity region and each of them lies beneath many hairpin-type vortices.
Above all, in this study, the results at least indirectly support that the bottom-up mechanism is one of the origins of VLSMs. It is confirmed here that the WMLES method has the ability to investigate the spatial features of LSMs and VLSMs in the TBL flow at high Reynolds numbers at neutrally stratified TBL, and the applicability of these results to a real atmospheric boundary layer flow that stability variation will be investigated in the future study.

Conclusions
In this paper, the spatial features for LSMs and VLSMs in the TBL flow under different surface roughness at very high Reynolds number, O (10 6 -10 7 ), were investigated using the WMLES approach. The VLSMs in TBL were studied progressively from a single point, in correlation with a field analysis. The statistical properties of the mean and fluctuating velocity for the nearly smooth case were compared with those found in the literature, including field observations, laboratory studies, and DNSs, and a close agreement was found using the inner and outer variables. These results, thus, provide a useful database for the turbulence statistics of TBL.
Our analysis of the streamwise velocity spectra demonstrates the tendency is that with increasing roughness, the scale of VLSMs decreases or even disappears. The (kδ) 0 decreases with height in an approximately log-linear manner. While, when height increases to a certain extent, (kδ) 0 begins to reach a constant value or even to increase. In addition, kS u decreases with height throughout all spectra.
For wall-normal velocity spectra, the crossover phenomenon occurs, and as roughness increases, the crossover point moves to larger wavenumber.
The VLSMs and LSMs carry a considerable portion of the TKE. While, the energy fraction of VLSMs in the present study is a little smaller than the field-measured data, but it is similar to experimental data, which can be explained by spectral characteristics.
The spatial correlation coefficient pertaining to five different roughness cases were calculated. The results have reached a close agreement with the wind spectra. For the ratio of the streamwise extent and wall-normal extent of the two-point correlations of the fluctuating velocity, the near-constant is 2.5. Moreover, the inclination angle of nearly smooth cases is consistent with field observation data at almost the same Reynolds number. The tendency of the inclination angle of different cases of roughness is similar to what has been shown in other studies.
Then, we conducted the instantaneous flow field analysis to illustrate the generation and development of VLSMs. It is clearly shown that both the VLSMs and LSMs were visualized based on the streamwise-elongated volumes associated with the negative streamwise velocity, and which are flanked on either side by similarly elongated positive streamwise velocity. In addition, the generation and extinction process of VLSMs is seen by three-dimensional flow field analysis.
In the end, Vortices field analysis is conducted to explain further the VLSMs locate in negative streamwise velocity region and each of them lies beneath many hairpin-type vortices.
Above all, in this study, the results at least indirectly support that the bottom-up mechanism is one of the origins of VLSMs. It is confirmed here that the WMLES method has the ability to investigate the spatial features of LSMs and VLSMs in the TBL flow at high Reynolds numbers at neutrally stratified TBL, and the applicability of these results to a real atmospheric boundary layer flow that stability variation will be investigated in the future study.