Optimal array of sand fences

TL;DR

This study uses Computational Fluid Dynamics to determine the optimal design of sand fence arrays for soil protection, balancing material cost and effectiveness in reducing wind shear velocity.

Contribution

It introduces a quantitative CFD-based approach to optimize sand fence arrays considering porosity, spacing, and height, advancing beyond empirical methods.

Findings

Optimal fence height is around 50 cm for minimal material use.

Fences of about 1.25 m height are most costly and less efficient.

Flow profiles and wind shear reduction depend on fence parameters.

Abstract

Sand fences are widely applied to prevent soil erosion by wind in areas affected by desertification. Sand fences also provide a way to reduce the emission rate of dust particles, which is triggered mainly by the impacts of wind-blown sand grains onto the soil and affects the Earth's climate. Many different types of fence have been designed and their effects on the sediment transport dynamics studied since many years. However, the search for the optimal array of fences has remained largely an empirical task. In order to achieve maximal soil protection using the minimal amount of fence material, a quantitative understanding of the flow profile over the relief encompassing the area to be protected including all employed fences is required. Here we use Computational Fluid Dynamics to calculate the average turbulent airflow through an array of fences as a function of the porosity, spacing and height of the fences. Specifically, we investigate the factors controlling the fraction of soil area over which the basal average wind shear velocity drops below the threshold for sand transport when the fences are applied. We introduce a cost function, given by the amount of material necessary to construct the fences. We find that, for typical sand-moving wind velocities, the optimal fence height (which minimizes this cost function) is around $50\,$cm, while using fences of height around $1.25\,$m leads to maximal cost.

Equations

$$u_0(z)=\frac{u_{\ast 0}}{\kappa }\log \frac{z-h}{\delta },$$

$${\tau }_0={\rho }_{\mathrm{air}}{u}_{\ast 0}^2,$$

$$\Delta p(z)=-\frac{1}{4{\Phi }^2}{\rho }_{\mathrm{air}}[u(z){]}^2\Delta m,$$

$$\mathcal{C}=\frac{S}{L_t}{h}_{\mathrm{f}},$$

Full text

Introduction

The transport of sand by wind and the concatenated erosion of sediment soils is one of the main causes for the propagation of desertification. Aeolian transport of sand particles is mainly due to saltation, i.e. particles move on approximately ballistic trajectories thereby ejecting new grains upon collision with the soil (splash) [1]. Moreover, the impacts of sand grains on the soil during saltation are a main factor for the emission of atmospheric dust particles [2] — which, once entrained, may travel thousands of kilometers in suspension thereby substantially affecting the Earth’s climate [3]. To prevent sand transport by wind is thus a concern of broad implication for the society .

Sand fences of various types have been constructed for centuries to control wind erosion and induce dune formation (Fig. 1). Typically, sand fences consist of lightweight wood strips, wire or perforated plastic sheets attached to regularly spaced stakes [4]. Indeed, the major pre-requisite for a sand fence is that its structure reduces the wind speed, but does not completely block the wind. Indeed, porous fences produce a longer area of leeward sheltered ground than solid fences do — the latter may also induce strong vortices that extend up to several barrier heights downwind [5]. Notwithstanding the large range of fence designs, all fences operate on the principle to create areas of low wind velocity both in front and behind the fence. In order to protect sand soils from wind erosion, often an array of sand fences is applied, where the fences are erected sequencially at a given spacing along the wind direction [4]. Depending on the area to be protected a large amount of material may be required to construct the fences. Moreover, the fences must be regularly maintained and replaced due to abrasion of fence’s material caused by wind-blown sand. In the present work, we address the problem of predicting the optimal array of fences, that is the array that uses the minimal amount of material necessary to protect a given area of sediment soil from wind erosion.

Many wind tunnel studies [6, 7, 8, 9, 10, 11, 12, 13, 14, 15], field works [16, 17] and numerical simulations [13, 18, 19, 20, 21, 22, 23, 24] have been performed in order to investigate the characteristics of the turbulent wind flow or sand flux around different types of fences. These studies showed that the amount of sand trapped depends on the fence height, its porosity, the number of fences, their spacing and the wind velocity (for a review see e.g. Ref. [25]). For sand fences, a porosity of $40\%$ or $50\%$ is recommended since it leads to optimal shielding while avoiding the formation of strong vortices. However, none of these studies focused on adjusting the design of the fence array to reduce building cost. Therefore, we investigate the shear velocity over an array of fences by means of Computational Fluid Dynamic modeling (described below) using the aforementioned porosity values as well as a representative sand-moving wind velocity (defined below). Moreover, we introduce a cost function (presented later in this manuscript), which depends on the fence height and spacing — to quantify the amount of material needed to construct an array of fences to protect the total area of soil. We will show how this function can be used to obtain the optimal height of an array of sand fences that minimizes the amount of material employed in the fences.

The schematic representation of the setup employed in our calculations is shown in Fig. 2.

The fences are placed on top of the bottom wall of a two-dimensional channel of height ${\Delta}z=10\,h_{\mathrm{f}}$ and width ${\Delta}x=80h_{\mathrm{f}}+10L$. Moreover, the soil level in the absence of the fences is considered constant and equal to zero (to simulate an approximately flat sand bed). To avoid that the results are affected by border effects, the dimensions of the box are such that the fences are far enough from the top and side walls of the channel. In other words, we have checked that the results of our calculations do not change significantly if the size of the box is increased. As depicted in Fig. 2, the wind velocity $u_{0}(z)$ at the inlet increases with the logarithm of the height $z$ above the bed level ($h$) [1, 4], In particular, $u_{0}(z)=0$ for $z-h=\delta$, where $\delta$ is the surface roughness, and increases with the height above the ground according to the following equation (valid for $z-h\geq\delta$)[1, 4]:

[TABLE]

where $\delta$ is the surface roughness, $\kappa=0.4$ the von Kármán constant and $u_{{\ast}0}$ the upwind shear velocity of the wind. The shear velocity $u_{{\ast}0}$, which gives the mean (upwind) flow velocity gradient with the height above the soil, is used to define the (upwind) shear stress,

[TABLE]

where ${\rho}_{\mathrm{air}}=1.225\,$kg$/$m3 stands for air density and $\delta=100\,{\mu}$m. We note that this value of $\delta$ has been obtained in Ref. [26] by fitting Eq. (1) to the steady-state wind profile within the numerical wind tunnel (Fig. 2). In Ref. [26] this wind profile has been generated by imposing a pressure difference between the in- and outlet of the simulation box, inducing different flow speeds [26]. Here, we have observed that using other values of surface roughness $\delta$ within the range between $10\,{\mu}$m and $1.0\,$mm [4] does not change much the shear velocity values obtained in our computations. The boundary conditions, discretization scheme and turbulence model are discussed in detail in Section Methods.

In the CFD simulation, each sand fence is modeled as a vertical, porous wall of height $h_{\mathrm{f}}$, which is varied in the present work from $10\,$cm to $2\,$m. Moreover, each fence consists of a special type of boundary condition which mimics a porous membrane of a certain velocity/pressure drop characteristics [27, 22, 28, 29]. Specifically, the pressure drop at height $z$ is given by the equation

[TABLE]

where $u(z)$ is the wind velocity normal to the fence, that is the horizontal wind speed at height $z$, ${\Delta}m$ is the fence’s thickness and $\Phi$ its porosity. In the present work, the fences’ thickness is set as ${\Delta}m=10^{-4}\,$m, while the effect of different values of porosity is investigated.

Results

Since most wind-tunnel and field experiments aimed to gain understanding on the effect of fences on sand transport have been performed using one fence, we have started our investigation using one single fence as well. Figure 3 shows results from calculations performed using a fence of height $h_{\mathrm{f}}=20\,$cm and different values of porosity. The upwind shear velocity in Fig. 3 is $u_{{\ast}0}=0.4\,$m$/$s, which gives an upwind shear stress of $\tau_{0}=0.196\,$kg$/$m2. As shown previously, this wind shear velocity is a representative value of $u_{\ast}$ above the threshold for saltation in real dune fields [30], although we note that on the field the wind strength has a strongly unsteady behavior and may vary substantially over the time [31, 32].

We see in Fig.3a the rescaled wind shear velocity $u_{\ast}/u_{{\ast}0}$ as a function of the rescaled downwind position $x/h_{\mathrm{f}}$, for different porosities $\Phi$. The fence is at the position $x=0$. As we can see, there is a strong decrease of $u_{\ast}$ as the wind approaches the fence from the upwind. This strong decrease is expected as the fence poses an obstacle to the wind thus extracting aeolian momentum. Moreover, the shear velocity decreases further in the fence’s wake until a minimum (denoted here $u_{{\ast}{\mathrm{min}}}$) is reached, whereupon upwind flow conditions are recovered later downwind. Similar results have been obtained before experimentally [33]. We see from Fig. 3a that very low porosities may lead to strong negative $u_{{\ast}{\mathrm{min}}}$ which means that backward flow occurs in the wake zone. We show in Fig. 3b the dependence of $u_{{\ast}{\mathrm{min}}}$ on $\Phi$ and $u_{{\ast}0}$. To the best of our knowledge, this is the first time that the three-dimensional diagram of Fig. 3b is computed. Such a diagram is useful for instance to predict under which conditions of $u_{{\ast}0}$ a fence of given porosity will lead to backward flow or which porosity is necessary to reduce the shear velocity to a pre-determined level below a given $u_{{\ast}0}$.

To get quantitative insight into how to use sand fences for large-scale soil protection, we extend the CFD calculation discussed above to investigate the airflow over an array of many fences. In what follows, we thus focus on the results regarding the array of 10 fences shown in Fig. 2.

Figure 4 shows the wind shear velocity as a function of the rescaled downwind position $x/L$, where $L$ is the spacing. The first fence is at the position $x=0$ and there are in total 10 fences at different spacing values (see legend). In the simulation of Fig. 4, the height of the fences is $h_{\mathrm{f}}=1.25\,$m and the porosity is $\Phi=40\%$. The wind shear velocity upwind of the fences is $u_{{\ast}0}=0.4\,$m$/$s.

As we can see from Fig. 4, after the strong decrease in the shear velocity upwind of the first fence, the behavior of $u_{\ast}(x)$ depends very much on the spacing. In particular, the downwind position at which $u_{{\ast}{\mathrm{min}}}$ is reached after the first fence varies strongly with $L$ (see Fig. 4). Moreover, due to the presence of fences downwind, upwind flow conditions are not achieved within the array. Instead, a maximal wind shear velocity $u_{{\ast}{\mathrm{max}}}$ is reached between each pair of neighbouring fences, whereas $u_{{\ast}{\mathrm{max}}}$ is smaller than the upwind shear velocity $u_{{\ast}0}$. We see that $u_{{\ast}{\mathrm{max}}}$ increases with $L$, which is expected since the larger the fences’ spacing the larger the fetch distance available for the wind flow to achieve higher speeds.

Moreover, we see in Fig. 4 that $u_{{\ast}{\mathrm{max}}}$ increases substantially from the first to the fourth fence for all values of $L$ investigated. However, further downwind $u_{{\ast}{\mathrm{max}}}$ increases much more slowly distance. We have found, that for a smaller fence porosity ($20\%$), the values of $u_{{\ast}{\mathrm{max}}}$ after the fourth fence are nearly constant with downwind position (that is, with fence number; see Fig. 5). In particular, Figs. 4 and 5 suggest that studies based on one to three fences cannot be used to understand the flow profile over arrays of more fences, as the shear velocity profile over the first three fences is very different from the profile further downwind. We see in Fig. 4 that $u_{{\ast}{\mathrm{max}}}$ increases approximately by a factor of two from the second to the last pair of neighbouring fences.

To investigate the characteristics of the flow over large-scale dune fields in the presence of fences, we now focus on the flow profile far downwind in the fences array. Specifically, we consider the results of the shear velocity between the ninth and tenth fences. Figs. 6 and 7 show the value of $u_{{\ast}\mathrm{max}}$, rescaled by the threshold shear velocity $u_{{\ast}\mathrm{t}}=0.25\,$m$/$s (consistent with medium sand [4]), as a function of $L/h_{\mathrm{f}}$ for different fences’ height, considering a porosity of $40\%$ and $20\%$, respectively. We note that $u_{{\ast}{\mathrm{t}}}$ may vary much from field to field depending on soil composition, grain size, the presence of non-erodible elements, humidity and the influence of moisture [1, 34]. From Figs. 6 and 7, it is possible to see for which range of fence height the shear velocity near the surface will not exceed the minimal threshold for transport, that is, ${u}_{{\ast}\mathrm{max}}/u_{{\ast}\mathrm{t}}<1$, which means total protection of the soil against erosion. The insets in Figs. 6 and 7 show ${u}_{\ast\mathrm{max}}/u_{{\ast}\mathrm{t}}$ as a function of $h_{f}$ for a fixed value of $L/h_{\mathrm{f}}=15$. As we can see, ${u}_{\ast\mathrm{max}}/u_{{\ast}\mathrm{t}}$ displays, in both insets, a maximal value at $h_{f}\approx 1\,$m and a minimum at $h_{f}=0.5\,$m, notwithstanding the different values of porosity associated with each case. Moreover, we note that in the regime of $h_{f}<0.5\,$m, surface effects become increasingly important as $h_{\mathrm{f}}$ decreases, thus affecting the behavior of ${u}_{\ast\mathrm{max}}/u_{{\ast}\mathrm{t}}$.

As we can see from Figs. 6 and 7, the critical value of $L/h_{\mathrm{f}}$ below which total protection against erosion is ensured, that is below which ${u}_{{\ast}\mathrm{max}}/u_{{\ast}\mathrm{t}}<1$, depends on the fence height. Based on the results of Figs. 6 and 7, we now investigate what is the maximal spacing between fences of given height $h_{\mathrm{f}}$ that guarantees no soil erosion by wind (${u}_{{\ast}\mathrm{max}}/u_{{\ast}\mathrm{t}}<1$). We call this maximal allowed spacing $L_{\mathrm{t}}$. The result for $L_{\mathrm{t}}$ is shown in Fig. 8a as a function of $h_{\mathrm{f}}$ for several porosities and different wind shear velocity thresholds for sustained transport, $u_{{\ast}{\mathrm{t}}}$. While the values $u_{{\ast}{\mathrm{t}}}=0.22\,$m$/$s and $0.25\,$m$/$s are consistent with fine and medium sand, respectively, we also performed calculations with a significantly larger $u_{{\ast}{\mathrm{t}}}$ ($0.32\,$m$/$s) to model enhanced resistance to mobilization due to stabilizing agents, such as moisture [34]. We see that $L_{\mathrm{t}}$ increases with $h_{\mathrm{f}}$ regardless of $\Phi$ and $u_{{\ast}{\mathrm{t}}}$, which can be understood by noting that the higher the fences the larger the sheltered distance.

We address the question of which is the most efficient fence array, that is which value of $h_{\mathrm{f}}$ should be used to protect an area of given size under the constraint of using the minimal amount of material to construct the fences. To address this question, we introduce the following cost function,

[TABLE]

where $S$ is the total downwind distance of the area to be protected. Note that $S/L_{\mathrm{t}}$ gives the length of the target field in units of numbers of fences. Fig. 8b shows the ratio ${\mathcal{C}}/S$. As we can see from this figure, for all studied values of $\Phi$ and $u_{{\ast}{\mathrm{t}}}$, ${\mathcal{C}}/S$ displays a minimum at $h_{\mathrm{f}}\approx 50\,$cm. This fence height is thus the optimal fence height to achieve total protection of a given soil area while ensuring minimal cost. Moreover, we see that, independent of the studied $\Phi$ and $u_{{\ast}{\mathrm{t}}}$, ${\mathcal{C}}/S$ has a maximum at around $1.25\,$m. This is a surprising result especially considering that the height of fences is often chosen to be $1.0\,$m [4]. Our result suggests the need for revisiting this choice. Note that our study concerns total protection (no erosion) of a large-scale dune field, where a serial array of multiple fences is used. Our results thus do not apply to an array with less than four fences, because, as shown in Figs. 4 and 5, the wind shear velocity profile over the first three fences is a transitional one. After the fourth fence, the maximal shear velocity between pairs of neighboring fences is much larger than in the transitional zone and increases only slowly downwind.

We also see in Fig. 8b that the value of ${\mathcal{C}}$, for a given $S$, becomes very large as $h_{\mathrm{f}}$ decreases down to values smaller than $50\,$cm. In particular, in this range of small $h_{\mathrm{f}}$, ${\mathcal{C}}$ increases with decreasing $h_{\mathrm{f}}$. This can be explained by the fact that such small fences are rather inefficient for soil protection as their wake region is too short. Moreover, we also see a decrease in ${\mathcal{C}}$ with $h_{\mathrm{f}}$ as the fence height exceeds $1.5\,$m. However, we note that using such large fences is not recommendable as it requires more effort to fixate them in the soil compared to the $50\,$cm ones. We conclude that there is an optimal fence height, which is around $50\,$cm, to guarantee total protection of soils against erosion with a large-scale array of multiple fences.

Discussion

We have shown, by means of CFD modeling, that the insights gained from studies using one to three fences cannot be applied to large-scale soil protection where an array of multiple fences is used. Our simulations show that the region corresponding to the first three fences is rather a transitional zone, in which considerable reduction of the shear velocity is achieved. However, after the third fence, the maximal shear velocity between pairs of neighbouring fences is much larger than in the transitional zone and increases only slowly downwind. Future investigations of the effect of sand fences on aeolian erosion for application in large-scale dune fields should thus consider at least four fences. Here, we have considered an array of 10 fences and investigated the shear velocity profile between the last two fences. Our calculations showed that minimal fabrication costs can be achieved if fences with height around 50 cm are used. This is the optimal fence height in an array of multiple fences to guarantee that no erosion occurs in the area to be protected. It is remarkable that this optimal fence height is independent of porosity $\Phi$ and, in particular, of typical threshold wind shear velocity $u_{{\ast}{\mathrm{t}}}$ — while both parameters affect $L_{\mathrm{t}}$ as shown in Fig. 8a. This means that the optimal array of fences applies both for mobile dune sand and for a terrain containing stabilizing elements or moisture.

The present work should be continued by computing the evolution of the sediment landscape in presence of fences. To this end, a morphodynamic modeling tool to simulate aeolian dune formation and migration [35, 36, 30, 37, 38, 39] should be coupled to the CFD simulation in order to model the erosion and deposition patterns resulting from the wind field over the terrain. As a matter of fact, here we have investigated the flow over a flat terrain covered with fences, but local topography evolves in time as sand is deposited in the areas between the fences. Therefore, the results of the present work apply to the initial surface conditions where the terrain is not covered with dunes. Moreover, fences are often applied in combination with the cultivation of vegetation, which acts as sand stabilizer thus helping fixate the soil [4]. It is thus important to include the effect of plants on the wind flow in future studies as well as to compute the topography resulting from the combined action of wind-blown sand and vegetation by extending the model presented in Refs. [40, 41, 42] to include the fences.

We remark that, while in our calculations we have assumed constant wind speed, in reality the wind velocity is varying over time which means that occasionally the assumed $u_{\ast}=0.4\,$m$/$s is exceeded. Calculations using unsteady winds would help to shed light on the characteristics of wind erosion under real conditions. Moreover, our calculations considered the two-dimensional soil profile in longitudinal direction, whereas three-dimensional flow effects [43, 44] are certainly important if the wind does not hit the fence perpendicularly. Three-dimensional CFD simulations should be thus performed to obtain quantitative insights into such effects. In particular, many different types of fence arrays and other obstacle geometries, such as placing the fences in zig-zag [24] or checkerboards [45] are in use. We thus hope that our CFD modeling will inspire future work to calculate the flow over such more complex (three-dimensional) geometries.

Methods

In the simulations, the fluid (air) is regarded as incompressible and Newtonian, while the average turbulent wind field over the soil is calculated as described in Refs. [46, 47]. The FLUENT Inc. commercial package (version 14.5.7) is adopted to solve the Reynolds-averaged Navier-Stokes equations, whereas in the computations the standard $\kappa-\epsilon$ model is applied to simulate turbulence.

The boundary conditions are modeled as follows. At the inlet of the channel, the logarithmic wind profile (Eq. (1)) is imposed, where the shear velocity $u_{{\ast}0}$ at the inlet is the only parameter of Eq. (1) varied in the calculations. A constant pressure ($P=0$) is applied at the oulet of the channel in order to produce a pressure gradient in flow direction. Moreover, we apply a non-slip boundary condition to the entire fluid-solid interface comprising the soil and the fences, while the shear stress of the wind at the top wall is set equal to zero (see also Refs. [46, 26, 48, 49]).

The time-averaged (or Reynolds-averaged) Navier-Stokes equations for the wind flow over the terrain are solved in the fully-developed turbulence regime. The standard $k-\epsilon$ model is used, and the default pressure-velocity coupling scheme (“SIMPLE”) of the solver is applied with its preselected values of parameters, as well as with the default option “standard wall functions” (see Ref. [47]). In particular, this option applies the wall boundary conditions to all variables of the $k-\epsilon$ model that are consistent with Eq. (1) along the channel’s bottom wall [50]. A second-order upwind discretization scheme is applied to the momentum, whereas for the tubulent kinetic energy and turbulence dissipation rate we apply a first-order upwind scheme [47]. A square grid with mean spacing of about $0.05h_{\mathrm{f}}$ is applied in the region close to the fence-fluid interface, as well as in the wake region at the front, upward and behind each fence close to the soil. However, this grid is larger in areas which are far away of the fences.

To solve the transport equations for the standard $k-\epsilon$ model [51], the following initial conditions are applied: The pressure and velocity are set to zero for all values of $x$ and $z$, while at the left wall ($x=0$), the logarithmic profile Eq. (1) is imposed. The convergence criteria for the numerical solution of the transport equations are defined in terms of residuals. These residuals provide a measure for the degree up to which the conservation equations are satisfied throughout the flow field. In the present work, convergence is achieved when the normalized residuals for both $k$ and $\epsilon$ fall below $10^{-4}$, and when the normalized residuals for both velocity components fall below $10^{-6}$.

Acknowledgements

This work was supported in part by CAPES, CNPq and FUNCAP (Brazilian agencies), the Brazilian Institute INCT-SC, by the German Research Foundation (DFG) Grant RI 2497/3-1 and by ERC Advanced grant FP7-319968 FlowCCS of the European Research Council.

Author contributions statement

I.A.L., A.D.A., E.J.R.P., J.S.A.Jr and H.J.H. designed the research, analyzed the data and wrote the paper.

Additional information

Competing financial interests: The authors declare no competing financial interests.