A Numerical Study of the Effect of Vegetative Windbreak on Wind Erosion over Complex Terrain

: Wind erosion is a typical issue for stone carvings in northwest China caves, and windbreaks such as shelterbelts have proven to be effective in mitigating wind erosion. This study has the main purpose of examining the effect of shelterbelts on alleviating the wind erosion degree of stone carvings. The applicability of the canopy model for reproducing the aerodynamic effects based on the realizable k– ε and LES model was examined by using a validation metric. The shelterbelt structure has been discussed with the goal of ﬁnding the optimum canopy structure to provide a guideline for designing shelterbelts. Compared with the LES model, the realizable k– ε model was adopted in this study based on its comprehensive performance. The results show that a canopy with porosity of ϕ = 30% and a width of 0.3 to 0.5 H has better sheltering efﬁciency. Compared to the case with no shelterbelt, the wind speed ampliﬁcation coefﬁcient decreased by 43%, and the signiﬁcant decrease in the value of the wind speed ampliﬁcation coefﬁcient in the primary-harm wind direction demonstrates the effectiveness of the shelterbelt. By exploring preventive protection technology in the context of historical stone carving, this study can promote the practice of scientiﬁc and technological protection of cultural relics. around shelterbelt over a real complex terrain were simulated, and the effect of shelterbelt on improving the wind environment in the Grottoes zone and alleviating the wind erosion degree of stone carvings was examined.


Introduction
Windbreaks such as shelterbelts and fences have been used around the world for centuries to protect soils from wind erosion [1][2][3][4][5][6][7]. A shelterbelt, as a kind of natural windbreak, creates resistance to the approaching airflow and forces the wind to reduce its speed. It has proven to be effective in mitigating wind erosion. The interaction between the shelterbelt and the airflow is complicated by the turbulent characteristics of the airflow and the complex structure of the shelterbelt [8]. Therefore, the prediction of the flow fields around the shelterbelt has been conducted by numerous researchers in recent decades [9][10][11][12][13][14].
There is a requirement for wind energy resource assessment [15] and pollutants dispersion [16], where the local wind fields are significantly affected by surrounding vegetation near mountain areas. Studies of simulated flow around shelterbelts over complex terrain have generated increasing interest in recent years [17][18][19][20]. Although researchers have made great progress in the studies of flow fields around shelterbelts over flat or simple terrains, there are few studies on flow fields over real complex terrain around shelterbelts. In addition, very limited studies have recently accounted for the relationship between the local flow characteristics and wind erosion of stone carvings and how to adjust the wind environment in the Grottoes by setting up a shelterbelt.
The shelter efficiency of shelterbelts is intrinsically related to the flow characteristics of the airflow around them, and the flow fields around the shelterbelt can be divided into  [21].
Numerous experiments have been conducted to measure flow around three-dimensional vegetation in a full-scale field [1,[22][23][24]. The advantage of this method is that the full complexity of the flow characteristics is considered under real atmospheric conditions. However, this method is limited by low spatial resolution, and it is difficult to understand the shelterbelt aerodynamics mechanism. Wind tunnel experiments can control the inflow conditions, and the stationary flow conditions can be maintained throughout the entire experiment process [25][26][27][28]. However, they are also limited by spatial measurement points and require adherence to similarity criteria. With the rapid improvement of computer power and the development of computing software programs, the CFD technique has been widely used to investigate flow fields around shelterbelts [29][30][31]. CFD has some advantages compared to previous methods since it obtains flow characteristics of each position in the flow fields, and there are no problems with violating similarity requirements because simulations can be performed at full scale [32,33]. However, the accuracy and reliability of numerical simulation largely depend on the use of closure schemes and the setting of the operation of the software. For these reasons, field measurements are also imperative to verify the accuracy of CFD numerical simulations.
Until now, there have been two types of approaches to modeling vegetation's effects on flow fields. One is to simulate vegetation with corresponding roughness, affecting the wind profile above the ground surface [34], and the other is to introduce a porous medium to represent a vegetation canopy model, creating resistance as the airflow passes through by including source terms for turbulent transportation equations and a drag term in the momentum equations [11,29,35,36]. The first approach is widely used in numerical studies on the flow fields around vegetation because of its convenience. However, the main problem with this method is that it can only roughly generate the general characteristics of airflow over the vegetation, and it can only provide the average velocity profile, which is insufficient for information regarding the turbulence characterization. The second approach is to use the canopy model, which is attracting increasing attention because of its ability to physically consider the drag effect from the vegetation.
In past simulation studies on the performance of windbreak, much research on canopy models with Reynolds Averaged Navier-Stokes (RANS) model has been carried out. They concluded that the performance of the RNG (Renormalization Group) and realizable Numerous experiments have been conducted to measure flow around three-dimensional vegetation in a full-scale field [1,[22][23][24]. The advantage of this method is that the full complexity of the flow characteristics is considered under real atmospheric conditions. However, this method is limited by low spatial resolution, and it is difficult to understand the shelterbelt aerodynamics mechanism. Wind tunnel experiments can control the inflow conditions, and the stationary flow conditions can be maintained throughout the entire experiment process [25][26][27][28]. However, they are also limited by spatial measurement points and require adherence to similarity criteria. With the rapid improvement of computer power and the development of computing software programs, the CFD technique has been widely used to investigate flow fields around shelterbelts [29][30][31]. CFD has some advantages compared to previous methods since it obtains flow characteristics of each position in the flow fields, and there are no problems with violating similarity requirements because simulations can be performed at full scale [32,33]. However, the accuracy and reliability of numerical simulation largely depend on the use of closure schemes and the setting of the operation of the software. For these reasons, field measurements are also imperative to verify the accuracy of CFD numerical simulations.
Until now, there have been two types of approaches to modeling vegetation's effects on flow fields. One is to simulate vegetation with corresponding roughness, affecting the wind profile above the ground surface [34], and the other is to introduce a porous medium to represent a vegetation canopy model, creating resistance as the airflow passes through by including source terms for turbulent transportation equations and a drag term in the momentum equations [11,29,35,36]. The first approach is widely used in numerical studies on the flow fields around vegetation because of its convenience. However, the main problem with this method is that it can only roughly generate the general characteristics of airflow over the vegetation, and it can only provide the average velocity profile, which is insufficient for information regarding the turbulence characterization. The second approach is to use the canopy model, which is attracting increasing attention because of its ability to physically consider the drag effect from the vegetation.
In past simulation studies on the performance of windbreak, much research on canopy models with Reynolds Averaged Navier-Stokes (RANS) model has been carried out. They concluded that the performance of the RNG (Renormalization Group) and realizable k-e turbulence model is generally better than the standard k-ε turbulence model in terms of the prediction of mean wind speed and turbulent kinetic energy [2,37], and the realizable k-ε model is appropriate for high-Reynolds number flows. In terms of modeling a vegetation canopy, some researchers have attempted to validate the performance of the LES model [25,[38][39][40]. It is considered that the LES model has good performance in the simulation of vortices in the wake region, and the importance of turbulent inflow conditions to the simulation accuracy of the LES model is emphasized.
The Xumishan Grottoes are located in northwest China. The Grottoes were initially built in the late period of the Northern Wei Dynasty (386-534) and were listed as a key statelevel cultural site in 1982. Windy weather frequently occurs in this area, and wind erosion is a typical issue for historical stone carvings in the caves. This study has the main purpose of examining the effect of shelterbelts on alleviating the wind erosion degree of stone carvings. The applicability of the canopy model for reproducing the aerodynamic effects based on realizable k-ε and LES model was examined for a comparison between numerical results and validation metrics. Then the shelterbelt with different canopy porosities and canopy widths were discussed with the goal of finding the optimum canopy structure to provide a guideline for designing shelterbelts. At last, flow fields around the shelterbelt over a real complex terrain were simulated with the purpose of examining the effect of the shelterbelt on improving the wind environment in the Grottoes zone.

Numerical Methods and Their Applications
In this study, the governing equations are derived from the Navier-Stokes equations, either filtered (LES) or time-averaged (RANS). It is assumed that air is incompressible and exists in a neutral atmospheric boundary layer (ABL), while Coriolis force was neglected due to the typical shelterbelt height, which is much smaller than the atmospheric boundary layer. Heat and mass transfer processes are also not considered in the present study. The governing equations for the conservation of mass and momentum can be expressed as follows: where x i (1 = x = streamwise, 2 = y = spanwise, 3 = z = vertical) are the Cartesian coordinates, u i is the wind velocity in the ith direction. p is pressure, ρ is density of the fluid, µ is the kinematic viscosity and S i is a source term added to the momentum equation. The over line denotes time-averaged mean value in the simulation with the RANS model, while in the LES model, it denotes a filtered value. Meanwhile, the expression of τ ij in Equation (A2) for the RANS model and the LES model is different. In RANS are the Reynolds stresses, τ ij = −ρu i u j , while τ ij in LES indicates the subgrid-scale stress, τ ij = −ρ u i u j − u i u j and accounts for contribution from unresolved smaller vortex to a large size vortex. u i = u i − u i is the fluctuating component of the instantaneous velocity in the x i direction. Both the Reynolds stresses and the subgrid-scale stress were modeled using an Eddy-viscosity assumption. A porous medium is introduced to represent a canopy within CFD modeling, in which a drag force term was added in the momentum equations, and source terms were added in turbulence equations.
Details of this numerical method (turbulence model, modeling of vegetation canopy, boundary condition, and solution scheme) are given in Appendix A.
In this paper, five cases are discussed. Firstly, the applicability of the canopy model for reproducing the aerodynamic effects based on realizable k-ε and LES model was examined for a comparison between numerical results and validation metrics. After sufficient validations, cases with different canopy porosities and canopy widths were discussed with the goal of finding the optimum canopy structure to provide a guideline for designing shelterbelts. Finally, flow fields around shelterbelt over a real complex terrain were simulated, and the effect of shelterbelt on improving the wind environment in the Grottoes zone and alleviating the wind erosion degree of stone carvings was examined.

Validation of Numerical Method
Flow fields around the canopy model over the flat terrain were simulated by the realizable k-ε and LES model separately. The numerical results of these two models are compared with field measurement data for a row of trees in Izumo in order to validate the numerical method [41]. The length, width, and height of the canopy are 74 m, 2 m, and 7 m, respectively. The distance from the bottom of the leaf to the ground is 1.2 m. The inflow conditions were measured at 35 m upstream of the canopy at the heights of 1 m, 3 m, and 9 m. The mean wind speeds were measured at 7 m, 14 m, 21 m, 28 m, and 35 m leeward of the canopy at four different heights. The inflow wind direction was perpendicular to the canopy. The measurement was taken under neutral atmospheric conditions, which means that measurement data were not affected by local thermal effects.  amined for a comparison between numerical results and validation metrics. After sufficient validations, cases with different canopy porosities and canopy widths were discussed with the goal of finding the optimum canopy structure to provide a guideline for designing shelterbelts. Finally, flow fields around shelterbelt over a real complex terrain were simulated, and the effect of shelterbelt on improving the wind environment in the Grottoes zone and alleviating the wind erosion degree of stone carvings was examined.

Validation of Numerical Method
Flow fields around the canopy model over the flat terrain were simulated by the realizable k-ε and LES model separately. The numerical results of these two models are compared with field measurement data for a row of trees in Izumo in order to validate the numerical method [41]. The length, width, and height of the canopy are 74 m, 2 m, and 7 m, respectively. The distance from the bottom of the leaf to the ground is 1.2 m. The inflow conditions were measured at 35 m upstream of the canopy at the heights of 1 m, 3 m, and 9 m. The mean wind speeds were measured at 7 m, 14 m, 21 m, 28 m, and 35 m leeward of the canopy at four different heights. The inflow wind direction was perpendicular to the canopy. The measurement was taken under neutral atmospheric conditions, which means that measurement data were not affected by local thermal effects.  The mean wind speed profiles obtained from numerical results were compared with field measurement data, as shown in Figure 3. In general, the simulation results based on the two models are in good agreement with the measurement data, validating the accuracy of the numerical method. This may be attributed to the fact that the fluid force formula for the two models is the same, thus providing the same momentum loss in the simulation. An S-shape of the wind profiles was observed in both simulation results of the two models. This is mainly due to: (a) the increase in wind velocity at the top of the canopy; (b) the decrease in wind speed behind the canopy; and (c) insignificant acceleration at the bottom of the canopy. Additionally, the three phenomena result in strong shear stress at the top of the canopy. It can be observed from Figure 3 that, at x/H = 4 and x/H = 5, the differences between the two simulated cases are not substantial. However, in the case of x/H = 3, the numerical results with the realizable k-ε model are slightly more accurate than the LES model. The reason for this may be that the accurate prediction of the LES model depends largely on turbulent inflow conditions and grid resolution. Furthermore, their complexity requires the proper calibration of additional input parameters to produce reliable results. The mean wind speed profiles obtained from numerical results were compared with field measurement data, as shown in Figure 3. In general, the simulation results based on the two models are in good agreement with the measurement data, validating the accuracy of the numerical method. This may be attributed to the fact that the fluid force formula for the two models is the same, thus providing the same momentum loss in the simulation. An S-shape of the wind profiles was observed in both simulation results of the two models. This is mainly due to: (a) the increase in wind velocity at the top of the canopy; (b) the decrease in wind speed behind the canopy; and (c) insignificant acceleration at the bottom of the canopy. Additionally, the three phenomena result in strong shear stress at the top of the canopy. It can be observed from Figure 3 that, at x/H = 4 and x/H = 5, the differences between the two simulated cases are not substantial. However, in the case of x/H = 3, the numerical results with the realizable k-ε model are slightly more accurate than the LES model. The reason for this may be that the accurate prediction of the LES model depends largely on turbulent inflow conditions and grid resolution. Furthermore, their complexity requires the proper calibration of additional input parameters to produce reliable results.  The normalized turbulent kinetic energy profiles obtained from numerical results were compared with field measurement data, as shown in Figure 4. It can be observed that the simulated profiles by the realizable k-ε model in the wake region of the canopy  The normalized turbulent kinetic energy profiles obtained from numerical results were compared with field measurement data, as shown in Figure 4. It can be observed that the simulated profiles by the realizable k-ε model in the wake region of the canopy are slightly underestimated in the range 2 < x/H < 4. However, in research on wind flow behind windbreaks, J.L. Santiago [2] concluded that the turbulence model could be considered acceptable under the condition that turbulent kinetic energy is underestimated. In addition, this paper aims to determine the optimum canopy structure giving rise to the largest shelter effect in order to protect the stone carvings, and the wind erosion mechanism discussed in Section 3.5 shows that the most effective way to mitigate the wind erosion of the historical stones is to reduce the average wind speed under the primary-harm wind direction. Therefore, this research mainly focuses on the mean flow characteristics around the shelterbelt. The average velocity profile simulated by the two models is similar to the measured data, so the performance of the realizable k-ε model is generally acceptable. The normalized turbulent kinetic energy profiles obtained from numerical results were compared with field measurement data, as shown in Figure 4. It can be observed that the simulated profiles by the realizable k-ε model in the wake region of the canopy are slightly underestimated in the range 2 < x/H < 4. However, in research on wind flow behind windbreaks, J.L. Santiago [2] concluded that the turbulence model could be considered acceptable under the condition that turbulent kinetic energy is underestimated. In addition, this paper aims to determine the optimum canopy structure giving rise to the largest shelter effect in order to protect the stone carvings, and the wind erosion mechanism discussed in Section 3.5 shows that the most effective way to mitigate the wind erosion of the historical stones is to reduce the average wind speed under the primary-harm wind direction. Therefore, this research mainly focuses on the mean flow characteristics around the shelterbelt. The average velocity profile simulated by the two models is similar to the measured data, so the performance of the realizable k-ε model is generally acceptable. Compared with the LES model, the realizable k -ε turbulence model can better balance the computational accuracy and computational resources. Hence, the flow modeling Compared with the LES model, the realizable k-ε turbulence model can better balance the computational accuracy and computational resources. Hence, the flow modeling through a windbreak in this study is based on realizable k-ε modeling of infinite dense canopy flow.

Influence of Canopy Internal Structural Parameter
In this section, the verified canopy model was applied to discuss the influence of canopy internal structure parameters on shelter efficiency. The most used parameter of canopy internal structure is porosity (α), which is a simple ratio of perforated area to total area. For thin artificial barriers, α is equivalent to optical porosity. However, optical porosity may not always satisfactorily describe the permeability of air flow because optical porosity shows only the two-dimensional gaps and cannot represent the three-dimensional spaces through which the wind flows across windbreaks with a certain width. The aerodynamic porosity, which describes the volume porosity, is a reasonable characterization parameter when the windbreak is a shelterbelt. Consequently, we have arbitrarily varied the value of C d A(z) and calculated the resulting aerodynamic porosity α. The aerodynamic porosity is defined as [42].
where u is the axial velocity component, u 0 is the axial velocity component far upstream, and the integration is performed over vegetation canopy area S. The aerodynamic porosity is the ratio of the flow flux through the canopy cross section to the flow flux at the same cross-section far upstream. To study the effect of the aerodynamic porosity on shelter efficiency, six cases with a value of C d A(z) = 0.4, 1.6, 3.2, 6.4, 9.6 and 25.6 were simulated. The respective calculated aerodynamic porosities α = 85, 57, 45,35,30,18. Except for the vegetation canopy aerodynamic porosity introduced in the numerical simulation, all the other configuration parameters are identical to those in Section 3.1. Figure 5 presents the contours of normalized mean wind speed around the canopy with different porosities in a cross-section of y/H = 0. As shown, the approach flow separates into two parts due to the canopy: displaced flow moves upward and passes over the canopy top, and bleed flow passes through the canopy pores. The mean wind speed decreases at the leeward of the windbreak because of the blockage effect created by the canopy. The larger the aerodynamic porosity, the less obvious the blockage effect is. As the porosity increases, the displaced flow speed gradually decreases. However, the bleed flow becomes stronger, inhibiting the formation of the quiet zone. For the sparse canopy ( Figure 5a,b), the ambiguous region of the quiet zone simulated by the porous model is in agreement with an earlier study on the simulations of Poggia et al. [17] for a canopy with large porosity. Hence, the performance of the porous model describing the flow is acceptable.
where u is the axial velocity component, 0 is the axial velocity component far upstream, and the integration is performed over vegetation canopy area S. The aerodynamic porosity is the ratio of the flow flux through the canopy cross section to the flow flux at the same cross-section far upstream.
To study the effect of the aerodynamic porosity on shelter efficiency, six cases with a value of C d A(z) = 0.4, 1.6, 3.2, 6.4, 9.6 and 25.6 were simulated. The respective calculated aerodynamic porosities α = 85, 57, 45,35,30,18. Except for the vegetation canopy aerodynamic porosity introduced in the numerical simulation, all the other configuration parameters are identical to those in Section 3.1. Figure 5 presents the contours of normalized mean wind speed around the canopy with different porosities in a cross-section of y/H = 0. As shown, the approach flow separates into two parts due to the canopy: displaced flow moves upward and passes over the canopy top, and bleed flow passes through the canopy pores. The mean wind speed decreases at the leeward of the windbreak because of the blockage effect created by the canopy. The larger the aerodynamic porosity, the less obvious the blockage effect is. As the porosity increases, the displaced flow speed gradually decreases. However, the bleed flow becomes stronger, inhibiting the formation of the quiet zone. For the sparse canopy (Figure 5a,b), the ambiguous region of the quiet zone simulated by the porous model is in agreement with an earlier study on the simulations of Poggia et al. [17] for a canopy with large porosity. Hence, the performance of the porous model describing the flow is acceptable.
(a) (b) As shown in Figure 5, maximum wind reductions are closely related to porosity, with low porosity producing high maximum reductions. The velocity along z/H = 1 has a Uref = 0.8 down to 0.2 for aerodynamic porosities of 85-18% respectively. Barriers with very low porosity create increased turbulence downwind compared to the medium-dense vegetation canopy. The higher turbulence level may result in recovery of mean horizontal wind speeds to upwind speeds closer to the low-porosity vegetation canopy, thus resulting in a shorter protected distance. very low porosity create increased turbulence downwind compared to the medium-dense vegetation canopy. The higher turbulence level may result in recovery of mean horizontal wind speeds to upwind speeds closer to the low-porosity vegetation canopy, thus resulting in a shorter protected distance.
The shelter efficiency of windbreaks are of great significance in their design. Over the years, many scholars have proposed several evaluation criteria, among which the most used are the distance d 20 over which wind speed is reduced by at least 20%, and the minimal wind speed U min/ U ref on the leeward side of the windbreak and its location x min [1]. The minimum relative wind speed U min/ U ref may be of most interest where the area requiring protection is small. Since this paper mainly focuses on the average wind speed in front of the cave, U min/ U ref is used as an index for the assessment of windbreak shelter efficiency. Figure 5 also shows that the canopy with a porosity of 18-30% has the same value of U min, whilst the canopy with a porosity of 30% has a larger value of d 20 . This means that it can efficiently decrease downwind mean wind speed. Therefore, the optimal porosity of the canopy is 30%.

Influence of Canopy External Structural Parameter
Natural shelterbelts, unlike artificial fences, have a certain width. In this section, the flow fields around canopies of different widths are numerically studied, and their influences on shelter efficiency are discussed.
As described in Equations (2) and (A7), drag force per unit air mass and drag coefficient per unit LAD are local drag force and local drag coefficient of the shelterbelt, respectively. While total drag force of the entire shelterbelt is usually used in the research. For a twodimensional artificial fence, the resistance coefficient or pressure-loss coefficient k r is the most used parameter to reflect aerodynamic properties of the entire windbreak. For the vegetation canopy, the value of k r may be roughly estimated as [43]: where w s is the width of the vegetation canopy.
We can see that k r is the integration of the product of local leaf-area density and the local dynamic property of leaf area along the i-direction, which relates the bulk aerodynamic property of the entire vegetation canopy to the local LAD and local drag coefficient. Therefore, for fixed values of k r , a vegetation canopy with different widths has the same porosity. As listed in Table 1, there are six cases of vegetation canopies of spatially uniform C d A(Z) with different widths ranging from 0.1 H to 3 H. Except for the canopy with different widths introduced in the numerical simulation, all the other configuration parameters are identical to those in Section 3.1.  Figure 6 presents the contours of normalized mean wind speed around the canopy with different widths in the cross-section of y/H = 0. The width of the canopy significantly influences the perturbed flow fields on the leeside of the canopy and then affects the mean wind speed and the recovery of wind speed. For canopy width less than 1 H, the mean wind speed on the leeward side decreases as the width increases. However, for canopy width greater than 1 H, the variation tendency of mean wind speed on the leeward side is opposite to the above situation. When the width is between 0.3 H and 0.5 H, the canopy has the largest quiet zone area on the leeward side. When the width of the canopy is equal to 3 H, the flow above the canopy may completely adjust to the roughness of the canopy, so the characteristic of flow above and inside the plant canopy is heterogeneous.
with different widths in the cross-section of y/H = 0. The width of the canopy significantly influences the perturbed flow fields on the leeside of the canopy and then affects the mean wind speed and the recovery of wind speed. For canopy width less than 1 H, the mean wind speed on the leeward side decreases as the width increases. However, for canopy width greater than 1 H, the variation tendency of mean wind speed on the leeward side is opposite to the above situation. When the width is between 0.3 H and 0.5 H, the canopy has the largest quiet zone area on the leeward side. When the width of the canopy is equal to 3 H, the flow above the canopy may completely adjust to the roughness of the canopy, so the characteristic of flow above and inside the plant canopy is heterogeneous. As shown in Figure 6, the shelter efficiencies expressed in the minimal wind speed Umin/Uref on the leeward side of the windbreak and its location xmin changes little for different widths. The principal reason for this is the compensation between the permeability and the downstream perturbed pressure of the canopy. In order to save forest resources, a width of 0.3 to 0.5 H for shelterbelts with a rectangular cross-section is suggested.

Flow Fields around Shelterbelt over a Real Complex Terrain
The Xumishan Grotto Zone is located at about 105°58'46"-105°59'21" E, 36°16'13"-36°17'18" N ( Figure 7). Complex terrain is a generic term used in the literature to describe any irregular terrain. Grottoes in northwest China are often located in complex terrains such as cliffs and valleys. The Xumishan Grottoes, which are composed of valleys and hills with different elevations, are characterized by high levels of complexity. The caves of the Xumishan Grottoes are mainly distributed within a 1 km radius of Cave 5. Therefore, this paper focuses on microscale wind fields over complex terrain. As shown in Figure 6, the shelter efficiencies expressed in the minimal wind speed U min/ U ref on the leeward side of the windbreak and its location x min changes little for different widths. The principal reason for this is the compensation between the permeability and the downstream perturbed pressure of the canopy. In order to save forest resources, a width of 0.3 to 0.5 H for shelterbelts with a rectangular cross-section is suggested.

Flow Fields around Shelterbelt over a Real Complex Terrain
The Xumishan Grotto Zone is located at about 105 • 58 46 -105 • 59 21 E, 36 • 16 13 -36 • 17 18 N (Figure 7). Complex terrain is a generic term used in the literature to describe any irregular terrain. Grottoes in northwest China are often located in complex terrains such as cliffs and valleys. The Xumishan Grottoes, which are composed of valleys and hills with different elevations, are characterized by high levels of complexity. The caves of the Xumishan Grottoes are mainly distributed within a 1 km radius of Cave 5. Therefore, this paper focuses on microscale wind fields over complex terrain. The Xumishan Grotto Zone is located at about 105°58'46"-105°59'21" E, 36°16'13"-36°17'18" N ( Figure 7). Complex terrain is a generic term used in the literature to describe any irregular terrain. Grottoes in northwest China are often located in complex terrains such as cliffs and valleys. The Xumishan Grottoes, which are composed of valleys and hills with different elevations, are characterized by high levels of complexity. The caves of the Xumishan Grottoes are mainly distributed within a 1 km radius of Cave 5. Therefore, this paper focuses on microscale wind fields over complex terrain. The numerical simulation of the flow fields is divided into two levels. At the higher level, in order to reflect the influence of upstream terrain on the wind flow characteristics [34], a terrain model with a horizontal resolution of 20 m within a 10 km ×10 km region is established based on a geographic information system (GIS) data ( Figure 8). According to guidelines for CFD practices in wind engineering [44,45], the height of the computational domain is five times the maximum elevation difference in the region. Note that the height of the computational domain in this paper is similar to the atmospheric boundary layer The numerical simulation of the flow fields is divided into two levels. At the higher level, in order to reflect the influence of upstream terrain on the wind flow characteristics [34], a terrain model with a horizontal resolution of 20 m within a 10 km × 10 km region is established based on a geographic information system (GIS) data ( Figure 8). According to guidelines for CFD practices in wind engineering [44,45], the height of the computational domain is five times the maximum elevation difference in the region. Note that the height of the computational domain in this paper is similar to the atmospheric boundary layer height; the inflow turbulent kinetic energy and dissipation rate are different from those in Sections 3.1-3.3 and are formulated using Equations (5)- (8): k(z) = (U(z)·I(z)) 2 (6) ε(z) = C 3/4 u k 3/2 /(κ v z) where z ref and U ref are the standard reference height and standard reference wind speed set to 10 m and 10 m/s in this paper, respectively. The fitting result for the profile index α is 0.26 according to the drone field test wind profile. The corresponding gradient wind height is 400 m, which means wind speed no longer increases with the distance from the ground.  The Grotto Zone in this study is 1 × 1 km 2 , and a horizontal resolution of 5 m within this region was established through drone oblique photogrammetry at the current level. The wind characteristics (including the velocity profile, turbulent kinetic energy profile, and dissipation rate profile) of the higher-level simulation are used as the inlet boundary conditions of the current-level simulation [46,47] (Figure 9).
where z ref and U ref are the standard reference height and standard reference wind speed set to 10 m and 10 m/s in this paper, respectively. The fitting result for the profile index α is 0.26 according to the drone field test wind profile. The corresponding gradient wind height is 400 m, which means wind speed no longer increases with the distance from the ground. The Grotto Zone in this study is 1 × 1 km 2 , and a horizontal resolution of 5 m within this region was established through drone oblique photogrammetry at the current level. The wind characteristics (including the velocity profile, turbulent kinetic energy profile, In order to validate the simulation accuracy, the simulated flow in the complex terrain around the caves was compared with local meteorological data. The Huangduobao meteorological station (MS1) is located northeast of this region, and it was located upstream of the prevailing wind direction as a reference. The meteorological station of the Xumishan Grottoes (MS2) is located in the center of the area and can accurately monitor the wind climate around Cave 5. The location of meteorological stations can be seen in Figure 10.
The simulated values are in black, the measured values are in red, and the wind direction standard deviation of the measured values is in red. Figure 11 shows a comparison between the simulated and measured values of wind speed ratio and wind direction at Xumishan Grottoes meteorological station. The numerical values of wind speed are generally within 25% of the corresponding measurements, and the maximum difference value of the wind direction is less than 20%. The results show that the method can accurately simulate the spatial wind field of complex terrain.
Pinus is a common, tall tree planted in the study area. It is a pine evergreen coniferous tree with an average height in its forests up to 30 m. This height was used as the canopy height. Based on the previous discussion, the optimum vegetation canopy structure was confirmed to provide a guideline for designing a shelterbelt. In this simulation, the porosity of the canopy was 30%, and the width was 0.4 H. For natural barriers, as the approach angle (the angle of average wind direction from perpendicular to the windbreak) increases, the effective porosity of the windbreak becomes smaller. The shelter efficiency of a canopy is a function of approach angle and porosity. Medium-dense (ϕ = 30%) canopy reaches a maximum shelter efficiency with the approach angle value of 90 • . The distance between windbreaks and shelterbelt is about 3 H. At this distance, the wind speed is reduced to the minimum value in the approach flow. Figure 12 presents the contours of normalized mean wind speed around the shelterbelt over complex terrain in a cross-section perpendicular to the canopy. It can be found that flow characteristics downwind of the shelterbelt are obviously different in various inflow directions. When the inflow direction is northeast, the quiet zone region on the leeward side of the canopy is significantly reduced, indicating that the disturbance from the downstream hills may still affect the flow. When the inflow direction is southwest, the downwind terrain of the shelterbelt is relatively flat, and the upstream hills have little influence on flow characteristics.

level. Elevation information of Domain 2 (b).
In order to validate the simulation accuracy, the simulated flow in the complex terrain around the caves was compared with local meteorological data. The Huangduobao meteorological station (MS1) is located northeast of this region, and it was located upstream of the prevailing wind direction as a reference. The meteorological station of the Xumishan Grottoes (MS2) is located in the center of the area and can accurately monitor the wind climate around Cave 5. The location of meteorological stations can be seen in Figure 10. The simulated values are in black, the measured values are in red, and the wind direction standard deviation of the measured values is in red.  Figure 11 shows a comparison between the simulated and measured values of wind speed ratio and wind direction at Xumishan Grottoes meteorological station. The numerical values of wind speed are generally within 25% of the corresponding measurements, and the maximum difference value of the wind direction is less than 20%. The results show that the method can accurately simulate the spatial wind field of complex terrain. Pinus is a common, tall tree planted in the study area. It is a pine evergreen conifero tree with an average height in its forests up to 30 m. This height was used as the cano height. Based on the previous discussion, the optimum vegetation canopy structure w belt over complex terrain in a cross-section perpendicular to the canopy. It can be found that flow characteristics downwind of the shelterbelt are obviously different in various inflow directions. When the inflow direction is northeast, the quiet zone region on the leeward side of the canopy is significantly reduced, indicating that the disturbance from the downstream hills may still affect the flow. When the inflow direction is southwest, the downwind terrain of the shelterbelt is relatively flat, and the upstream hills have little influence on flow characteristics.    Figure 13 shows the streamline behind the shelterbelt over complex terrain. The flow fields are significantly affected by the geomorphic features of the terrain. When the inflow direction is in the northeast, due to the blocking effect of the downstream hills, the airflow accelerates between the shelterbelt and the hills. The wind speed on the leeward side of the shelterbelt increases gradually as the slope increases and reaches the maximum value at the ridge, which also explains why the quiet zone region on the leeward side of the canopy decreases significantly in Figure 13a. When the inflow direction is southwest, the flow characteristics downwind of the shelterbelt are less affected by topography. The flow separation phenomenon occurs in the area near the leeward side, and a vortex is also generated. direction is in the northeast, due to the blocking effect of the downstream hills, the airflow accelerates between the shelterbelt and the hills. The wind speed on the leeward side of the shelterbelt increases gradually as the slope increases and reaches the maximum value at the ridge, which also explains why the quiet zone region on the leeward side of the canopy decreases significantly in Figure 13a. When the inflow direction is southwest, the flow characteristics downwind of the shelterbelt are less affected by topography. The flow separation phenomenon occurs in the area near the leeward side, and a vortex is also generated.
It can be seen from the simulation results that terrain features are an important factor affecting the flow characteristics, and the downstream hills have a large influence on the flow characteristics in the quiet zone of the shelterbelt. In order to make sure that the simulated flow fields' flow characteristics are closer to the real, it is necessary to take turbulent flow, influenced by the leeside of the actual hill and the presence of the shelterbelt, into comprehensive consideration.

Relationship between Local Wind Conditions and Wind Erosion
Cave 5, which this paper focuses on, is representative of the Tang Dynasty statues in the Xumishan Grottoes. The feet of the Buddha statue in the Grottoes have been eroded by strong winds into the honeycomb caves (Figure 14a Figure 15 shows the contours of the simulated wind speed ratio in the horizontal plane at Cave 5 height under eight wind directions without a shelterbelt. When the inlet wind direction is northeast (Figure 15a), the inflow direction is close to the trend of the valley channel, and the wind speed around Cave 5 is less blocked by the mountain. An- It can be seen from the simulation results that terrain features are an important factor affecting the flow characteristics, and the downstream hills have a large influence on the flow characteristics in the quiet zone of the shelterbelt. In order to make sure that the simulated flow fields' flow characteristics are closer to the real, it is necessary to take turbulent flow, influenced by the leeside of the actual hill and the presence of the shelterbelt, into comprehensive consideration.

Relationship between Local Wind Conditions and Wind Erosion
Cave 5, which this paper focuses on, is representative of the Tang Dynasty statues in the Xumishan Grottoes. The feet of the Buddha statue in the Grottoes have been eroded by strong winds into the honeycomb caves (Figure 14a,b). Figure 15 shows the contours of the simulated wind speed ratio in the horizontal plane at Cave 5 height under eight wind directions without a shelterbelt. When the inlet wind direction is northeast (Figure 15a), the inflow direction is close to the trend of the valley channel, and the wind speed around Cave 5 is less blocked by the mountain. Another reason for this is that as the airflow from the open area to the channel terrain, it accelerates through and forms a high-speed canyon wind. For these reasons, the wind speed amplification coefficient reaches the maximum value of 2.1. When the angle between the inlet wind direction and the river channel increases, the wind speed in front of the cave decreases significantly, and the mountain barrier effect is obvious. When the wind direction is northwest (Figure 15g), the wind speed in front of the cave on the leeward side of the mountain reaches the lowest value due to the blocking effect of the mountain on the northwest side and the wind speed amplification coefficient is 0.3.
It can be seen from the simulation results that terrain features are an important factor affecting the flow characteristics, and the downstream hills have a large influence on the flow characteristics in the quiet zone of the shelterbelt. In order to make sure that the simulated flow fields' flow characteristics are closer to the real, it is necessary to take turbulent flow, influenced by the leeside of the actual hill and the presence of the shelterbelt, into comprehensive consideration.

Relationship between Local Wind Conditions and Wind Erosion
Cave 5, which this paper focuses on, is representative of the Tang Dynasty statues in the Xumishan Grottoes. The feet of the Buddha statue in the Grottoes have been eroded by strong winds into the honeycomb caves (Figure 14a,b).    accelerates through and forms a high-speed canyon wind. For these reasons, the wind speed amplification coefficient reaches the maximum value of 2.1. When the angle between the inlet wind direction and the river channel increases, the wind speed in front of the cave decreases significantly, and the mountain barrier effect is obvious. When the wind direction is northwest (Figure 15g), the wind speed in front of the cave on the leeward side of the mountain reaches the lowest value due to the blocking effect of the mountain on the northwest side and the wind speed amplification coefficient is 0.3.   Figure 16 shows the angle of attack (AOA) around Cave 5 when the inlet wind direction is northeast (Figure 16a), and the calculated values of AOA range between −5 • and 10 • , illustrating that wind characteristics are associated with lower turbulence level. If the inlet wind direction is northwest (Figure 16g), the calculated values of AOA range between −15 • and 20 • . The significant change in AOA indicates that turbulence increases locally. In general, the angle of attack is positive on the upslope and negative on the downslope. Combined with the terrain slope, it is found that the position with a larger angle of attack corresponds to a larger slope. Additionally, an angle of attack close to or above 15 • reflects flow separations, especially in downslope flow. The angle of attack changes suddenly for the northwest flow attributed to Cave 5 located at the leeward wall of the valley, resulting in a recirculation zone around Cave 5. The flow near mountain surfaces tends to travel vertically, and near-ground flow travels opposite to the inlet wind direction.

OR PEER REVIEW
15 of 24 Figure 16 shows the angle of attack (AOA) around Cave 5 when the inlet wind direction is northeast (Figure 16a), and the calculated values of AOA range between −5° and 10°, illustrating that wind characteristics are associated with lower turbulence level. If the inlet wind direction is northwest (Figure 16g), the calculated values of AOA range between −15° and 20°. The significant change in AOA indicates that turbulence increases locally. In general, the angle of attack is positive on the upslope and negative on the downslope. Combined with the terrain slope, it is found that the position with a larger angle of attack corresponds to a larger slope. Additionally, an angle of attack close to or above 15° reflects flow separations, especially in downslope flow. The angle of attack changes suddenly for the northwest flow attributed to Cave 5 located at the leeward wall of the valley, resulting in a recirculation zone around Cave 5. The flow near mountain surfaces tends to travel vertically, and near-ground flow travels opposite to the inlet wind direction.
(e) (f) The simulation results showed that Cave 5, being located in low wind speed areas, is associated with increased 3D turbulence of the northwest flow, and strong wind speed is associated with lower turbulence of the northeast flow. When the wind approaches from the northwest, the wind velocity and turbulence level in this characteristic wake region are significantly affected by the presence of the mountain. The shelter effect results in a region of low-speed flow beneath the mountain height and an increased turbulence level experienced at approximately the height of the mountain. When the wind approaches from the northeast, the effect of the mountain on the wake region is reduced, and flow with lower turbulence and high speed is observed. Figure 17a shows the contours of the hill surface pressure distribution in the absence of a shelterbelt: the flow separations lead to a negative pressure value at the ridge and the leeside pressure increases gradually. The value of pressure around Cave 5 is −95 Pa. Figure 17b shows the contours of the hill and canopy surface pressure distribution with the influence of the shelterbelt. We can see that the pressure is significantly affected by the presence of the shelterbelt. The shelterbelt behaves in a similar way to an obstacle, increasing the pressure substantially at the windward side of the shelterbelt, forming a large area of positive pressure area. During the displaced flow through the shelterbelt, a larger negative pressure area is generated on the leeward side, and a region characterized by low turbulence is generated due to the reversed flow. The value of pressure around Cave 5 is −150 Pa. The reduction of pressure on the stone-surface is beneficial to the protection of carving to a certain extent. The simulation results showed that Cave 5, being located in low wind speed areas, is associated with increased 3D turbulence of the northwest flow, and strong wind speed is associated with lower turbulence of the northeast flow. When the wind approaches from the northwest, the wind velocity and turbulence level in this characteristic wake region are significantly affected by the presence of the mountain. The shelter effect results in a region of low-speed flow beneath the mountain height and an increased turbulence level experienced at approximately the height of the mountain. When the wind approaches from the northeast, the effect of the mountain on the wake region is reduced, and flow with lower turbulence and high speed is observed. Figure 17a shows the contours of the hill surface pressure distribution in the absence of a shelterbelt: the flow separations lead to a negative pressure value at the ridge and the leeside pressure increases gradually. The value of pressure around Cave 5 is −95 Pa. Figure 17b shows the contours of the hill and canopy surface pressure distribution with the influence of the shelterbelt. We can see that the pressure is significantly affected by the presence of the shelterbelt. The shelterbelt behaves in a similar way to an obstacle, increasing the pressure substantially at the windward side of the shelterbelt, forming a large area of positive pressure area. During the displaced flow through the shelterbelt, a larger negative pressure area is generated on the leeward side, and a region characterized by low turbulence is generated due to the reversed flow. The value of pressure around Cave 5 is −150 Pa. The reduction of pressure on the stone-surface is beneficial to the protection of carving to a certain extent. The simulation results showed that Cave 5, being located in low wind speed areas, is associated with increased 3D turbulence of the northwest flow, and strong wind speed is associated with lower turbulence of the northeast flow. When the wind approaches from the northwest, the wind velocity and turbulence level in this characteristic wake region are significantly affected by the presence of the mountain. The shelter effect results in a region of low-speed flow beneath the mountain height and an increased turbulence level experienced at approximately the height of the mountain. When the wind approaches from the northeast, the effect of the mountain on the wake region is reduced, and flow with lower turbulence and high speed is observed. Figure 17a shows the contours of the hill surface pressure distribution in the absence of a shelterbelt: the flow separations lead to a negative pressure value at the ridge and the leeside pressure increases gradually. The value of pressure around Cave 5 is −95 Pa. Figure 17b shows the contours of the hill and canopy surface pressure distribution with the influence of the shelterbelt. We can see that the pressure is significantly affected by the presence of the shelterbelt. The shelterbelt behaves in a similar way to an obstacle, increasing the pressure substantially at the windward side of the shelterbelt, forming a large area of positive pressure area. During the displaced flow through the shelterbelt, a larger negative pressure area is generated on the leeward side, and a region characterized by low turbulence is generated due to the reversed flow. The value of pressure around Cave 5 is −150 Pa. The reduction of pressure on the stone-surface is beneficial to the protection of carving to a certain extent.   The deterioration of the Xumishan Grottoes is the long-term result of this, and wind blowing and rotary grinding are the main causes of wind erosion. Under the conditions of blowing and rotary grinding, the surface of the rock mass is mainly affected by pressures and wind velocity. Woodruff, N.P, first proposed a wind erosion equation that is widely used around the world [48]. The wind erosion climate factor index C, one of the five variables of the equation, is used to represent the influence of climate conditions on the wind erosion degree. The wind erosion climate factor index is given by [49]: where C is the wind erosion climatic factor index, u is the monthly average wind speed (m/s) at 2 m, P i is the monthly precipitation (mm), ETP i is the monthly potential evaporation (mm), and d is the number of days in a month. Figure 18 shows the contours of the simulated wind speed ratio with the effect of the shelterbelt. The shelterbelt generates a clear wake region behind itself and shows a quite different mean wind speed compared to the upstream region. Regarding the wind velocity, when the inlet direction is the primary-harm wind direction (Figure 18a), the wind speed amplification coefficient around Cave 5 decreases to 1.2. Compared to the case with no shelterbelt, the wind speed amplification coefficient decreased by 43%. The significant decrease in this value in the primary-harm wind direction demonstrates the effectiveness of the shelterbelt. Since the shelterbelt is modeled as a porous medium, the compensation between the effects of permeability and perturbed pressure results in the lack of recirculating eddies in the immediate lee, with decreased turbulence. Thus, Cave 5 can be located in areas with decreased 3D turbulence and low wind speed. It is, therefore, effective to use the shelterbelt to decrease the degree of wind erosion on stone carvings in this area. This conclusion can be used to guide how to adjust the wind field near the caves by setting up reasonable passive technologies such as the shelterbelt in the Grotto Zone and provides guidance for solving the problem of stone carving erosion caused by strong wind erosion in semi-open Grottoes to achieve preventive protection of stone carvings in these areas. The deterioration of the Xumishan Grottoes is the long-term result of this, and wind blowing and rotary grinding are the main causes of wind erosion. Under the conditions of blowing and rotary grinding, the surface of the rock mass is mainly affected by pressures and wind velocity. Woodruff, N.P, first proposed a wind erosion equation that is widely used around the world [48]. The wind erosion climate factor index C, one of the five variables of the equation, is used to represent the influence of climate conditions on the wind erosion degree. The wind erosion climate factor index is given by [49]: where C is the wind erosion climatic factor index, u ̅ is the monthly average wind speed (m/s) at 2 m, P i is the monthly precipitation (mm), ETP i is the monthly potential evaporation (mm), and d is the number of days in a month. Figure 18 shows the contours of the simulated wind speed ratio with the effect of the shelterbelt. The shelterbelt generates a clear wake region behind itself and shows a quite different mean wind speed compared to the upstream region. Regarding the wind velocity, when the inlet direction is the primary-harm wind direction (Figure 18a), the wind speed amplification coefficient around Cave 5 decreases to 1.2. Compared to the case with no shelterbelt, the wind speed amplification coefficient decreased by 43%. The significant decrease in this value in the primary-harm wind direction demonstrates the effectiveness of the shelterbelt. Since the shelterbelt is modeled as a porous medium, the compensation between the effects of permeability and perturbed pressure results in the lack of recirculating eddies in the immediate lee, with decreased turbulence. Thus, Cave 5 can be located in areas with decreased 3D turbulence and low wind speed. It is, therefore, effective to use the shelterbelt to decrease the degree of wind erosion on stone carvings in this area. This conclusion can be used to guide how to adjust the wind field near the caves by setting up reasonable passive technologies such as the shelterbelt in the Grotto Zone and provides guidance for solving the problem of stone carving erosion caused by strong wind erosion in semi-open Grottoes to achieve preventive protection of stone carvings in these areas.

Conclusions
Based on the RANS equation and realizable k-ε model, flow fields over complex terrain are investigated in this paper. The main conclusions are summarized as follows: The canopy model combined with the realizable k-ε turbulence model can accurately simulate the flow fields around vegetation. The results have shown that the canopy model has a good performance on parameterizations of aerodynamic effects of trees via CFD simulations.
The interaction of penetrating flow with the perturbation pressure and displace flows over the canopy creates a maximum wind speed reduction downstream of the canopy. The results show that a canopy with porosity of φ = 30% and a width of 0.3 H to 0.5 H has a better shelter efficiency. Therefore, the resulting flow fields around the canopy also provided quantitative information to help in the assessment of designing shelterbelt and the effect of vegetative windbreak on alleviating the wind erosion degree of stone carvings.
Terrain features are an important factor affecting flow characteristics, and the downstream hills have a large influence on the flow characteristics in the quiet zone of the shelterbelt. In order to make sure that the simulated flow fields are similar to the real flow characteristics, it is necessary to take turbulent flow that is influenced by the leeside of the actual hill and the presence of the shelterbelt into comprehensive consideration.
The significant decrease in this value in the primary-harm wind direction demonstrates the effectiveness of the shelterbelt. Since the shelterbelt is modeled as a porous medium, the compensation between the effects of permeability and perturbed pressure results in the lack of recirculating eddies in the immediate lee, with decreased turbulence. Thus, Cave 5 can be located in areas with decreased 3D turbulence and low wind speed. It is, therefore, effective to use the shelterbelt to decrease the degree of wind erosion on stone carvings in this area. This conclusion can be used to guide how to adjust the wind field near the caves in the Grotto Zone and provides guidance for solving the problem of

Conclusions
Based on the RANS equation and realizable k-ε model, flow fields over complex terrain are investigated in this paper. The main conclusions are summarized as follows: The canopy model combined with the realizable k-ε turbulence model can accurately simulate the flow fields around vegetation. The results have shown that the canopy model has a good performance on parameterizations of aerodynamic effects of trees via CFD simulations.
The interaction of penetrating flow with the perturbation pressure and displace flows over the canopy creates a maximum wind speed reduction downstream of the canopy. The results show that a canopy with porosity of ϕ = 30% and a width of 0.3 H to 0.5 H has a better shelter efficiency. Therefore, the resulting flow fields around the canopy also provided quantitative information to help in the assessment of designing shelterbelt and the effect of vegetative windbreak on alleviating the wind erosion degree of stone carvings.
Terrain features are an important factor affecting flow characteristics, and the downstream hills have a large influence on the flow characteristics in the quiet zone of the shelterbelt. In order to make sure that the simulated flow fields are similar to the real flow characteristics, it is necessary to take turbulent flow that is influenced by the leeside of the actual hill and the presence of the shelterbelt into comprehensive consideration.
The significant decrease in this value in the primary-harm wind direction demonstrates the effectiveness of the shelterbelt. Since the shelterbelt is modeled as a porous medium, the compensation between the effects of permeability and perturbed pressure results in the lack of recirculating eddies in the immediate lee, with decreased turbulence. Thus, Cave 5 can be located in areas with decreased 3D turbulence and low wind speed. It is, therefore, effective to use the shelterbelt to decrease the degree of wind erosion on stone carvings in this area. This conclusion can be used to guide how to adjust the wind field near the caves in the Grotto Zone and provides guidance for solving the problem of stone carving erosion caused by strong wind in order to achieve preventive protection of stone carvings in the Grottoes. In the present simulation, the turbulent flow is simulated using the realizable k-ε turbulence model. The turbulence closure equations for the realizable k-ε turbulence model are given as follows: where G k is the turbulent kinetic energy production, and S k and S ε are source terms for turbulent kinetic energy and its dissipation rate. σ k and σ ε are turbulent Prandtl number for k and ε, respectively. C 1 = max 0.43, η η+5 , η = S k ε , S = 2S ij S ij , and the turbulence model constants are C 2 = 1.9, σ k = 1.0, σ ε = 1.2 [37].
The source terms for the turbulent kinetic energy and the dissipation rate can be expressed as [50]: where β p , β d , C ε4 and C ε5 are model constants, and the optimal values of these parameters depend on the study case. Based on the previous literature, β p = 1, β d = 4 and C ε4 = C ε5 = 1.5 were used in this study [51].

Appendix B
According to guidelines for CFD practices in wind engineering [44], the size of the computational domain for both simulations with the realizable k-ε and LES model is shown in Figure A1, in which H is the height of vegetation canopy. This shows that the profiles of the mean wind speed and turbulent kinetic energy would be valid in the current simulation box. Boundary layer elements (prismatic cells) parallel to the walls can better meet the flow of wind in the boundary layer. In the vertical plane, the height of the first layer near the ground is set as 0.05 m with an expansion factor of 1.2 in accordance with the standard wall function requirements. Tetrahedral unstructured grids are filled to adapt the geometry. The geometry was refined in two zones. Zone 1 coincided with the volume of the porous subdomain representing the vegetation canopy. Zone 2 coincided with the remaining volume of the domain representing the flow fields. A uniform grid size was used in Zone 1, and a non-uniform grid size was adopted in Zone 2. To ensure the accuracy of the simulation results and to reduce the use of computing resources, three grid meshing schemes are adopted for the grid sensitivity analysis. Details of the maximum element size applied to these zones for each mesh can be found in Table A1. Similar profiles are obtained by the medium and fine mesh scheme, indicating that medium grid is fine enough.   Simulations with the realizable k-ε model and the LES model have the same boundary conditions at outlet, top, and side of the computational domain. The sides and top surfaces of the computational domain are symmetry boundary conditions. At the outlet, a zero-pressure outlet boundary is set. A wall function with an aerodynamic roughness of 0.01 m is used on the ground. The aerodynamic roughness length zo is converted into the corresponding wall function parameters and Cs [52]. The canopy surface is set as the interior to implement the transition between the"porous scale" and the microscale. As the calculation domain height in Sections 3.1 to 3.3 is much smaller than the atmospheric boundary layer, a logarithmic velocity profile was specified to fit the field measured data as the inflow velocity profile with the realizable k-ε model. The inflow profiles of mean wind speed U(z), turbulent kinetic energy k, turbulence dissipation rate ε and aerodynamic roughness parameters are defined by Equations (A12)-(A15), respectively [53]: Simulations with the realizable k-ε model and the LES model have the same boundary conditions at outlet, top, and side of the computational domain. The sides and top surfaces of the computational domain are symmetry boundary conditions. At the outlet, a zero-pressure outlet boundary is set. A wall function with an aerodynamic roughness of 0.01 m is used on the ground. The aerodynamic roughness length z o is converted into the corresponding wall function parameters ks and Cs [52]. The canopy surface is set as the interior to implement the transition between the "porous scale" and the microscale. As the calculation domain height in Sections 3.1-3.3 is much smaller than the atmospheric boundary layer, a logarithmic velocity profile was specified to fit the field measured data as the inflow velocity profile with the realizable k-ε model. The inflow profiles of mean wind speed U(z), turbulent kinetic energy k, turbulence dissipation rate ε and aerodynamic roughness parameters are defined by Equations (A12)-(A15), respectively [53]: where z is the height from the ground, κ is the von Karman constant with a value of 0.41. u_* is friction velocity with a value of 0.2, which was deduced from inflow conditions measured at 35 m upstream of trees with heights of 1 m, 3 m, and 9 m. In the simulations with the LES model, a random flow generation (RFG) technique is used for the inlet boundary. Firstly, an inlet boundary the same as in the realizable k-ε model was used to simulate the steady flow. Then, the stationary flow field was used as the initial flow field by LES to simulate unsteady flow. A non-dimensional time step ∆tUref/H = 0.005 is used in this study. When the time period reached 600, the relative error at the reference is less than 3%, which means that the flow reached a stable stage. The convergence criterion is that the residuals stay at a relatively low value (continuity: 10 −3 , k and ε: 10 −4 , and x-y-and z-momentum: 10 −5 ), and the wind speed of key points basically does not change. The commercial software package FLUENT (Ansys Inc.), based on finite volume method was used to solve the governing equations. The semi-implicit method for the pressure-linked equation (SIMPLE) algorithm is used to solve the pressure velocity coupling. The pressure term is discretized in the second-order scheme, and the convective term and viscous term are discretized in the second-order quadratic upstream interpolation for convective kinematics (QUICK) scheme.
The canopy parameters used in CFD modeling have the same values in the literature [54], with the equivalent leaf drag coefficient C d = 1.6 and A(z) = 1.17 m −1 in the drag force term for vegetation. It is noteworthy that the value of the drag coefficient is set to reflect an average value instead of a value constant with height in this research. The simulated and measured profiles of U(z) and k at position x/H = −5 are shown in Figure A2. boundary. Firstly, an inlet boundary the same as in the realizable k-ε model was used to simulate the steady flow. Then, the stationary flow field was used as the initial flow field by LES to simulate unsteady flow. A non-dimensional time step ΔtUref/H = 0.005 is used in this study. When the time period reached 600, the relative error at the reference is less than 3%, which means that the flow reached a stable stage.
The convergence criterion is that the residuals stay at a relatively low value (continuity: 10 −3 , k and ε: 10 −4 , and x-y-and z-momentum: 10 −5 ), and the wind speed of key points basically does not change. The commercial software package FLUENT (Ansys Inc.), based on finite volume method was used to solve the governing equations. The semi-implicit method for the pressure-linked equation (SIMPLE) algorithm is used to solve the pressure velocity coupling. The pressure term is discretized in the second-order scheme, and the convective term and viscous term are discretized in the second-order quadratic upstream interpolation for convective kinematics (QUICK) scheme.
The canopy parameters used in CFD modeling have the same values in the literature [54], with the equivalent leaf drag coefficient = 1.6 and A(z) = 1.17 m −1 in the drag force term for vegetation. It is noteworthy that the value of the drag coefficient is set to reflect an average value instead of a value constant with height in this research. The simulated and measured profiles of U(z) and k at position x/H = −5 are shown in Figure A2.