The field experiment was conducted in an oasis of Minqin County, Gansu Province, China (38° 34′ 28″N, 102° 59′ 05″E; Fig. 1). The area is situated in the lower reaches of the Shiyang River Basin. Mean annual precipitation is 116.36 mm and mean annual evaporation is 2383.7 mm. The prevailing wind comes from northwest and the mean occurrence of gale-force events (velocity≥ 17 m/s) is 27.4 d/a with an annual mean wind velocity of 2.4 m/s and an annual mean sandstorm occurrence of 25 d/a.

Fig. 1

	Figure Option ViewDownloadNew Window
	Fig. 1 Location of the study area in Minqin County, Gansu Province, China

We monitored the wind velocity of gale events from February to April 2006 at the heights of 1-49 m from two 50-m high observation towers constructed by the Gansu Desert Control Research Institute. One of the two observation towers is located at the edge of the Minqin Oasis and the other inside the oasis. The tower at the edge of the oasis is situated directly at the upstream of the tower inside the oasis with reference to the prevailing wind direction (northwest). The distance between the two towers is 3.5 km and multiple arrays of sand-fixing and windbreak forests are distributed between these two towers. The sand-fixing forests mainly consist of Haloxylon ammodendron, with a mean height of 1.6 m, a mean crown dimension of 0.89 m, a mean density of 570 trees/hm² and a mean array spacing of 9 m. The windbreak forests mainly consist ofPoplar sp., with a mean height of 18 m, a mean crown dimension of 3.5 m, varying density and varying array spacing. Wind velocity was observed at 13 different heights (1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45 and 49 m) by three-cup anemometers (DES1, Changchun, China) at each of these two towers. When the prevailing wind from the northwest reached a certain velocity (e.g., ~5 m/s at the height of 1 m), the mean wind velocity at different heights was recorded at an interval of 5 min by a data logger. All the data were analyzed using OriginPro 7.5 (OriginLab Corp, Northampton, MA, USA).

In this study, the FLUENT software was used, along with the gird-generating GAMBIT 2.2 software for pre-modeling configuration.

2.3.1 Basic assumptions

The following assumptions were made: incompressible flow, steady state, flat terrain, and homogeneous sand-fixing and windbreak forests.

2.3.2 Model selection and governing equations

There are many turbulence models to describe airflow patterns in the atmospheric boundary layer. Among these, the k-ε model has a proven ability to simulate airflow patterns above a forest (Kobayashi et al., 1994; Liu et al., 1996; Sanz, 2003; Katul et al., 2004; Alinot and Masson, 2005, Li et al., 2006; Blocken et al., 2007; Dalpé and Masson, 2009). FLUENT software, the adopted software for this study, provides three kinds of k-ε models: standard, renormalization-group (RNG), and realizable k-ε models. Santiago et al. (2007) assessed the performance of these three turbulence models using wind tunnel experiments and found that the RNG and realizable k-ε models seemed to perform better than the standard k-ε model. In this study, we selected the RNG k-ε model to carry out the numerical simulations.

The simulations were based on Reynolds-averaged Navier-Stokes (RANS) equations. Four types of equations were solved: continuity Equation 1, RANS Equation 2, and two turbulence closure equations for turbulent kinetic energy (k), Equation 3 and the dissipation rate of turbulent kinetic energy (ε ), Equation 4. Besides, the Reynolds stress was calculated by Equation 5.

Where, u_j is the jcomponent of velocity (m/s), x_j thej coordinate; u_i (m/s) is the icomponent of velocity, t (s) the time, x_i thei coordinate, ρ (kg/m³) the density of the air, μ (m²/s) the dynamic viscosity, the Reynolds stress, g_i (m/s²) the gravitational body force; μ _eff is the effective turbulent viscosity, G_k represents the turbulence kinetic energy resulted from the mean velocity gradients, α _k and α _ε are the inverse effective Prandtl number for kand ε , respectively, η is the ratio of Skto ε (note: S is the scalar measure of deformation tensor), δ _ij the Kronecker symbol. The default values of the model constants η ₀, β , C_{ε 1}, C_{ε 2}, and C_μ are listed in Table 1 (Shih et al., 1995; Fluent, 1996; Chan et al., 2002; Santiago et al., 2007).

Table 1

Table 1

Table 1 Default values for the constants of the k-ε turbulence model Parameter η 0 β Cε 1 Cε 2 Cμ Value 4.38 0.012 1.42 1.68 0.0845

Table 1 Default values for the constants of the k-ε turbulence model

2.3.3 Modeling sand-fixing and windbreak forests and their interaction with wind

The CFD model is two-dimensional one: 4300-m long and 200-m high. It consists of multiple arrays of sand-fixing forests (a mean height of 1.6 m and an equally-distanced array spacing of 9 m) and 12 arrays of windbreak forests (a mean height of 18 m and a varying array spacing of 72-410 m), being similar to the field reality (Fig. 2).

Fig. 2

Figure Option ViewDownloadNew Window

Fig. 2 Sketch map of the computational domain and the study locations (A-M). Note 1: the horizontal scale was greatly compressed in comparison with the vertical scale. Note 2: the distances both between the locations A and B and between the locations B and C within the sand-fixing forest were 500 m. Note 3: the distances between two neighboring study locations within the windbreak forests are varying, ranging from 72 to 410 m.

2.3.4 Domain and grid

In this study, a two-dimensional grid was generated using the grid generator GAMBIT. Taking into account of the computer capacity and also the expectation for the result accuracy, we set the distance between the horizontal and vertical grid points around the sand-fixing and windbreak forests at less than 0.4 m. The grid space increased with height up to 3 m (the top of the calculation domain). A total of 13 study locations (A-M) were chosen in this study (Fig. 2).

2.3.5 Boundary conditions

At the inlet boundary, the inflow velocity profile along a direction perpendicular to the flat surface obtained from the tower at the edge of the oasis followed a logarithmic equation (Eq. 6)

defined by the UDF function in FLUENT software. The turbulent kinetic energy of velocity inlet boundary (k_in) and the dissipation rate of turbulent kinetic energy of velocity inlet boundary (ε _in) were obtained by Equations 7 and 8, respectively (Fluent, 1996).

Where, u_(z) (m/s) is the mean wind velocity at height z (m), u_* (m/s) the friction velocity, and z₀(m) the aerodynamic roughness length, u_* and z₀were calculated from the velocity profile obtained from the tower at the edge of the oasis, Kis the Karman's constant, being equal to 0.4; u_avg (m/s) is the mean flow velocity, I (%) the turbulence intensity that was calculated from the field experiment data; C_μ is an empirical constant specified in the turbulence model (approximately 0.09), and l (m) is the turbulence length scale.

At the outlet boundary, the flow was assumed to be a fully developed flow. Because the top boundary was sufficiently high above the ground surface, symmetrical boundary conditions were applied. At the ground boundary, a no-slip condition was used. Each array of sand-fixing forests or windbreak forests was simplified as an idealized porous media. The physical effect of the presence of a porous medium on the flow was equivalent to a pressure drop through it, creating a momentum sink. The pressure change (Δ p) was defined as a combination of Darcy’ s Law and an additional inertial loss term, which could be expressed as Equation 9 (Fluent, 1996).

Where, α (m²) is the permeability of the medium, C₂ (1/m) the pressure-jump coefficient, v(m/s) the velocity normal to the porous face, and Δ m (m) the thickness of the medium.

It general, the mean wind velocity in a neutral turbulent boundary layer with uniform air density and with a small temperature gradient follows the logarithmic equation (Eq. 6). However, over large rough surfaces such as trees or shrubs, it should be calculated by Equation 10.

Where, u_(z) (m/s) is the mean wind velocity at height z (m), and d (m) the zero-displacement height, being generally dependent on the height of trees or shrubs.

At the edge of the oasis, the wind velocity of gale events at the heights of 1-49 m exhibited a logarithmically linear relationship with the height (Fig. 3a). We calculated u_* and z₀ based on the observed wind velocities. The mean value ofu_* was 0.80 m/s and the mean value of z₀ was 13.7 mm. However, the wind velocity profiles of gale events at the heights of 1-49 m inside the oasis differed significantly from those at the edge of the oasis (Fig. 3b). To study the effects of sand-fixing and windbreak forests on wind flow characteristics, we divided the wind velocity profiles inside the oasis at the heights of 1-49 m into two portions based on the mean height of windbreak forests: upper portion above 18 m and lower portion below 18 m. The upper portion profiles well fit the curves defined by Equation 10, in which d was 9 m (half of the mean height of windbreak forests), mean u_* was 1.16 m/s and mean z₀ was 444 mm. Our comparison shows that both u_* and z₀ were higher in the upper portion inside the oasis than the counterparts at the edge of the oasis, manifesting the important influences of sand-fixing and windbreak forests on airflow in the upper portion. The lower portion profiles do not fit the curves defined either by Equation 6 or by Equation 10. Figure 3b shows that the wind velocity increases faster with increasing height in the lower portion than in the upper portion inside the oasis.

Fig. 3

	Figure Option ViewDownloadNew Window
	Fig. 3 Observed wind velocity of seven gale events in 2006 at the heights of 1-49 m at the edge of the oasis (a) and inside the oasis (b)

To validate the simulation results, we used the wind velocity profiles observed at the edge of the oasis during gale events as the inflow velocity profiles and we then simulated the wind velocity profiles for both inlet (i.e., the tower at the edge of the oasis) and outlet (i.e., the tower inside the oasis) using the CFD model. Our comparison between the observed data and the simulated results (Fig. 4) shows that the simulated results well match the observed data, boosting our confidence on using the CFD model to simulate the wind velocity profiles.

The simulated results of wind velocity profiles at different locations (i.e., locations A, B, C, ..., M in Figure 2) between the two observation towers are shown in Fig. 5. With increasing number of sand-fixing forest and windbreak forest arrays, the velocity generally decreased along the prevailing wind direction and the wind velocity profiles within the considered heights (i.e., 1-49 m) varied significantly. Figure 5 also shows that the lower-standing sand-fixing forests (i.e., at locations A, B and C in Figure 2) had considerably less effects than the higher-standing windbreak forests (i.e., at locations from D to M in Figure 2). Specifically, over the 1000-m horizontal distance from locations A to C where there existed many arrays of sand-fixing forests, the wind velocity seemed to be influenced by the sand-fixing forests only within the heights of 1-30 m. However, once the airflow passed through the third array of windbreak forests, the wind velocity was greatly reduced from the ground level (i.e., 1 m) to the top of the considered layer (i.e., 49 m) and the velocity reduction was most significant within the heights of 1-40 m.

In addition, just as the aerodynamic roughness length (z₀) is dependent on several characteristics of the surface features (e.g., height, shape, spacing, etc.) (Wolf and Nickling, 1993), the wind velocity is also dependent on several characteristics of the surface features and the surface features here refer to sand-fixing forests and windbreak forests. It means that the wind velocity does not necessarily always decrease with increasing arrays of sand-fixing and windbreak forests. For example, although the study locations E, F, G and H were situated at the upstream of the study location J along the prevailing wind, the wind velocities at these locations (i.e., E-H) within 1-20 m heights were lower than the velocity at location J because of the shorter array spacing (Fig. 6). It should also be noted that the aerodynamic roughness lengths (z₀) of these locations (i.e., E-H) above 20 m, calculated with Equation 6, was averagely 3381 mm, being much higher than that of location J (1560 mm).

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 Observed and simulated wind velocities at the heights of 1-49 m inside the oasis and at the edge of the oasis. The observed wind velocity data were obtained on 3 March 2006 from the two 50-m high observation towers: one inside the oasis and another at the edge of the oasis (tower locations in Fig. 1).

Fig. 5

	Figure Option ViewDownloadNew Window
	Fig. 5 Simulated wind velocity profiles at locations A-M at the heights of 1-49 m (see Fig. 2 for A-M locations)

Fig. 6

	Figure Option ViewDownloadNew Window
	Fig. 6 Simulated wind velocity profiles at locations E, F, G and J below 20-m height (see Fig. 2 for E, F, G and J locations)

To detect the influences of the surface features on wind velocity along a horizontal direction, we plotted the horizontal variations of wind velocity at two heights: 9.0 and 22.5 m. Figure 7 shows that at both the heights, the wind velocity decreased generally with increasing number of sand-fixing forest and windbreak forest arrays with dramatic variations around the windbreak forests, being quite similar to the experimental results by Bradley and Mulhearn (1983) and also by Santiago et al. (2007). It is rather noticeable that those short-spaced windbreak forest arrays (i.e., locations E, F, G and H) are extremely efficient in reducing the wind velocity at the height of 9.0 m. At the height of 22.5 m, those short-spaced windbreak forest arrays (i.e., locations E, F, G and H) are still effective in reducing the wind velocity, but the efficiency is considerably lowered in comparison with that at the height of 9.0 m.

Fig. 7

	Figure Option ViewDownloadNew Window
	Fig. 7 Horizontal variations of wind velocity at the heights of 9.0 and 22.5 m from locations B to K (see Fig. 2 for E, F, G and H locations)

It is notable that the horizontal variations in wind velocity at the height of 9.0 m are consistent with the experimental results of Nelmes et al. (2001) and Torita and Satou (2007) in showing that the wind velocity decelerates as airflow approaches to the windbreak forests and that the wind velocity accelerates as airflow passes over the windbreak forests. Moreover, our results in the case of 22.5 m height are consistent with the experimental results of Bradley and Mulhearn (1983), i.e., the wind velocity accelerates as airflow approaches to the windbreak forests and then decelerates as airflow passes over the windbreak forests.

To further identify the combined influences of sand-fixing and windbreak forests on wind velocity, we eliminated all sand-fixing forests behind the first array of windbreak forest (location D) in the CFD model and calculated the new wind velocity profiles. Figure 8 shows the wind velocity profiles with the sand-fixing forests and without the sand-fixing forests accompanying the windbreak forests, and Figure 9 shows the discrepancy of wind velocity at the height of 1.0 m from locations D to L.

Fig. 8

	Figure Option ViewDownloadNew Window
	Fig. 8 Simulated wind velocity profiles with and without the sand-fixing forests within the windbreak forests from locations D to L at the heights of 1-49 m (see Fig. 2 for locations D-L)

Fig. 9

	Figure Option ViewDownloadNew Window
	Fig. 9 Wind velocity discrepancy between with and without the sand-fixing forests within the windbreak forests at the height of 1.0 m from location D to location L (see Fig. 2 for locations of D-L)

It is obvious that the sand-fixing forests existed within the windbreak forests plays a detectable role in reducing wind velocity and also in reinforcing the sheltering effect of shelterbelt system. However, the efficiency in reducing wind velocity and thus the efficiency in wind-protection effects varied with the array spacings of windbreak forests where the sand-fixing forests existed. For example, from locations D to G, the discrepancy in wind velocity at the height of 1.0 m decreased with increasing array spacings of windbreak forests. And, from locations G to L, the discrepancy in wind velocity at the height of 1.0 m increased with increasing array spacings of windbreak forests (see Fig. 9). These results demonstrate that the wind protection effects of the sand-fixing forests within the entire forest protection system are the most efficient when the arrays of the windbreak forests are adequately spaced (i.e., from location I to location L). However, if the array spacing of windbreak forests is smaller than seven times of the height of windbreak forests (i.e., from location E to location H), the effects of sand-fixing forests in reducing wind velocity will be diminished or completely masked by the effects of windbreak forests.

Our simulated results confirmed the finding of Judd et al. (1996) and the finding was that the array spacing is a controlling physical parameter for the wind flows in a windbreak system. It means that the numerical simulation base on the experimental parameters is an effective method for studying the interactions of wind flow with array spacing of windbreak forests and thus a useful tool to develop strategies in arranging the patterns of sand-fixing and windbreak forests for preventing wind erosion in the most convenient and the cheapest ways.