**Soil Erosion Dust Control and Sand Stabilization**

Itzhak Katra Edited by

Printed Edition of the Special Issue Published in *Applied Sciences*

www.mdpi.com/journal/applsci

**Soil Erosion**

## **Soil Erosion**

## **Dust Control and Sand Stabilization**

Editor

**Itzhak Katra**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Itzhak Katra Department of Geography and Environmental Development, Ben Gurion University Israel

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/soil erosion dust sand aeolian).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-03943-889-1 (Hbk) ISBN 978-3-03943-890-7 (PDF)**

© 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

**Itzhak Katra** received his PhD degree at the Department of Geography, Bar Ilan University, Israel. He was a postdoctoral fellow at the Division of Earth Sciences and Ecosystem Sciences, Desert Research Institute (DRI), Nevada, USA. Katra is a Professor at the Department of Geography and Environmental Development, Ben-Gurion University of the Negev, Israel, and the Head of the Aeolian Simulation Laboratory. The overriding aim of his research is to understand the dynamic processes of soil erosion in arid areas, sand transport, dust emission, and related air pollution. Katra has published papers in journals of earth sciences, environmental sciences, multidisciplinary sciences, and applied sciences.

## **Preface to "Soil Erosion"**

Soil erosion by wind is significant to Earth systems and human health. Soil-derived dust particles with origins in various source areas constitute one of the major components of global aerosols. There is a strong interest in understanding the factors and processes of soil erosion by wind as well as in developing and applying methods to control dust emission from soils and to stabilize active sands. The Special Issue contains information on applications of natural and synthetic materials to reduce soil erosion, development of materials and methods, experimental methods and modeling, impacts on the soil quality and the environments, and quantification of the efficiency in dust control and sand stabilization applications. Eight papers were accepted for publication, namely six research papers, one review paper, and one technical note.

> **Itzhak Katra** *Editor*

### *Editorial* **Soil Erosion: Dust Control and Sand Stabilization**

#### **Itzhak Katra**

Department of Geography and Environmental Development, Ben Gurion University, Beersheba 8410501, Israel; katra@bgu.ac.il

Received: 4 November 2020; Accepted: 9 November 2020; Published: 13 November 2020 -

**Abstract:** This Special Issue on soil erosion invites novel and original articles based on physical and chemical theories, field and laboratory experimental, soil analyses, and/or statistical and mathematical modeling that advance our knowledge on dust control and sand stabilization.

**Keywords:** aeolian processes; arid areas; dust emission; dust sources; environmental pollution; infrastructures; human activities; particle size distribution; polymers; sand dune; sand transport; soil erosion; soil quality

#### **1. Soil Erosion by Wind**

Soil erosion by wind is significant to Earth systems and human health, e.g., [1–3]. Soil-derived dust particles with origins in various source areas constitute one of the major components of global aerosols. Annual global dust emissions from soils into the atmosphere are estimated to be as high as 3000 million tons, including particulate matter (PM) that is less than 10 micrometers in diameter (PM10) [4]. Climate change of drier conditions is associated with desertification and, thus, increased dust emission from soils and sand-dune transport. Moreover, many soils throughout the world are subjected to the impacts of rapid population growth and extensive land uses, including agricultural fields, grazing areas, unpaved roads, mines and quarries, waste soils, active sand dunes and sand sheets, and more. There is a strong interest in understanding the factors and processes of soil erosion by wind as well as in developing and applying methods to control dust emission from soils and to stabilize active sands. This Special Issue on soil erosion invites novel and original articles based on physical and chemical theories, field and laboratory experimental, soil analyses, and/or statistical and mathematical modeling that advance our knowledge on dust control and sand stabilization.

#### **2. Diverse Impacts and Solutions in Soil Erosion**

This Special Issue was introduced to collect the latest research on relevant topics to address present challenging issues in dust control and sand mobilization. The Special Issue contains information on applications of natural and synthetic materials to reduce soil erosion; development of materials and methods; experimental methods and modeling; impacts on the soil quality and the environments; quantification of the efficiency in dust control and sand stabilization applications. Eight papers were accepted for publication; six research papers, one review paper, and one technical note. The review paper of Lal [5] provides us with a complete picture of soil erosion and gaseous emissions. The large magnitude of annual erosion of soil organic carbon has severe adverse impacts on soil quality and functionality, and emission of multiple greenhouse gases (GHGs) into the atmosphere. Three papers focus on sand dunes. The paper of Bird et al. [6] was aimed to investigate the temporal trends of four taxonomic groups to determine the effect of vegetation removal on dune assemblages over a 12-year period. They show that fixed dune treatment had very little effect, while a stronger response was found in semi-fixed treatments in particular for mobile dune indicator species. The paper of Yang et al. [7] is about the characteristics of the aeolian dune, wind regime and sand transport in

the Hobq Desert, China. Their work provides a scientific basis for the prevention and treatment of regional sand disasters. The paper of Wang et al. [8] on wind tunnel measurements of surface shear stress on an isolated dune downwind a bridge highlights the possible impacts on sand dune transport due to civil infrastructures. The other papers in this Special Issue focus on arid and semi-arid soils. The paper of Cheng et al. [9] is about the shearing behavior of the loess and post-harvest waste (PHW) mixture using small-scale and large-scale direct shear tests. Their work provides us with information on the ability of the loess-PHW mixture to resist seepage force and thus soil erosion on slopes. In the paper of Raveh-Amit and Tsesarsky [10], the biostimulation in desert soils for microbial-induced calcite precipitation (MICP) is a soil amelioration technique to prevent desertification and soil erosion. The results of their work demonstrate that biostimulated MICP is feasible in low-carbon mineral topsoils. The paper of Katra [11] was aimed to fill a clear gap in the efficiency of common product applications for reducing dust emission in quarry roads. The results of the wind tunnel experiments indicate that Hydrous magnesium chloride (Brine) was the most efficient compared with synthetic and organic polymers. The paper of Hanegbi and Katra [12] is a technical note on the development of a clay-based geopolymer for dust control and soil stabilization in semi-arid loess.

#### **3. Future Advances in Dust Control and Sand Stabilization**

Dust emission and sand mobilization can be reduced by using various applications. Products for dust control are based mainly on synthetic or natural polymers, which are applied by wetting the soil surfaces. A wide range of the products has been tested for dust emission by human activity such as mining and vehicles traveling on unpaved roads. Yet, there is a lack of fundamental research examining the efficiency of diverse products in the suppression of dust emission by wind. Moreover, further study is needed to investigate the possible environmental impacts of the diverse dust suppression substances, including the toxicity of atmospheric particulate matter when dust is emitted from the treated soils and/or soil-groundwater pollution because of vertical fluxes of the applied solutions on the surfaces. We aim to report on the future advances in the Special Issue "Soil Erosion: Dust Control and Sand Stabilization", Volume 2.

**Funding:** This research received no external funding.

**Acknowledgments:** This Special Issue is a result of a long-term work by all the authors, the reviewers, and the editorial team. A special thanks to Karena Pan, Section Managing Editor, Applied Sciences, MDPI.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Soil Erosion and Gaseous Emissions**

#### **Rattan Lal**

Carbon Management and Sequestration Center, the Ohio State University, Columbus, OH 43210, USA; lal.1@osu.edu

Received: 27 March 2020; Accepted: 14 April 2020; Published: 17 April 2020

**Abstract:** Accelerated soil erosion by water and wind involves preferential removal of the light soil organic carbon (SOC) fraction along with the finer clay and silt particles. Thus, the SOC enrichment ratio in sediments, compared with that of the soil surface, may range from 1 to 12 for water and 1 to 41 for wind-blown dust. The latter may contain a high SOC concentration of 15% to 20% by weight. The global magnitude of SOC erosion may be 1.3 Pg C/yr. by water and 1.0 Pg C/yr. by wind erosion. However, risks of SOC erosion have been exacerbated by the expansion and intensification of agroecosystems. Such a large magnitude of annual SOC erosion by water and wind has severe adverse impacts on soil quality and functionality, and emission of multiple greenhouse gases (GHGs) such as CO2, CH4, and N2O into the atmosphere. SOC erosion by water and wind also has a strong impact on the global C budget (GCB). Despite the large and growing magnitude of global SOC erosion, its fate is neither adequately known nor properly understood. Only a few studies conducted have quantified the partitioning of SOC erosion by water into three components: (1) redistribution over land, (2) deposition in channels, and (3) transportation/burial under the ocean. Of the total SOC erosion by water, 40%–50% may be redistributed over the land, 20%–30% deposited in channels, and 5%–15% carried into the oceans. Even fewer studies have monitored or modeled emissions of multiple GHGs from these three locations. The cumulative gaseous emissions may decrease at the eroding site because of the depletion of its SOC stock but increase at the depositional site because of enrichment of SOC amount and the labile fraction. The SOC erosion by water and wind exacerbates climate change, decreases net primary productivity (NPP) and use efficiency of inputs, and reduces soils C sink capacity to mitigate global warming. Yet research information on global emissions of CH4 and N2O at different landscape positions is not available. Further, the GCB is incomplete and uncertain because SOC erosion is not accounted for. Multi-disciplinary and watershed-scale research is needed globally to measure and model the magnitude of SOC erosion by water and wind, multiple gaseous emissions at different landscape positions, and the attendant changes in NPP.

**Keywords:** global carbon budget; soil organic carbon erosion; deposition; gaseous emissions; enrichment ratio; soil depletion; preferential removal

#### **1. Introduction**

As a natural geological process, soil erosion over eons has created the world's most fertile alluvial and aeolian (loess) soils. Acceleration of the natural erosion process by human activities, ever since the dawn of settled agriculture ~12 millennia ago, has caused the most severe environmental problems of the 21st century. Soil erosion, involving breakdown and transport of soil particles, requires energy, and a specific type of erosion depends on the source of energy (Figure 1). Water and wind are among the principal sources of energy, and thus major factors of erosion. Being a selective process, soil erosion removes and transports fine (i.e., clay and silt) and light (soil organic carbon or SOC) fractions. These three constituents (i.e., clay, silt, and SOC) are also key determinants of soil quality, and its capacity to provide numerous ecosystem services (ESs). However, these essential constituents

are depleted over time in soils prone to accelerated erosion. The latter has plagued the Earth and humanity for millennia. The data based on the analysis of sediments from 600 lakes worldwide show that anthropogenic activities accelerated global soil erosion 4000 years ago [1]. Many once-thriving civilizations vanished because they treated their soil like dirt [2,3]. The current problem of accelerated soil erosion is driven by a rapid and an indiscriminate expansion of agroecosystems for feeding the growing population. Further, the problem of soil erosion is also exacerbated by anthropogenic global warming [4,5]. In addition to adversely impacting the wellbeing of 3.2 billion people [6], accelerated soil erosion is also polluting the environment (i.e., soil, water, and air). It affects and is affected by the present and will be aggravated by the projected climate change.

**Figure 1.** Types of soil erosion driven by source of energy.

Under natural ecosystems, SOC stock is a sink of atmospheric carbon dioxide (CO2), and is protected against microbial processes through the formation of organo-mineral complexes and stable structural units or aggregates. Conversion of natural to managed ecosystems disrupts aggregates, exposes the hitherto protected SOC, and increases its vulnerability to transport by erosion and decomposition by microbial processes. Preferentially removed light SOC fraction is redistributed over the landscape, deposited in channels and transported to aquatic ecosystems and depressional sites (Figure 2). The labile SOC fraction is exposed to microbial processes when being transported, and following after redistribution and deposition phases of the erosion process. Furthermore, the historic land use based on extractive farming practices also mined off the SOC stock as a source of plant nutrients. Thus, soils of most agroecosystems are depleted of their original SOC stock. Consequently, soil quality is degraded, the capacity to perform ESs is impaired, and the environment (i.e., soil, water, air, and biodiversity) jeopardized.

**Figure 2.** The fate of soil organic carbon transport by erosion. About 40%–50% may be redistributed over the land, 20%–30% may be deposited in channel, 5%–15% may be carried into the ocean, and about 15%–20% may be emitted into the atmosphere. However, the exact partition may vary among soil, climate, land use, and other site-specific factors. Whereas the cumulative emission of CO2 may decrease at the eroded site, it may increase at the transported and depositional zones.

The global magnitude of historic depletion of SOC by all processes may be as much as 135 Pg C [7]. Consequently, degraded and depleted soils also have a large carbon (C) sink capacity to reabsorb atmospheric CO2 into SOC stock upon conversion to a restorative land use and adoption of conservation-effective practices.

It is this potential of restoring the global SOC stock, for advancing food and climate security and strengthening soils' capacity to provide ESs, that sustainable soil management is receiving the attention of policymakers. Ever since the launch of the 4 Per Thousand (4P1000) initiative at COP 21 in Paris in 2015 [8], world soils have been on the global agenda as an option to sequester C and mitigate global warming. Such initiatives are aimed at achieving greenhouse gas (GHG) neutrality through low-carbon farming [9].

Transport of C by accelerated soil erosion at a global scale is one such process that impacts the emission of CO2, methane (CH4), and nitrous oxide (N2O). The drastic increase in SOC erosion by anthropogenic activities poses a daunting challenge of assessing its impact on the global C budget (GCB) and GHG emissions. Therefore, it is important to credibly assess the mean annual flux of GHGs from soils during different erosional phases so that the magnitude of the carbon dioxide equivalent (CO2 eq) can be estimated. Whereas the soil C transported by erosional processes comprises of SOC and soil inorganic C (SIC), the fate of SOC transported by water and wind erosion that impacts the emission of GHGs [10] is not understood. Therefore, the objectives of this article are to describe the effects of erosion on the emission of GHGs into the atmosphere, explain processes affecting gaseous emissions by soil erosion, describe generic options that can reduce risks of soil erosion and minimize the emission of GHGs, and identify researchable priorities. This article is based on the hypothesis that accelerated soil erosion is a source of major GHGs including CO2, CH4, and N2O during all three phases of the erosional process.

#### **2. Materials and Methods**

The literature is replete with articles on soil erosion by water and wind. Thus, the literature search was specifically focused on available information on the magnitude of SOC transported by water and wind erosion was collated from the Web of Science, Google, and other sources. The literature search involved journals dealing with basic and applied sciences. The focus included journals dealing with: (a) earth sciences such as Global Change Biology, Global Biogeochemical Cycles, Biogeosciences, Geomorphology, J. Geophysical Research, Earth Surface Processes and Landforms, Geochemistry, J. Geophysical Res., J. Hydrology, Aeolian Research, (b) popular journals such as Science, Nature, Philosophical Transactions of Royal Society, (c) environmental sciences including Env. International, Climatic Change, Ecosphere, (d) journals devoted to soil science including Soil Research, Geoderma, Soil Sci. Soc. Amer. J., Catena, European J. Soil Sci., Australian J. Soil Res., J. Soil and Water

Conservation, and (e) those dealing with policy issues such as Land Use Policy, Science Policy, and Land Degradation and Development. Only those articles were selected for discussions in the present review which contained quantitative data on the magnitude of SOC or total carbon (TC) transported by erosional processes, and information on gaseous emissions at different landscape positions within an eroding landscape. While the literature searched is global, most of the articles addressing this theme were those published from the research done in the U.S.A., Europe, East Asia, Australia, and South America.

#### *2.1. Soil Erosion by Water: Transport, Redistribution, and Deposition of Soil Organic Carbon Over the Landscape*

Water erosion affects as much as 1.1 billion hectares (B ha) of the land area [11]. Available data on the magnitude of sediment load transported by world rivers are more credible [12] than that for the amount of soil moved by aeolian processes. The global land–ocean flux of sediment has reportedly increased from 14.0 Pg/yr. (Pg = peta gram = 1 billion metric ton) during the pre-human era to the contemporary flux in the absence of reservoir trapping to 36.6 Pg/yr. [12]. Sediments are enriched in SOC, and the global increase in sediment load may cause a strong increase in the transport of SOC, whose fate must be understood in relation to emissions of GHGs. Soil erosion on U.S. cropland increased by ~17% over the 20th century through the expansion of the land area under agriculture [13]. The SOC fraction entrained in the shallow runoff is moved and redistributed over the landscape. Erosion of soil and SOC stock has direct and indirect effects on soil and environment quality, net primary productivity (NPP), and efforts to achieve land degradation neutrality or LDN (Figure 3). The magnitude of the effect of emissions of GHGs is governed by the pathways of SOC erosion. The fate of SOC being redistributed depends on how it is being moved by the fluvial processes and on the temperature and moisture regimes at the redistribution and depositional positions (Figure 2). Quantitative assessment of the movement of SOC over the landscape is essential to establishing the watershed level C budget [14] that can be scaled up to the river basin and eventually to regional, national, or global scale. The magnitude of SOC erosion by fluvial processes varies widely (Table 1) depending on a range of factors. Important among these are climate [10,13], soil [15–17], terrain [18] and land use [13–15,19–23]. On the basis of some empirical data from 240 runoff plots studied over the entire rainy season from diverse global ecoregions, Mueller-Nedebock and Chaplot [18] estimated that the total amount of SOC displaced by sheet erosion from its source would be 1.32 ± 0.20 Pg C, or about 11.4% of the annual anthropogenic emission of 11.5 Pg C in 2019 [21]. Integrating all C fluxes for the EU agricultural soils, Lugato et al. The author of [24] estimated a net C loss or gain of −2.28 Tg CO2 e/yr. and +0.79 Tg CO2 e/yr., and they argued that strong agricultural policies are needed to prevent or reduce soil erosion. Assessing and accounting for all the additional feedback and C fluxes due to displacement by erosion, Lugato et al. [15] estimated a net source of 0.92 to 10.0 Tg C/yr. from agricultural soils in the European Union to the atmosphere over the period of 2016–2100.

**Figure 3.** Effects of accelerated soil erosion on the global carbon cycle and the increase in the daunting challenge of achieving land degradation neutrality.


**Table 1.** Examples of regional, national, or global terrestrial soil organic carbon (SOC) erosion by water and other processes.

The SOC being eroded is either deposited in the landscape, in the channel, or carried into the ocean (Figure 2). Some of the SOC being transported is emitted into the atmosphere as CO2 or CH4, depending on the degree of wetness or anaerobiosis. In China, Fang et al. [17] observed that 42% of the eroded SOC was redeposited within the catchment. The mean residence time (MRT) of the deposited C depends on a range of site-specific factors, and the fraction composition (labile, intermediate, passive) of the eroded SOC. Wang et al. [30] reported that cumulative emission of soil CO2 decreased slightly at the erosion site but increased by 56% and 27% at the transport and depositional zones, respectively, in comparison to non-eroded sites. Wang and colleagues concluded that overall, CO2 emissions contributed 90.5% of the total erosion-induced C loss over the 4-month experiment. Whereas buried SOC at depositional sites may have a higher MRT even for the fast and intermediate turnover pools [31], susceptibility to decomposition may be much higher for the labile fractions redistributed within the landscape.

Examples of the magnitude of SOC erosion by water from different regions are shown in Table 1 and may range from 1.1 to1.3 Pg C per year. The preferential removal of SOC by water erosion is indicated by a high enrichment ratio of SOC (and clay) in sediments compared to that of the surface soil from which the sediments originated. Consequently, the enrichment ratio for SOC in alluvial sediments is >1 and may be as high as 12 (Table 2). Erosional processes lead to a preferential transport of SOC because it has a low bulk density and is concentrated in the surface soil layer. In cases where sediments are derived from subsoil (i.e., gully erosion), the enrichment ratio can be less than 1 [32].

There are a few studies involving techniques of quantitative measurement of SOC/TC transported by erosion from a watershed or a well-demarcated area. In Australia, Chappell et al. [25] estimated the magnitude of 137Cs -derived redistribution of SOC by all processes (water, wind, and tillage) at 4 Tg SOC/yr for 1950–1990s. This represents an average loss of 2% of TC stock, assuming that total C is mineralized as CO2, and would represent a net national flux of 15 Tg CO2 eq/yr from all C pools in Australia. In a follow-up study, Chappell et al. [33] estimated the global terrestrial SOC erosion at 0.3–1.0 Pg C/yr. For the Seyhan River Basin in the Mediterranean region of Turkey, the rate of SOC erosion was estimated at ~0.2 Mg C/ha·yr [20].



#### *2.2. Wind Erosion*

Wind erosion, affecting about 550 million hectares of the global land area [11,39], is caused by aeolian (or eolian) processes. The term "aeolian" is derived from the Greek god "Aeolus", the keeper of the wind. Wind erosion is strongly affected by soil texture. Soils most susceptible to wind erosion may have <5% clay and <3% silt, and >50 cm deep surface layer [22]. Wind erosion may create 500–5000 Tg (million tons) of dust annually with a strong impact on soil properties, air quality, and human health [39–41]. The environmental impacts of wind erosion during the Dust Bowl Era of 1931 through 1939 are described by Steinbeck [42].

Accelerated erosion affects critical biotic and abiotic processes governing the soil/ecosystem C cycle. The magnitude of the loss of SOC by wind erosion is related to that of the fine soil fraction [43]. The loss of C-enriched fine soil particles depletes its SOC and reduces its future potential to restore the SOC pool. The aeolian erosion process affects both progressive and regressive pedogenesis in dry eco-regions. On agricultural lands, erosion degrades soil quality by removal of silt, clay, and SOC fractions through effective sorting processes that leave behind only coarse sand and gravels [23]. The loss of NPP reduces the plant feedback and aggravates the SOC loss [44].

The wind-blown dust is also enriched in SOC, which may also depend on soil texture. The global estimate of SOC erosion by wind may be as much as 0.3 to 1.0 Pg C/Yr (Table 3) and some highly vulnerable soils may lose 3.6 Mg C/ha per year [32]. Losses of SOC by wind erosion in Northern China are estimated at 0.9 Tg C/yr [45]. The loss of PM10 (particles of <10 micrometer) adversely impacts soil nutrient reserves [40]. Wind-blown dust is also enriched in SOC and has a high enrichment ratio (Table 4). In Niger, Sterk et al. [46] assessed nutrient and C losses in saltation and suspension transport by conducting chemical analysis of the trapped material at 0.05, 0.26, 0.5, and 2 m. The sediments were three times richer than topsoil at 0.5 m and 17 times at 2 m. In Australia, Webb et al. [47] observed that

the SOC-enrichment ratio ranged from 2.1–41.9 for a sandy and 2.1 for clayey soil (Table 4). Webb and colleagues hypothesized that in addition to particle size, distribution, and the degree of aggregation, size-selective sorting of SOC during transport may enhance the enrichment of SOC dust emissions. The SOC concentration in two of the dust samples was 15%–20% by weight. A study in China by Ravi et al. [48] documented that an increase in the particulate matter emissions (e.g., black earth) from biochar-amended soils may counteract the negative emissions potential of biochar. Magnitude of dust emitted is aggravated by human activities [40].


**Table 3.** Examples of regional, national, or global terrestrial SOC erosion by wind.

**Table 4.** The enrichment ratio of carbon in sediments derived from wind erosion.


Both direct and indirect effects of accelerated erosion, especially in arid and semi-arid regions, exacerbate the risks of desertification and drastically increase the already daunting challenge of achieving land degradation neutrality or LDN by 2030 [54,55]. Soil degradation impacts of accelerated erosion by wind, and its positive feedback to the process of desertification, have strong adverse consequences on Earth systems and human environments [40], as well as on NPP, the input of biomass-C into the soil, and on the disruption of the global C cycle. Thus, achieving LDN would necessitate the global adoption of conservation-effective measures that reduce risks of both aeolian and alluvial processes of soil erosion [56]. Soil restoration strategies must be directed towards increasing the input of biomass-C into the soil. Increase in NPP, and the attendant increase in the input of biomass-C into soil, would restore SOC stock [57]. Dust emission caused by wind erosion may be aggravated by the projected climate change. Thus, Duniway et al. [58] recommended multidisciplinary and multijurisdictional approaches and perspectives to understand the complex processes of dust emission and identify strategies of its mitigation.

#### *2.3. Gaseous Emissions from Eroded Sediments and the Fate of Carbon Transported and Deposited over the Landscape*

Soil C stock is an important component of the global C cycle. The historic C loss from soil may have emitted as much as 537 Pg C or 27% of the amount present before the onset of agriculture about 10 millennia ago [59]. Erosion and redistribution disturb a large quantity of soil C in managed and natural landscapes. The fate of soil C impacted by erosion may differ among the sites of erosion, redistribution, and deposition (Figure 2). Some examples of gaseous emissions from eroded and depositional sites are shown in (Table 5) [60,61]. Assuming an average flux of 300 mg CO2 eq/m2·h based on literature review, Oertel et al. [62] estimated the global annual net soil emissions at ≥350 Pg CO2 e, as compared with the 2018 anthropogenic emission of 42.1 Pg CO2 e [21]. However, the large emissions from soil C transported by aeolian and alluvial processes are not considered in the global C budget

(GCB), which creates a lot of uncertainty. Nonetheless, understanding, managing, and reducing the erosion-induced gaseous flux of CO2, CH4, and N2O (Figure 2) is an important researchable priority to reduce uncertainty in the GCB. It is also critical to identify hot spots (vulnerability, resilience, and action), and plans of targeted interventions for managing the flux [63]. Several pedological processes impacted by erosion also affect NPP through alterations in availability and uptake of water, nutrients, and photosynthesis. Credible assessment of C dynamics in agricultural and other landscapes is important to addressing global issues [19]. In a Mediterranean Seyhan river basin, Cilek [29] estimated SOC erosion of 0.163 Mg C/ha yr. (total SOC loss of 349, 850 Mg C/yr. over a total watershed area of 21,485 km2). Based on the assessment of nine river basins in China, Wang et al. [30] found that total SOC erosion was 68.4 and 77.3 Tg C/yr. for 1995–1996 and 2010–2012, respectively. Of this, 57% and 47% were redistributed over land, 25% and 44% was deposited in channels, and 18% and 8% were delivered into the sea, respectively. However, how much and which gases were emitted were not determined. For the period A.D. 1850–2005, Naipal et al. [27] estimated global SOC flux of 47 ± 18 Pg C, of which 79%–85% occurs on agricultural and grasslands.

**Table 5.** Examples of gaseous emission from eroded sediments and disrupted/broken aggregates by erosional processes.


#### *2.4. Implications of Ignoring Erosion Induced Transport of Carbon in Estimating Global Carbon Budget*

Erosion-induced transport of soil C (SOC and SIC) is a large and growing component with a strong impact on the GCB, but which is now being omitted [21]. However, the erosion-induced impact on soil C stock and flux, a large component comprising of multiple gases (CO2, CH4, N2O) and multiple processes (e.g., water, wind, gravity, and tillage), must be dually considered. High-resolution models [20] must be developed to improve methodological protocols to account for this serious omission. By credibly accounting for the effects of erosion on net C exchange between the soil and the atmosphere, it may be possible to identify global hot spots of undertaking targeted interventions to mitigate the erosion-induced positive feedback to global warming. In Australia, Chappell et al. [25] found SOC erosion by all processes at 4 Tg/yr. (or 2% of total C stock in 10-cm depth). Assuming that most of this is mineralized, Chappell and colleagues estimated a flux of ~15 Tg CO2 e/yr. representing 12% emissions from all C pools in Australia and concluded that it was an important source of uncertainty in the national carbon budget. By extending this study globally, Chappell et al. [50] estimated global terrestrial SOC erosion of 0.3–1.0 Pg C/yr, highlighted the significance of ignoring it in the accounting of the GCB, and suggested that accounting for SOC erosion would reduce uncertainty in the GCB.

#### **3. Conservation: E**ff**ective Measures for Reducing SOC Erosion and A**ff**orestation of Eroded Lands for Sequestration of Atmospheric CO2 in the Terrestrial Biosphere**

During the 1950s to 1990s, the objective of erosion control was to conserve soil and water, reduce the loss of soil fertility, and minimize the risks of non-point source pollution. Since the 1990s, two among important objectives of soil conservation and effective erosion control measures are to: (1) promote low-C agriculture (9), reduce risks of water runoff and soil erosion so that the attendant emission of GHGs can also be reduced from SOC erosion, and (2) to sequester atmospheric CO2 in soil and vegetation through the restoration of eroded soils and desertified ecosystems. Examples of technological options to accomplish these include the adoption of conservation agriculture including no-till farming with retention of crop residue mulch [4,64–67] and use of cover cropping during the off-season [68] that reduce the risks of water runoff [69], boost SOC stock for food and climate [8,67], and

curtail transport of SOC and nutrient-enriched sediments [70] for accomplishing objective 1. Similarly, promoting afforestation [71] and adopting the concept of "Reducing Emissions from Deforestation and Forest Degradation or REDD [72] through afforestation of degraded soils [73], and establishments of shelterbelts in areas prone to wind erosion [74–76] would be pertinent to advancing objective 2 of sequestration of atmospheric CO2 in the terrestrial biosphere.

#### **4. Conclusions**

The synthesis and a critical review of the literature presented above indicate that the hypothesis of the study is proven. Over and above the offsite effect of sedimentation and non-point source pollution, accelerated soil erosion is also an important source of major GHGs (CO2, CH4, and N2O) emitted during all three stages of the erosion process.

The discussion presented has also led to accomplishing the major objectives of the study:


**Funding:** This research received no external funding.

**Acknowledgments:** Support received from Maggie Willis and Kyle Sklenka, the staff of the Carbon Management and Sequestration Center of the Ohio State University, is appreciated and acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Can Vegetation Removal Successfully Restore Coastal Dune Biodiversity?**

#### **Tania Leah Fairfax Bird 1,\*, Amos Bouskila 2, Elli Groner <sup>3</sup> and Pua Bar Kutiel <sup>1</sup>**


Received: 16 January 2020; Accepted: 17 March 2020; Published: 28 March 2020

**Abstract:** Coastal dune habitats have been declining globally over the last several decades due to rapid urbanization. Within remaining dune systems, dune fixation has resulted in further losses of mobile dunes with negative impacts on their associated species. Some studies suggest vegetation removal can initially promote habitat heterogeneity, and increase availability of suitable habitats for psammophile, xeric and endemic mobile dune species, but longer-term responses are generally unknown. We investigated the temporal trends of four taxonomic groups to determine the effect of vegetation removal on dune assemblages over a 12-year period at an LTER site. Three different forms of removal are investigated here—removal in a grid form on fixed dunes, removal of the wind-facing slope vegetation on semi-fixed dunes and opportunistic off-road driving on disturbed dunes. Results were varied across taxa, highlighting the need for multi-taxa monitoring in conservation and restoration management. Overall, fixed dune treatment had very little effect, while a stronger response was found in semi-fixed treatments in particular for mobile dune indicator species, which showed evidence of recolonization within a few years following treatment. Disturbed dunes were most similar to mobile dunes for animal taxa indicating that pulse removal may not be as effective as continuous press disturbance. Nevertheless, a less destructive form of disturbance such as re-introduction of grazing might be preferable and requires further investigation.

**Keywords:** dune stabilization; restoration; coastal dune; vegetation removal; multi-taxa; biodiversity; LTER; temporal dynamics; shrub encroachment

#### **1. Introduction**

In the last 70 years there have been extensive losses of coastal habitats and naturally functioning coastal dune systems globally [1–7]. This is due to dramatic population expansion and urban development along coastal plains, coupled with changes in land-use practices and climate-related changes in aeolian processes. [1,8–10]. Globally this has led to substantial increases in dune fixation (stabilization) and a significant loss of naturally mobile (shifting) dunes [11–16]. Restoration strategies for coastal dune systems are complex because sandy dunes are exposed to multiple dynamic processes that operate on different temporal and spatial scales [17,18]. Moreover, dune restoration methods are varied and sometimes conflicting. These range from re-establishment of vegetation, to removal of vegetation and eliminating exotics; the strategy usually depends on the management objectives (reviewed by [2,19]). Many dune management regimes around the world have focused on promoting dune fixation through the introduction of (sometimes non-native) grasses and nourishment in order to enhance flood protection and wave erosion defences. These activities have further added to the loss of dune mobility in coastal systems, and their associated specialist species [20–25]. In this study we test the effect of vegetation removal to see whether it can be used as a tool to restore mobile dune biodiversity.

Flora and fauna found on mobile dunes are often psammophilic, endemic species, highly adapted to harsh, xeric sandy conditions, and often form unique assemblages [6,26–31]. Several researchers recognize that dune fixation, whether through succession or intentional re-establishment, can result in adverse consequences for mobile dune species [14,24,25,32,33]. In fact, a review of coastal management actions reported that dune fixation always had a negative impact on coastal species diversity, either through species loss or increase (e.g., invasive species), or changes in composition [19].

Most open habitats in Europe are semi-natural, including dry sandy grasslands and coastal dune systems, which have been traditionally maintained by anthropogenic activities for thousands of years [34]. Some argue that pastoral human activity simply replaced natural dynamic dune characteristics such as wind erosion, wildlife grazing and trampling, which had existed prior to the advancement of civilization [1,35]. Contemporary coastal researchers now advocate actively restoring mobile dunes by removing vegetation and introducing disturbance, in order to reactivate sand mobility and imitate natural disturbances [4,6,7,15,33,36–44].

A disturbance can generally be defined as any event that disrupts the assemblage or population structure, or changes the resource pool, substrate availability or physical environment [45,46]. It can therefore refer to either natural or anthropogenic events. Three main types of disturbances are commonly recognized; i) pulse disturbances are short-term and sharply delineated in time (e.g., natural flood, wildfire, or deforestation); ii) press disturbances may arise sharply and then reach a constant level that is maintained (e.g., grazing, persistent pollution, or off-road vehicle use), and iii) accumulating ramp disturbances, which occur when the strength of the disturbance steadily increases over time (e.g., drought or build-up of toxic waste; [47,48]. We investigate two pulse disturbances and one stress disturbance in this study.

Mobile coastal dunes are naturally exposed to high degrees of press disturbance in the form of wind erosion and sand burial, as well as pulse disturbances, such as blow-outs and sea storms [15,49,50]. As dunes shift inland, weakening of these processes allows perennial species, such as *Artemisia monosperma*, to establish. Once established, nutrients build up in the sand around their bases, facilitating the establishment of successional species, and in turn, leading to dune fixation [5,51,52]. New embryonic dunes forming at the beach would usually replace fixed dunes in a continuous cycle. Established dunes can also alternate between mobile and stable states following climatic perturbations or other disturbances in longer cycles [24,53]. Removal of disturbances from the system, such as removal of human activity or changes in wind and sand deposition, can lead to loss of this cycle. In Wales, rapid dune fixation occurred over the last 60 years after rabbits that had maintained dune perturbation were nearly all wiped out due to extermination and disease [15]. In New Zealand and California, shrub encroachment by invasive species resulted in extensive dune fixation [25,32].

Removal of perennial vegetation has been shown to be an effective tool for reinstating the geomorphological traits of naturally occurring mobile dunes [15,24,40,54–56]. This restoration approach may seem unconventional because a) it involves the intentional introduction of an anthropogenic disturbance, b) it aims to restore an earlier successional state, as opposed to passively allowing natural vegetation succession to continue. However, it is argued that sand dynamics in most coastal dune ecosystems no longer function naturally due to restrictions of sand movement along the beachfront and inland [57]. Therefore, human intervention and introduction of disturbance is needed to restore the dunes.

In general, ecological restoration projects are widely reported, but are rarely evaluated systematically, particularly on mid- and long-term temporal scales [2,58,59]. Lack of clearly defined and measurable restoration targets also makes evaluation of manipulation responses difficult [2,60–62]. Restoration can involve abiotic and biotic manipulations; if the goal of a restoration project is to restore

biodiversity and improve ecosystem functioning, the recovery of biotic components of ecosystem processes across multiple trophic levels must be understood [63–66]. Restoration of coastal dune systems should therefore consider the response of both faunal and floral species.

Biotic responses are often not reported in landscape restoration, and rarely for animal taxa [59,67]. Most coastal restoration projects only consider geomorphological and vegetation components, and usually are only monitored over a short time frame. Only a handful of coastal restoration studies evaluate more than one taxon concurrently, and a recent review [2] found that 88% of coastal restoration studies only considered plant responses, usually only in terms of vegetation cover or species richness, not composition. There is a clear need for a coastal restoration studies involving vegetation removal to consider responses of species composition over a long timeframe and across multiple taxa (particularly faunal taxa). We found only one study that investigated composition responses of multiple taxa (arthropods and plants) to disturbance in coastal dunes [49], while Kutiel et al. [68] are the only authors to consider animal and plant responses to vegetation removal concurrently. The Nizzanim Dune Nature Reserve (NDNR), on the southern Mediterranean coast of Israel, provides a rare opportunity to conduct a restoration study considering all these components.

Livestock grazing and firewood collection by nomadic tribes was prevalent for hundreds of years in what is now NDNR [69–71]. Exclusion of this anthropogenic presence since the middle of the last century has led to substantial shrub encroachment and dune fixation across the reserve [69,72]. It is expected that in the absence of disturbance, the dunes in NDNR will eventually homogenize and become fixed, which will result in the loss of mobile dunes and their associated species, and an overall loss of β- and γ-diversity across the reserve [26].

To mitigate the threat of dune fixation in NDNR, the Nizzanim Long Term Ecological Research (LTER) project implemented two restoration experiments involving perennial vegetation removal as a form of pulse disturbance. The goal of these trials was to restore the heterogeneity of dune states and increase γ-diversity by providing a suitable habitat for mobile dune species [69]. This is the largest landscape scale restoration program undertaken in Israel, and the long-term responses have been monitored as part of the LTER program since 2005.

Parallel to the trials, off-road vehicles (ORV) have been utilizing several dunes in the reserve providing an opportunity to examine the effect of a continuous press disturbance on species composition. In this study, we investigate the effect the two controlled pulse removals and the un-controlled ORV press disturbance, in order to determine whether these different types of vegetation removal recreate suitable habitat for mobile dune species. Restoration can take several years, even decades, and can rarely precisely replicate an original status [2,73]. Therefore, it was expected that by removing vegetation and creating disturbance, the dunes would become more suitable for mobile-dune obligate species, and less suitable for fixed-dune species.

Employing the terminology of Lake [73], we used a M-RCI experimental design; comparing multiple repeated replicates in Reference (mobile), Control (un-treated) and Intervention (treated) dunes. We addressed the following questions:

i. How do species compositions of different taxa respond to different removal treatments?

ii. Can indicator species represent assemblage level responses?

iii. Can the different forms of removal contribute to successful coastal dune management?

#### *Hypotheses Tested*

First, we tested the hypothesis that in control and reference dunes, fixed dunes would support a greater abundance (*H1)* and species diversity (*H2)* under productivity theory [74].

Next, we tested the responses to treatment predicting that:

*H3*: Reduction of primary productivity through mechanical removal will shift a) abundance, and b) diversity towards reference dune levels

*H4:* Composition of treated dunes will shift away from the control dune composition towards reference (mobile) dune composition for each taxon

Over time, responses in composition could therefore be expected to either:

*H4a*: Gradually become more similar to the reference dune state (turnover)

*H4b*: Present an initial shift and then gradually return to control state (resilience)

or

*H4c:* Remain the same and not shift away from the control state (resistance)

Natural coastal dunes are dynamic habitats exposed to a variety of press (wind erosion/sand burial, grazing) disturbances, which creates assemblages of highly adapted mobile-dune species. Therefore, we predicted that:

*H5:* Mechanical press disturbance (ORV disturbance) will be more effective than a single mechanical pulse disturbance (removal by bulldozer), at creating community composition similar to the reference dune.

Finally, we considered the impacts of treatment on individual Indicator Species (IS). Since IS will tend to be highly affiliated to their specific habitat type, species abundance should increase with habitat availability.

*H6:* Treatment of semi-fixed and fixed dunes will result in an increase in abundance of species that are affiliated to mobile dune habitat, while species that prefer semi-fixed and fixed habitats will decrease.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The Nizzanim LTER is an on-going collaborative project established in 2004 to monitor plants, arthropods, small mammals and reptiles [26]. The study was conducted at the Nizzanim LTER site in Nizzanim Dunes Nature Reserve (NDNR), Israel (31◦42'–31◦44'N, 34◦35'–34◦36'E) covering an area of 20 km2 along the Mediterranean coastline. Annual average temperature is 20 ◦C, and rainfall is 400-500 mm per annum, falling mainly during winter (November–April). The common wind direction is south-west with a very low Drift Potential Index [70].

Sand dunes in NDNR are a globally distinct ecosystem situated at the intersection between desert and Mediterranean coastal dune systems [69]. The edaphic conditions and the physiognomy of the sand dunes enable numerous Saharan taxa to range northward in spite of the temperate climate that favors Mediterranean species [75], creating an unusual and highly diverse assemblage compared to other Mediterranean coastal systems [76]. Three distinct fixation states of dunes are recognized in NDNR, each state with its own characteristic species diversity and composition [26,77]. Following the declaration of the reserve, livestock grazing and wood harvesting by local communities (predominantly Bedouin) were prohibited and rapid dune fixation occurred, increasing the shrub cover from 4% to 20% of the landscape over a period of 60 years [27,72,78].

The Nizzanim LTER site consists of mobile, semi-fixed, and fixed dunes separated by inter-dune depressions [79]. The three fixation states can be categorized based on the gradient of perennial percentage vegetation cover (PPC), aeolian sand movement, and visual indicators such as dune geomorphic structure, plant distribution, dominant perennial species and soil color [27,56,79–81].

Mobile dunes have 5–15% PPC, mainly *Ammophila arenaria*, and make up approximately 20% of the reserve [27]. Wind-driven shifting sands result in frequent burial and exposure of perennial plants on mobile dunes. Semi-fixed dunes have 16–30% PPC, are dominated by wormwood (*Artemesia monosperma*) and desert broom (*Retama raetam*), and cover approximately 70% of the area. Lastly, fixed dunes vary between 30–50% PPC, and compose approximately 10% of the area of the reserve.

#### *2.2. Removal Experiments*

Three separate treatments were tested (Figure 1). For each set of treated dunes, there were paired untreated control dunes and mobile dunes were considered as the restoration reference for all treatments (Figure 1e) (other dunes were monitored under the LTER but were not included in this study). A summary table of the overall experimental design has been included in the supplementary information (Table S1).

**Figure 1.** Dune treatments investigated in Nizzanim Long Term Ecological Research (LTER); (**a**) Fixed untreated (control) dune; (**b**) Fixed treated (**c**) Semi-fixed untreated (control); (**d**) Semi-fixed treated; (**e**) Mobile dune (reference) dune; (**f**) Off-road vehicle disturbance. Semi-fixed and Fixed Treated dunes were paired with un-treated Semi-fixed and Fixed untreated control dunes respectively. All treated dunes were compared to their un-treated control pair and Mobile dunes as the restoration reference. Disturbed dunes were exposed to illegal, un-quantified off-road vehicle disturbance. They were compared to both Mobile and Semi-fixed untreated dunes to examine their relative composition. Frames of the different dune types shown are colour consistent throughout all figures.


design was chosen in order to leave some vegetation needed for animal survival, while more closely emulating the distribution of natural vegetation found on mobile dunes [see 77]. Four untreated semi-fixed dunes were selected as controls (Figure 1c).

iii. Disturbed Dunes Treatment: Off-road vehicle (ORV) disturbance has been occurring illegally in the reserve for many years (Figure 1f), which began prior to the initiation of the experimental plots. We began monitoring these dunes in 2012, in order to study the impact of ORVs on dune biodiversity. We consider these dunes as 'treated' with an unquantified press disturbance. The exact natural state of these dunes prior to the disturbance is unknown; on the one hand, one might assume they were similar to mobile dunes based on historical aerial photos from 1965 and 1999 [69]. On the other hand, most dunes in the immediate vicinity are semi-fixed. Disturbed dunes were therefore compared with both mobile and un-treated semi-fixed dunes as the reference dunes. Annual plants are almost entirely absent from disturbed dunes and as such this taxon was not monitored on these dunes.

#### *2.3. Data Collection*

#### 2.3.1. Rodents

Rodents were monitored in autumn (September–October) from 2005–2016, using 27 Sherman traps per dune in three rows of nine traps each, placed alternately in open and bush patches spaced approximately 10 m apart. Trapping took place over two consecutive nights, using mark recapture to calculate the estimated abundance for each species. We used the Chapman estimator, which is an adjusted form of the Lincoln–Petersen Index [82].

#### 2.3.2. Reptiles

Reptiles were sampled in autumn 2005–2016 using a combination of several methods, including ten pitfall traps (10 L buckets), two track transects (100 m long), four activity transects (100 m long) and opportunistic sightings. Ranked abundance for each species was then calculated using a max-log algorithm across all sampling methods with a range of 0 (absent) to 5 (highly abundant). For full calculation methods see Shacham and Bouskila's work [28].

#### 2.3.3. Beetles

Ground-dwelling beetles were monitored in spring (March–May), between the years 2006–2016 (not all years were sampled). Beetles were collected using dry pitfall traps of 12 cm depth and 10 cm diameter, in a regularized pattern of 10 traps on each dune section (windward slope, crest and slip-face), alternating between open and shrub patches. Pitfall traps remained open for 36 hours. The 30 traps were then pooled to give one dune sample.

Beetles were trapped live, identified to morphospecies (*sensu* [83]), and released where possible or collected and later identified by experts. Some species remained as morphospecies or were identified to genus or family, i.e., to the best Recognizable Taxonomic Unit (RTU) possible. Since data was collected using only one method, we were able to use a direct measure of abundance for beetles (compared to the rank abundance used with the mixed methods for reptiles, or estimated abundance for mark-recaptured rodents).

#### 2.3.4. Annual Plants

Annuals were monitored in spring (March–May) 2006–2016, using 40 × 40cm quadrats. As for beetles, all three sections of the dune were sampled, resulting in a total of 72 quadrats per dune. In each slope a 100 m2 plot was marked out, and within each plot 12 open and 12 bush quadrats were randomly placed. Total dune cover was then standardized based on the area contribution of each slope to the whole dune.

#### *2.4. Data Analysis*

For all analyses we tested each taxon separately. First, species with a total abundance of <5 individuals (or 5% cover for annuals) across all years and all dunes were removed. We also removed all species that had been found in less than three years, irrespective of their total abundance. Each experimental treatment type was then examined separately; for each we compared the treated dunes to the respective control dunes (e.g., treated and un-treated fixed dunes) and the reference dunes (mobile). In a similar manner, for the disturbed dunes we consider un-treated semi-fixed dunes as the control pair, and mobile dunes as the reference. All analyses were performed in R version 3.4.3 [84].

#### 2.4.1. Assemblage Abundance and Richness

To understand the effects of treatment on assemblages, we examined the response of assemblage abundance, species richness and composition. First, we compared control dunes (untreated semi-fixed and fixed dunes) with the restoration reference (mobile dunes). Then we used these differences to set the expectations for restoration effects. For example, if there were no differences between control and reference states, we expected no difference between treatment and control, or treatment and reference states (effectively *H1* would be the same as *Ho*). Conversely, if there were differences between control and reference, then treated dunes were expected to become significantly different from the control, and would either become more similar to the reference, or different from both control and reference.

For each treatment we tested the treated dunes against their untreated control and against the mobile dune reference. Estimated species richness, abundance of rodents, and percentage cover of annuals are all continuous, numerical variables. Therefore, to test differences in these parameters we used linear mixed-effect models as described in Laird and Ware [85]. Since dune sampling was not balanced across all years, we set year as a random intercept effect without interactions for any models, using the *nlme* package [86].

Generalized Linear Mixed Models (GLMM) allow the introduction of random effects in hierarchical, non-linear, count data [87]. Species richness of all taxa, total abundance of beetles, and rank abundance of reptiles are count data. Therefore, we tested these parameters for differences among dune states using GLMM with Penalized Quasi-Likelihood (GLMM-PQL) with the *MASS* R package [88], with *dune state* specified as a fixed effect and *year* as a random intercept effect, under a Poisson distribution.

#### 2.4.2. Composition

To test the effect of treatments on dune composition we used Redundancy Analysis (RDA) in the R package *vegan* [89], in a constrained model with year specified as a conditional factor:

$$\text{rda(species matrix)} \sim \text{dure.treat}^\* \text{ year} + \text{Conditional(year)} \tag{1}$$

This allowed us to examine the effect of dune type (i.e., treated, reference and control dunes) after the effect of year had been removed as a covariable. We then used the scores along the first RDA axis (RDA1) to calculate the treatment effect as a proportion of total distance between control and reference dunes, shown by directional arrows. ). Finally, we used pairwise Permanova to test the effects of treatment on the centroids of treated dunes compared to their control pairs only (i.e., without mobile dunes), with years as replicates. Multivariate methods were redundant for examining rodents (since there are effectively only two species), but we included ordination plots for consistency in comparison with other taxa.

To visualize the multivariate changes over time, we plotted the RDA ordination as a function of year and used the mean scores for each dune type on the first RDA axis (RDA1). This is conceptually related to the Principal Response Curves approach (PRC) [90], but here our analysis is conducted without standardizing to a baseline. We did this because we wanted to be able to see the fluctuation in all dune types in response to environmental stochasticity.

Since the data were unbalanced, we could not run a standard PRC permutation test, so to test for temporal effects on treated dunes, we used Permutational Multivariate Analysis of Variance (Permanova). These tests were conducted for each year separately, with the treated dune type set as the reference dunes with which to compare both the control pair and the reference (mobile) dunes simultaneously.

#### 2.4.3. Indicator Species

We identified Indicator Species (IS) for each taxon in two stages. First, we identified species that are highly affiliated to a single dune state (mobile, semi-fixed or fixed dunes) using the *IndVal* method. This method was proposed by Dufrêne and Legendre [91] to identify the ecological preference of a given species among a set of alternative site groups. An Indicator Value (IV) is a measure of association ranging from 0% (individuals are spread equally between the dune types) to 100% (all individuals of a species are observed in all dunes of only one dune state i.e., the indicated dune state). IVs were calculated for each species with the *indicspecies* package in R [92]. To avoid rare species that were unlikely to be captured or useful for monitoring, we only included common species with total abundance of 20% or more of the average species abundance within each taxon. Note this was in addition to the species that were already removed from the whole study as described earlier in Section 2.4. IVs were tested against a null distribution using a permutation test [93]. IVs that were significantly different from the permuted random distribution (*p* < 0.01) were considered as strong candidates to be indicator species.

The value of IV is the degree of affinity of each species to a given dune state, but it does not indicate the vegetation cover that the species is associated with. To understand this association, we calculated a second measure which reflects the percentage of vegetation in which a species is most commonly active; this is effectively PPC weighted by abundance *(wPPC)* such that:

$$wPPC\_i = \frac{\sum\_{1}^{i} \left(n\_{ij} \times PPC\_i\right)}{\sum\_{i}^{j} \left(n\_{ij}\right)}\tag{2}$$

where *n* is the total abundance over time of species *i* in dune *j*, and *PPCj* is the average perennial plant cover in dune *j.* Thus, a highly obligate mobile dune species would have a high IV and a low *wPPC*, while a highly obligate fixed dune species would have a high IV and a high *wPPC*.

Theoretically, *wPPC* can range from 0% to 100%, but the dunes in the Nizzanim TER varied from 10% to 36% cover, so these are the minimum and maximum possible values in this study, respectively. Equation (2) forms the preliminary steps for calculating the Species Sandiness Index (SSI) in Rubinstein et al. [77], but the range of SSI is dependent on the range of plant cover included in the study and is designed to provide a scale where sandiness in mobile dunes has higher values than in dunes with high vegetation cover. The value of *wPPC* for each species is comparable across other habitats, and the positive direction reflecting increasing PPC is more informative in the context of this study.

For each taxon, we objectively selected the two IS with the lowest and highest *wPPC*, to reflect the species most affiliated to mobile and fixed dune states respectively. This stage was done using data from control dunes only. These IS are expected to be the most specialized, adapted or obligate species for the respective dune states. That said, on the fixed dunes they could also be generalist species that are at the edge of their natural (e.g., inland) range and unable to inhabit the more dynamic dune states.

We then tested whether effect of the treatments on the eight selected IS was according to the expectations, i.e., away from the control dunes towards the reference dunes. According to the H1 hypothesis, we should expect that mobile IS would increase in abundance in the treated dunes, while fixed IS would decrease.

The abundance of each IS was tested using GLMM-PQL (reptiles and beetles) and LME (rodents and annuals) models similar to those used for assemblage abundances [26]. We assessed the temporal trends of each IS using a permuted Empirical Cumulative Distribution Function (*ecdf*) in the *R stats* package [84]. This simple permutation approach generates a *p*-value for the null hypothesis of zero difference in abundance of a species between two habitats. It allows us to consider unbalanced data with zero-inflated distributions that are not testable under ANOVA. We conducted each test on pairwise treatment and control, or treatment and reference dunes for each year separately.

#### **3. Results**

Overall, across a 12-year study in 21 dunes, we analyzed data from 19 reptile species (*n* = 1265 records), five rodent species (*n* = 824), 48 beetle morphospecies (*n* = 6207), and 66 annual plant species, after removal of rare species.

#### *3.1. Abundance & Richness*

We tested assemblage-level abundance and species richness to determine whether these parameters were affected by vegetation removal. We considered any result with *p* < 0.05 as significant Full results and *p*-values are given in Supplementary Tables S2 and S3. We plotted the results as boxplots where the upper and lower edges of the box correspond to the 25% and 75% quartiles (Figures 2 and 3). The upper whiskers extend from the edge of the box to the highest value, excluding outliers. Data beyond 1.5 of the 25–75% inter-quartile range are outliers and plotted as points. Note that for many of the rodent samples, either one or two species was found. Therefore, for some dune types, the boxplot appears as a single line with no variation. On its own, boxplots for the rodent data would normally not be the most appropriate way to present these results. However, we felt it pertinent to maintain the same graphics across taxa to make them more easily comparable.

**Figure 2.** Boxplots of average abundance for each taxon (see taxonomic icons by row) for (**a**–**d**) control dunes, (**e**–**h**) semi-fixed dune treatment, (**i**–**l**) fixed dune treatment, and (**m**–**o**) off-road sport vehicle disturbance. Each treated dune type was tested against its untreated pair (control) and mobile (reference) dunes. Disturbed dunes were not monitored for annuals, and untreated semi-fixed dunes were used as the control. Significant differences were tested with generalized mixed models for beetles and reptiles (Poisson), and linear mixed models for rodents and annuals. *p* < 0.05 for treated dunes compared to \* control and + reference dunes.

**Figure 3.** Boxplots of species richness for each taxon (see taxonomic icons for each row for (**a**–**d**) control dunes, (**e**–**h**) semi-fixed dune treatment, (**i**–**l**) fixed dune treatment, and (**m**–**o**) off-road sport vehicle disturbance. Each treated dune type was tested against its untreated pair and reference (mobile) dunes. Disturbed dunes were not monitored for annuals and un-treated semi-fixed dunes were used as the control. Significant differences were tested with generalized mixed models for beetles and reptiles (Poisson), and linear mixed models for rodents and annuals. *p* < 0.05 for treated dunes compared to \* control and + reference dunes, respectively.

We first tested un-treated fixed control and semi-fixed control dunes against the reference mobile dunes, to set up the expectations for treatment effects. There was no difference between control and reference dunes for any of the faunal taxa (Figure 2a–c). This meant we rejected *H1* that abundance of animals increases with increasing primary productivity. For annual plants, increasing perennial cover led to increased annual cover as expected (Figure 2d).

These findings in control dunes established the expectations for responses to treatment (*H3*): For annuals, we expected treatment to result in a decrease in cover, while for animal taxa, no difference was found between control and reference dunes, so no difference would be expected in response to treatment.

For the semi-fixed treatment, there was no significant difference in abundance between treated and untreated semi-fixed dunes or mobile dunes for all three faunal taxa (as expected based on the expectations set up in the previous statement (Figure 2b,e–g). In the fixed dune treatment, rodents maintained similar (not significantly different) abundance to un-treated control and mobile dunes, as expected (Figure 2i), while reptile abundance was significantly lower in treated fixed dunes than their control pairs (Figure 2j). Conversely, beetle abundance was significantly higher in treated fixed dunes than both control and reference dunes (Figure 2k). Thus, for reptiles and beetles, treatment did have an effect, despite there being no expectation for any treatment effect.

Disturbed dunes exposed to off-road vehicle disturbance were also expected to be more similar to mobile and less similar to untreated semi-fixed control dunes. There was no effect of vehicular disturbance on rodent abundance (Figure 2m). However, reptile abundance was significantly lower on these dunes compared to both mobile and semi-fixed control dunes (Figure 2n). In contrast, beetle abundance was significantly higher in disturbed dunes than in mobile and semi-fixed control dunes (Figure 2o).

Cover of annuals increased with increasing dune fixation so, treatment was expected to reduce % cover of annuals, making the treated dunes more similar to mobile dunes. However, on semi-fixed dunes, annual % cover did not significantly decrease in response to treatment, with % cover remaining similar to semi-fixed controls and significantly higher than in mobile dunes (Figure 2h). In treated fixed dunes, the % cover was significantly higher than both mobile and fixed control dunes, which meant the direction of response was opposite to expectations (Figure 2l).

As with abundance, we tested the control dunes for differences in species richness (Figure 3a–d). For rodents and beetles there was a gradient of declining species richness with increasing perennial cover (fixation state) in control dunes (Figure 3a,c), which opposes the productivity-diversity theory (*H2*). Thus, the expectations for treatment would be an increase in richness following removal, if treatment was successfully replicating mobile dune conditions and therefore biodiversity structure. Note that for beetles, the decline in richness was only significant between mobile and fixed dunes, so the semi-fixed treatment was expected to remain similar to both control and reference. Conversely, reptile and annual plant richness increased with increasing dune fixation (Figure 3b,d), supporting the theory that productivity increases diversity. The expectations for these latter taxa were thus that species richness would decrease in response to treatment.

In semi-fixed dunes, no response of richness was seen in reptiles (Figure 3f) and annuals (Figure 3h) which meant an absence of expected effect, while beetle richness remained unchanged in response to treatment as expected (Figure 3g). In fixed dunes, rodent and annual richness showed no difference between treated and untreated fixed dunes, both remaining significantly different from mobile dunes (Figure 3i,l) and contradicting the expectations for a decrease (rodents) and increase (annuals) in response to treatment. The apparent reduction in reptile richness in treated dunes became similar to (i.e., not significantly different from) to mobile dunes (Figure 3j), but the decline was not large enough to become significantly different from the untreated control dunes (*p* = 0.10) so the expected result was only half met. Beetles showed a small increase in treated fixed dunes towards mobile dunes compared to untreated fixed dunes, and thus beetle species richness did not differ significantly from the reference or the control (Figure 3k).

Lastly, dunes disturbed by ORV had significantly higher rodent species richness compared to untreated semi-fixed dunes, similar to mobile dunes (Figure 3m). Reptiles were found to have significantly lower species richness on disturbed dunes compared to semi-fixed control dunes, but were not significantly lower than mobile dunes (Figure 3n). Beetle richness in treated semi-fixed dunes remained similar to richness in both untreated semi-fixed dunes and mobile dunes (Figure 3o).

#### *3.2. Composition*

To examine the effect of treatments on composition we used constrained multivariate ordination. Firstly, we conducted a redundancy analysis (RDA) of all dune states together, in order to compare the relative composition of all dunes in the study, and to visualize the relative magnitude of the effect across press and pulse treatment types (Figure 4a–d). We then tested the treatment effect in each dune state separately, in order to meet the assumptions for the statistical tests *(*Figures 5–7*)*. For all taxa, composition was a good indicator of dune state as expected under hypothesis *H2*, as seen by the strong clustering and gradient from fixed control to mobile dunes in all taxa along the first RDA axis (RDA1). Explained variation was much higher in RDA1 than RDA2 for all taxa (Figure 4).

Rodent assemblages in Nizzanim coastal dunes were heavily dominated (98%) by two gerbil species *Gerbillus pyramidum* and *G. andersoni allenbyi.* Therefore, ordination formed an arc effect along the first axis (Figure 4a) as the ratio between the two species shifted from only *G. allenbyi* in fixed dunes, to only *G. pyramidum* in disturbed dunes. The disturbed dune centroid was located beyond (further left) the mobile dune's along RDA1 (yellow arrow, Figure 4a), since both species are present in the latter.

**Figure 4.** Overall summary of redundancy analysis (RDA) analyses with all dune types included for (**a**) rodents; (**b**) reptiles (**c**) beetles and (**d**) annual plants. Arrows depict direction of treatment effect along RDA1. Significance testing was performed in each dune type separately (see Figures 5–7).

**Figure 5.** Redundancy analysis for (**a**) rodents, (**b**) reptiles, (**c**) beetles and (**d**) annual plants in treated semi-fixed dunes (pale blue triangle) compared to un-treated semi-fixed (control) dunes (dark blue square) and Mobile dunes (orange square) which are the treatment reference. Distance bars along the 1st RDA axis show the treatment effect (blue arrow) as a percentage of total distance between control and reference (black dashed). Significant differences were tested using pairwise Permanova. \* and + represent *p* < 0.05 for treated dunes compared to control and reference dunes, respectively.

**Figure 6.** Redundancy analysis for (**a**) rodents, (**b**) reptiles, (**c**) beetles and (**d**) annual plants in treated fixed dunes compared to untreated fixed (control) and mobile dunes (reference). Distance bars along the 1st RDA axis show the treatment effect (green arrow) as a percentage of total distance between control and reference (black dashed). Significant differences were tested using pairwise Permanova. \* and + represent *p* < 0.05 for treated dunes compared to control and reference dunes, respectively.

Both semi-fixed and fixed dunes showed minor shifts towards mobile dunes (the expected direction for our restoration goals), as a result of treatment for faunal taxa (blue and green arrows in Figure 4a–c), which supported hypothesis *H4*. A small and insignificant shift in the opposite direction was seen for annuals (green arrow; Figure 4d). Annuals in both semi-fixed control and treated dunes were similar to mobile dune in RDA1 (Figure 4d). ORV disturbed dune composition was virtually identical to mobile dune composition for beetles (Figure 4c), and was similar for reptiles along RDA1, but was distant along RDA2 (yellow arrow; Figure 4b). Overall, press treatment (ORV) appeared to be more similar to reference dunes, compared to pulse treatments, in accordance with hypothesis *H5.*

Changes in response to treatment were tested for each dune type separately, comparing each treatment to the reference (mobile) dunes, and their control dunes for semi-fixed (Figure 5), fixed (Figure 6) and disturbed dunes (Figure 7). Pairwise Permanova results are given in Table S4. Rodent composition shifted 30% along RDA1 towards the reference as a response to treatment (Figure 5a). Their composition was significantly different between treated and control dunes for semi-fixed dunes. In other words, the ratio between the two rodent species changed significantly as a result of the semi-fixed treatment. Reptiles and beetles shifted along the 1st RDA by 17% and 12% respectively, towards the reference centroid in response to treatment in semi-fixed dunes (Figure 5b,c). The annual plants' centroid shifted away from the reference centroid by 9% along RDA1 (Figure 5d). The centroid shifts detected for these three taxa were not significant in pairwise Permanovas. The 2nd RDA axis for all taxa explained a very low percentage of the variance (0.7–3.3%) compared to the first axis (11.3–49.5%).

**Figure 7.** Redundancy analysis for (**a**) rodents, (**b**) reptiles and (**c**) beetles in disturbed dunes (off-road vehicle disturbance) compared to semi-fixed untreated (control) and mobile dunes (reference). Distance bars along the 1st RDA axis show the treatment effect (yellow arrow) as a percentage of total distance between control and reference (black dashed); the % is calculated under the scenario that pre-disturbance state was semi-fixed. Significant differences were tested using pairwise Permanova. \* and + represent *p* < 0.05 for treated dunes compared to control and reference dunes, respectively.

Rodent composition on control fixed dunes consisted of just one species (*G. a. allenbyi)* and the second species did not reappear on treated dunes within the timeframe of the study (Figure 6a). Fixed dune treatment produced a minor shift in the intended direction towards the reference for reptiles (2%) and beetles (6%) (Figure 6b,c), while a greater shift (17%) away from mobile dunes was detected for annuals (Figure 6d). Treated dunes' centroids were not significantly different from controls in any taxa in pairwise Permanovas (see Table S4*).* The composition of rodents in disturbed dunes shifted from two species found on both mobile and untreated semi-fixed dunes, to one psammophilic species forming an arch-effect in ordination (Figure 7a). Species composition was significantly different to both mobile and un-treated semi-fixed dunes for rodents and reptiles (Figure 7a,b). The centroids for both of these taxa on disturbed dunes were beyond (left of) the centroids of mobile dunes.

The pre-disturbance state for ORV dunes could have been mobile or semi fixed. Under the latter conditions, one could infer a 155% and 261% shift from control to reference along the 1st RDA axis for rodents and reptiles respectively. For beetles, the composition shift along the 1st RDA axis was 88% of the distance between control and reference dunes and was only significantly different from semi-fixed controls in Permanova (Figure 7c). Note that for reptiles, significant differences between disturbed and mobile dunes were found along the primary RDA axis (Figure 7b), while these differences were reflected in the second RDA axis in the overall ordination shown in Figure 4b. This is because the relative difference in composition between disturbed and mobile dunes is much smaller than the difference between mobile and untreated semi-fixed or fixed dune states, irrespective of treatment type.

#### *3.3. Temporal Trends in Composition*

To visualize the multivariate changes over time, we plotted the mean scores for each dune type on the first RDA axis (RDA1) from Figures 5–7 as the *y*-axis, as a function of year along the *x*-axis (Figure 8). Rodents on treated semi-fixed dunes were significantly different from their untreated pairs in 2014 (Figure 8a). Beetle composition in treated semi-fixed dunes changed incrementally over time towards the mobile dune centroid, and was eventually significantly different from untreated semi-fixed dunes in 2016 (Figure 8g). A small shift was seen for reptiles and rodents from untreated to treated semi-fixed dunes in 2016, but it was not significant (Figure 8a,d). These findings could be seen as partial support for hypothesis *H4a.*

**Figure 8.** Temporal Redundancy Analysis (RDA) curves for each treatment type (see column titles) for (**a**–**c**) rodents, (**d**–**f**) reptiles, (**g**–**i**) beetles, and (**j**–**k**) annual plants. Scores on the first axis (RDA1) in constrained redundancy analysis with year as a covariable are plotted over time. Data points are the scores for centroids of each dune type in a given year for each taxon. *p* < 0.05 for all significant results; \* Treated semi-fixed dunes were significantly different to untreated semi-fixed control dunes in specified years. \*\* Disturbed dunes were significantly different to un-treated semi-fixed control dunes in all years for all taxa. + Disturbed dunes were significantly different to reference Mobile dunes in specified years. ++ Treated dunes were significantly different to the reference in all years for all taxa.

All faunal taxa had significantly different composition in treated semi-fixed dunes and reference (mobile) dunes across all years (Figure 8a,d,g). The same was true for treated fixed dunes (Figure 8b,e,h). In summary, none of the treated dunes fully reached the reference state in any year for faunal taxa.

For annuals, the semi-fixed treatment had no effect on RDA1 (Figure 8j). Fixed dune treatment appeared to have no effect on RDA1 for any taxa (Figure 8b,e,h,k) supporting the resistance hypothesis *H4c*. Conversely, disturbed dunes were significantly different from untreated semi-fixed dunes across all years (Figure 8c,f,i) suggesting turnover (*H4a*). In addition, disturbed dunes were significantly different from mobile dunes in all years for reptiles (Figure 8f), in 2012 and 2015 for rodents (Figure 8c) and in 2015 for beetles (Figure 8i). In all other years, rodents and beetles were similar to (i.e., not significantly different from) mobile dunes, in RDA1.

Lastly, we conducted a final RDA ordination for each taxon, this time including all dune types across all years (Figure 9). This allowed us to visualize the relative differences in composition across all dune sites within the study area. Irrespective of treatment, semi-fixed dune composition was more similar to fixed dunes than mobile dunes for rodents and beetles (Figure 9a,c), and relatively similar to mobile dunes compared to fixed dunes, for annuals and reptiles (Figure 9b,d). While disturbed dunes were significantly different to mobile dunes in many samples (Figure 8c,f,i), these dunes were still the most similar to mobile dunes relative to other treatments for all faunal taxa (Figure 9a–c).

**Figure 9.** Temporal Redundancy Analysis (RDA) curves (**a**) rodents, (**b**) reptiles, (**c**) beetles and (**d**) annual plants, for all dune types together showing relative contribution to the first axis from 2006–2016, with year as a covariable. Data points are the scores for centroids of each dune type in a given year.

Figure 10 summarizes the overall findings for abundance, richness and community level responses across all treatment described above.

**Figure 10.** A summary of the responses for three dependent parameters to the three independent treatments, across four different taxa. Note some parameters showed no significant change, which could still be the 'desired' response if there was no difference between control and reference dunes. The comparative strengths of the responses are not depicted.

#### *3.4. Indicator Species*

Twenty-eight species were identified as suitable candidates to be indicator species (IS) based on the *IndVal* analysis, including two rodents, seventeen annuals, nine beetles, and four reptiles (Table 1). The Indicator Value or *IV* (a measure of habitat affinity) ranged from 61% to 98% across all taxa. Weighted PPC (*wPPC)* ranged from 11% for the ground beetle *Scarites striatus* to 34% for the annual *Bromus rigidus.* Only one reptile (*Stenodactylus sthenodactylus)* and five annual species were identified as potential IS for semi-fixed dunes. Most annuals (*n* = 11) were fixed dune affiliates with high SSI. Only one annual plant species (*Cutandia memphitica*) was identified as a potential IS for mobile dunes, although it had a relatively low IV of 64% and relatively high *wPPC* of 18%.


**Table 1.** Candidate species identified as potential indicators, using the *IndVal* function, the dune type they are affiliated to, the Indicator Value (IV) and the weighted Perennial Percentage Cover (wPPC) for each species. Species highlighted with the lowest and highest wPPC were chosen for further analysis as \* Mobile dune Indicator Species (MIS) and \*\* Fixed dune Indicator Species (FIS) respectively.

The species with the highest and lowest *wPPC* were selected as IS for further analysis, such that every taxonomic group had one Mobile dune Indicator Species (MIS) and one Fixed dune Indicator Species (FIS) (species in bold in Table 1). The FIS for rodents is actually an indicator for semi-fixed dunes, as identified by the IV analysis, but we refer to it as an FIS for convenience to compare with other taxa since no true FIS was identified. Full model outputs for both MIS and FIS are given in Supplementary material Tables S5 and S6.

As expected, the abundances of MIS were much lower (or absent) on untreated semi-fixed dunes compared to the mobile dunes (which is what defines them as good IS with high IV). No change in abundance was visible for the gerbil *G. pyramidum* (Figure 11a) in treated semi-fixed dunes. A small increase in abundance for the lizard *A. scutellatus* was also observed in treated semi-fixed dunes, becoming more similar to levels found in mobile dunes, but the increase was not great enough to become significantly different from untreated semi-fixed dunes *(*Figure 11b). Meanwhile a significant increase in abundance was observed for *S. striatus* in treated compared to untreated semi-fixed dunes, although it did not reach reference levels (Figure 11c). The annual MIS, *C. memphitica* had a large range of cover across samples, and cover in treated semi-fixed dunes was not significantly different to mobile or untreated semi-fixed dunes (Figure 11d).

**Figure 11.** Boxplots of (**a**–**d**) Mobile dune Indicator Species (MIS) and (**e**–**h**) Fixed dune Indicator Speices (FIS) for each taxon (see icons by row), showing the 25% and 75% percentiles of the measure of abundance (Y axis units e.g., rank abundance/ % cover) on each dune type for rodents, reptiles, beetles and annuals. Dune types were tested for significant differences between treated dunes, their control pair and reference mobile dunes, using mixed effect models. Disturbed dunes were tested against semi-fixed control dune and mobile dunes. Significant results (*p* < 0.05) are shown for treated dunes compared to \*control and <sup>+</sup>reference dunes (dune type abbreviations refer to the Nizzanim LTER coding system for continuity).

While not significantly different from control, we calculated that both the rodent and reptile MIS abundance increased by 15% of the difference between control and reference levels overall (Table 2i a and b). This can be considered as 15% of the desired response, if the treated dunes were expected to become identical to mobile dunes (*H6*). For beetles the increase was 27%, and for annuals there was overcompensation in the MIS abundance, with a 114% increase compared to reference (Table 2).


**Table 2.** Effectiveness of treated dunes on abundance of Mobile and Fixed-dune Indicator Species (IS). Calculated as the change (either increase or decrease) in abundance (faunal taxa) or % cover (annuals) from control to treated dune, as a percentage of the difference in abundance between control and reference dunes for each species.

The abundances of the MIS *G. pyramidum* and *A. scutellatus* on disturbed dunes were also significantly higher than their abundances in semi-fixed control dunes and were similar to the abundances found on mobile dunes (Figure 11a,b), increasing by 126% and 68% of the expected abundance respectively (Table 2). The abundance of *S. striatus* was significantly higher on disturbed dunes compared to semi-fixed control dunes (42% increase), but it was still significantly lower than on mobile dunes (Figure 11c). None of the MIS responded to fixed dune treatments rejecting H6; they remained absent on treated fixed dunes (Figure 11a–d).

FIS was not found to show any significant response to fixed dune treatment in any taxa (Figure 11e–h) leading us to reject *H6*. Although the abundance of *G. a. allenbyi* significantly decreased in treated semi-fixed dunes compared to their untreated pair, the response was not enough to become similar to mobile dunes (Figure 11e). Lastly, disturbed dunes had a significantly lower abundance of this rodent than both mobile and untreated semi-fixed dunes, while the beetle FIS was only found to be significantly lower than untreated semi-fixed dunes, being almost entirely absent from mobile dunes (Figure 11g).

Mean abundance of rodent and beetle FIS declined by 30% and 33% of the desired reduction respectively (Table 2). Annual FIS % cover declined by 21%, and there was no change in fixed dunes for reptile FIS. In addition, mean abundance decreased for all FIS in semi-fixed treated dunes by 30%, 36%, 30% and 42% in rodents, reptiles, beetles and annuals respectively, towards the reference.

We plotted the abundance of each species (mean for each dune type) over time to show populations fluctuations (Figure 12). These abundances were tested yearly in a pairwise permuted distribution test. The significance of these tests should be interpreted with caution due to the small sample sizes and multiple testing. Nevertheless, some interesting trends were apparent across the different taxa.

All MIS abundances fluctuated stochastically on the mobile dunes over time (Figure 12a–d). They were all observed with relatively low to zero abundance on semi-fixed control dunes, with the exception of the lizard *A. scutellatus*, with intermediate abundance in semi-fixed dunes. All MIS were effectively absent from fixed control dunes. In semi-fixed treated dunes, the abundances of both the rodent *G. pyramidum* and the beetle *S. striatus* were similar to untreated semi-fixed dunes (i.e., low abundance) and were significantly different from AC, until 2016. In the final year of monitoring, the abundances of both these species increased and they became similar (not significantly different from) mobile dunes and significantly different from untreated semi-fixed dunes (Figure 12a,c). The same was partially true of the annual *C. memphitica*, which was similar to untreated semi-fixed dunes until the final year 2016, but was not consistent across years in comparison to mobile dunes (Figure 12d). The lizard *A. scutellatus* showed a small but similar trend in the final year of monitoring in treated semi-fixed dunes, but it remained significantly different from mobile dunes for all years, and similar to untreated semi-fixed dunes in all years except 2013 (Figure 12b).

**Figure 12.** Temporal changes in abundance of (**a**–**d**) Mobile Indicator Species (MIS) and (**e**–**h**) Fixed Indicator Species (FIS) on different dune types for rodents, reptiles, beetles, annuals (see icons by row). Dune types were tested for significant differences between treated dunes, their un-treated control pair and the reference mobile dunes each year. \* and + represent p<0.05 for treated dunes compared to control and reference dunes, respectively. Not all significant results are shown given for intelligibility. See Supplementary Tables S7 and S8 for full results.

In the fixed dunes, none of the MIS were found to have differences in abundance between treated and un-treated fixed dunes (i.e., they remained absent or in very low abundance), and all the faunal MIS remained significantly different from mobile dunes in all years (Figure 12a–c).

In the disturbed dunes, *G. pyramidum* had similar abundances to mobile dunes in all years except 2014 (Figure 12a). The abundance of *S. striatus* was significantly lower in disturbed than in mobile dunes in all years except 2014, although there did appear to be a trend of increasing abundance in disturbed dunes over time (Figure 12c). In addition, the latter species was observed in disturbed dunes in all monitored years, while it remained absent from un-treated semi-fixed dunes throughout the study. In fact, both *G. pyramidum* and *S. striatus* had significantly higher abundances in disturbed dunes, compared to untreated semi-fixed dunes in all years (Table S7). Abundance of *A. scutellatus* on disturbed dunes was most similar to mobile dunes in all years except 2015, and was significantly different to untreated semi-fixed dunes in 2013 and 2014 (Figure 12b).

Next, we considered the effect of each treatment on FIS (Figure 12e–h). In semi-fixed dunes, the rodent *G. a. allenbyi* initially had similar abundances in treated and untreated semi-fixed dunes, which were significantly higher than in mobile dunes, but later in the study, the abundance of this species declined and became more similar to mobile dunes and less similar to untreated semi-fixed dunes (Figure 12e). The abundance of the lizard *A. schreiberi* was similar across mobile, treated and untreated semi-fixed dunes (i.e., low to medium abundance), except in 2012 where it was absent in mobile dunes, but remained present in both treated and untreated semi-fixed dunes (Figure 12f). The abundance of ground beetle *G. sharonae* was initially similar (low abundance) in mobile, treated and untreated semi-fixed dunes. Then in later years, abundances in treated semi-fixed dunes became different to mobile dunes, but remained similar to untreated semi-fixed dunes (both increasing slightly) in 2015

and 2016 (Figure 12g). The annual plant *B. rigidus* had generally very low cover in all dune types (Figure 12h). The cover in treated semi-fixed dunes was similar to mobile dunes in all years (i.e., absent) and similar to untreated semi-fixed dunes in all years except 2014, when some presence was observed on untreated semi-fixed dunes, while remaining absent from treated semi-fixed dunes.

In the treated fixed dunes, abundance of *G. a. allenbyi* was initially different to mobile dunes but became similar in the later years of the study (Figure 12e). However, abundance of this species in DC also showed declines until 2009, and abundance in treated and untreated fixed dunes was similar in all years except 2005 and 2008. Unfortunately, data is lacking for the period between 2009–2016 on treated fixed dunes, and trends in the interim years are unknown. The other FIS, *A. schreiberi*, *G. sharonae*, and *B. rigidus* had similar abundances on treated and untreated fixed dunes throughout all years, despite stochasticity (Figure 12f–h). A marked increase in cover was observed in the last year of the study for *B. rigidus*, but this was detected in both treated and untreated fixed dunes, so overall, no difference was found. Both *G. sharonae* and *B. rigidus* were absent from mobile dunes in all years. All three FIS were significantly different from mobile dunes in all years (Figure 12e–g), excluding years where *B. rigidus* was absent in both treated and untreated fixed dunes and mobile dunes (Figure 12h).

In the disturbed dunes, all FIS were effectively absent, with the exception of *G. a. allenbyi*, whose abundance was low to medium on both disturbed and mobile dunes in most years (Figure 12e). This species was not detected in disturbed dunes in 2013, while it was still observed in mobile dunes. The abundance of this species in untreated semi-fixed dunes was significantly higher than in disturbed dunes for all years. *A. schreiberi* was present occasionally in low abundances on mobile and untreated semi-fixed dunes over the course of the study, but remained absent in disturbed dunes in all years monitored (Figure 12f). *G. sharonae* was effectively absent on both disturbed and mobile dunes in all years, and only occasionally appeared in low abundances on untreated semi-fixed dunes. Thus, abundance of this FIS in disturbed dunes was found to be similar to both mobile and untreated semi-fixed dunes in all years (Figure 12g).

#### **4. Discussion**

Coastal dune habitats have declined globally over the last several decades due to urbanization. Within remaining dune systems, substantial dune fixation has resulted in further losses of mobile dunes with negative impacts on their associated species [15,19,25,27,36,49,78,94]. Some studies suggest vegetation removal can initially promote habitat heterogeneity, and increase availability of suitable habitats for psammophile, xeric and endemic mobile dune species, but longer term responses are generally unknown [15,24,28,40,54]. We investigated the temporal trends of four taxonomic groups (including fauna) to determine the effect of vegetation removal on dune assemblages over a 12-year period. Detecting community level responses to restoration can be difficult. Indeed, two studies in Californian coastal dunes were unable to show significant assemblage level responses for terrestrial arthropods [32] or plants [24] despite intensive vegetation removal. Our study found only marginal restoration effects in most cases.

This large-scale project across multiple taxa and multiple treatments meant that our samples sizes were relatively small, and some caution is needed to interpret the results. Nevertheless, some inferred conclusions can be made. Comparison of control dunes only supported *H1* for annual plants (% cover increases with perennial cover) and *H2* for annuals and reptiles (richness increases with perennial cover). Composition level parameters showed only small shifts away from the control dune, suggesting high levels of resistance at the community level. Although many of these changes were not statistically significant, they were usually in the direction desired by the manipulation, towards mobile dune status. Semi-fixed dune responses were more pronounced than fixed dune responses. The ORV-disturbed dunes shared more aspects with mobile dunes for faunal taxa compared to both semi-fixed and fixed treated dunes, while annuals were effectively absent. Most mobile dune indicator species also showed some trends in the intended direction in response to most treatment types suggesting that treatment did increase habitat availability for these species, while not removing suitable habitat for fixed dune

species. Our results suggest that more pronounced responses might occur given a longer time frame, and that press disturbance is more effective than pulse disturbance for restoring mobile dune habitat.

Overall, abundance and richness measures were not appropriate indicators for treatment success (since there were few differences between control dune states). Richness did correlate with the gradients of increasing PPC across control dune dunes, but the effects of treatment detected for this parameter were marginal. Both abundance and richness are important to monitor in restoration projects, but they are more appropriate for assessing whether any unexpected artefacts resulted from the intervention. Composition is known to a better parameter for detecting changes [26,28]. The indicator species that were selected also were highly suitable for detecting changes. Both composition and indicator species showed a temporal trend in treated semi-fixed dunes, away from the control dunes towards the reference dunes over time; this trend was more consistent in the faunal taxa. Over a period of 12 years, we found that the response at the assemblage level was between −9% to 30% in semi-fixed treatment, compared to just −17% to 6 % in fixed dunes.

#### *4.1. E*ff*ect of Treatment on Semi-Fixed Dunes*

Establishing dune dynamics on the wind facing slope can lead to a rejuvenation of the landscape through reactivation of sand, leading to burial of any newly developing vegetation [31]. This intervention enables pioneer species found in mobile dunes to recolonize [4]. Rodents responded well to treatment in semi-fixed dunes with a significant shift in composition detected, achieving 30% of the desired Euclidean ordination distance towards reference, as expected by Hypothesis *H4*. In addition, though the change was not significant for reptiles or beetles, the response to treatment was consistently in the desired direction towards mobile-dune composition. For annuals the shift was minor, but away from the desired composition, becoming less similar to mobile dunes, and *H4* was rejected.

At the temporal scale, rodents showed significant differences in some years. The change in reptile and beetle composition was more pronounced in the final years of the study, and became significant for beetles in the last year. Similar temporal responses were also detected at the indicator species (IS) level, particularly for "mobile-dune indicator species" (MIS). These results indicate some support for Hypothesis *H4a*, where the composition turnover is due to recolonisation by MIS.

The IS were selected objectively based first on their affinity to fixed and mobile dunes, and secondly as the species with the lowest and highest weighted PPC. Three MIS (all except the reptile), as well as the rodent "fixed-dune indicator species" (FIS) appeared to respond desirably in the direction intended for restoration. Given that the beetle *S. striatus* was almost entirely absent from control dunes and increased over time in treated dunes, we can infer that its increased abundance was due to colonization rather than competitive release. The small increase detected for rodents and reptiles was not significant when looking at averages across years, but the gerbil *G. pyramidum* did respond favorably over time, showing significant increases in the final year of the study. This increase may be due to the overall decline in abundance of the subordinate *G. a. allenbyi* populations, since these two species compete sympatrically at the patch scale [95–97]. Further investigation is needed, since the trends of the two species were not consistent over time. The annual MIS *C. memphitica* also had relatively high cover in some years on treated dunes, while remaining much less prevalent in control dunes. This was most likely a response to competitive release following the removal of perennial vegetation. Over time, this species was erratic in both mobile and treated dunes and the temporal trends were not synchronized, making it a poor IS, as expected due to the relatively low IV.

Given similarities between control and reference dunes in terms of reptile and beetle species richness, there was no expectation for treated dunes to elicit any response. However, control and reference dunes were different for annual plants, so for this taxon, a lack of response can be seen as lack of desired treatment effect. Only rodents increased in richness as expected under Hypothesis *H3b* (from one to two species) and *H3a* was rejected for all taxa in semi-fixed treated dunes. These results highlight the importance of examining composition rather than univariate measures when studying restoration effects on biodiversity.

The composition shift in semi-fixed treated dunes may continue towards a mobile state with more species responding positively, but without a press disturbance, we may observe a return to the control state. In the Netherlands, aeolian activity increased after an intervention, up to a climax in the third year, and then started to decline because of vegetation development on bare spots [4]. Arens et al. [40] suggested that follow-up removal of roots is needed for three to five years and it is not clear whether the dunes will remain mobile in the long term.

#### *4.2. E*ff*ect of Fixed-Dune Treatment*

A recent review of coastal restoration projects revealed that fixed dunes are the most frequently restored dune type (54%), followed by semi-mobile dunes (30%), mobile dunes (27%) and dune slacks (16%; [2]. However, the most frequently used restoration mechanism was revegetation (42%), while destabilization or vegetation removal only accounted for less than 5% of restoration projects. Removal methods varied among removal experiments.

The extensive and intrusive removal of perennial vegetation in a grid formation on fixed dunes had virtually no effect on most of the parameters monitored, rejecting hypothesis *H3b*, *H4a*, and *H6 for* this treatment type. The only consistent effect was a reduction in the average abundance of all taxa (*H3a*); richness and composition were not affected. That said, at the temporal scale, some weak responses of reptiles were observed both for the IS and for composition in the first two years after treatment, suggesting some support for resilience under Hypothesis *H4b*. All other taxa showed very limited response to fixed dune treatment, indicating strong resistance to change (supporting Hypothesis *H4c*). Detachment from mobile dune assemblages may explain this lack of response. In addition, the dunes from which vegetation was removed are surrounded with dense interdune vegetation, including tall sycamore trees, which reduce wind speed, preventing sand migration and potential recolonization from neighboring areas. Thus, sand deposition in the treated dunes was minimal, and sand enriched with nutrients and organic matter remained present after the vegetation removal. Throughout Europe, several coastal dune management schemes involving vegetation removal achieved partial success in restoring dune mobility, from a geomorphological perspective [6,15,31,40].

Evidence suggests that soil processes in fixed dunes have changed the soil irreversibly [19]. In experimental dune restorations, removal of vegetation and topsoil in the slacks and inner dunes often fail to restore aeolian activity due to the groundwater table and the usually moist conditions at the surface, which prevent saltation [4]. In Mediterranean coasts, wind-related variables may only be of minor importance compared to soil properties [98]. Encroachment and intended establishment by grasses, shrubs and woody plants have accelerated in response to high rates of nitrogen deposition [33,99–102]. Top-soil inversion has become a popular method for restoration of inland dunes in recent years, where topsoil is buried under a layer of subsoil, so that the original layers are intact but their position in the profile is changed [103,104]. Results indicate that inversion can successfully lower surface soil fertility and reduce competition, creating a nutrient-poor space for the establishment of early-successional stage species in calcareous grasslands and sand dunes [104,105]. Ploughing and topsoil inversion enhanced lichens and annual plants, lowered organic matter and increased rabbit activity compared to control plots in grass heath habitats [106,107]. Thus, topsoil removal within the grid formation applied on Nizzanim fixed dunes may have created a better response at the community level.

#### *4.3. E*ff*ect of ORV Disturbance*

Off-road vehicle (ORV) disturbance is an anthropogenic press disturbance that is increasingly threatening coastal ecosystems [108,109]. In general, recreational disturbances result in a decrease of species diversity [19,110–112]. Motorbikes, 4 × 4 vehicles and quad bikes are driven on beaches and dunes for recreation purposes [113]. These vehicles can cause more widespread damage than human trampling [109,114,115], another disturbance that is a widely reported issue for fragile coastal biodiversity [110–112,115–117]. The negative impacts ORVs cause include damage of the physical properties and stability of the substrate, destruction of vegetation, and disturbing, injuring or killing

fauna [108,114,118–121]. In California, dunes exposed to ORV disturbance showed significant decline in Coleoptera diversity and species evenness [122].

The effect of disturbance on species composition depends on whether the disturbance exacerbates or reduces environmental harshness, and the conditions that favor specialization (Attum 2006). Since recreational driving in the dunes reduces perennial shrub cover and facilitates sand remobilization, this disturbance could theoretically benefit dunes. The situation is complex because we were not able to quantify the start date, frequency, or intensity of vehicular usage. It is speculated that dunes used by ORVs are most likely to have been mobile dunes, but we do not know what trends would have occurred naturally without disturbance over this time period. However, by comparing the disturbed dunes with both mobile and control (un-treated) semi-fixed dunes, we were able to quantify the outcomes of ORV disturbance on species composition.

While disturbed dunes were significantly different to mobile dunes in many samples, these dunes were still the most similar to mobile dunes relative to other treatments for all faunal taxa. Nevertheless, there was undoubtedly a negative impact of ORV disturbance in Nizzanim reserve; firstly, annuals were almost entirely absent due to the intensity of disturbance preventing seedling establishment. Secondly, a significant reduction in total abundance and overall richness of reptile species was seen compared to both mobile and semi-fixed states.

Reptiles are known to be particularly sensitive to the presence of vehicular disturbance, affecting the relative abundance and proportional use of preferred microhabitats in coastal systems of Argentina [123]. Rocha and Bergallo [124] found that a gradual reduction in beach vegetation over a decade was accompanied by population declines of a lizard in coastal Brazil. Luckenbach & Bury [120] reported heavy impacts of ORVs on lizards, mammals and arthropods in California. We also found that the populations of the ground beetle *S. striatus* were significantly lower than on mobile dunes.

If the state of disturbed dunes prior to disturbance was more similar to semi-fixed dunes, then the ORV disturbance achieved a much greater shift in composition on these dunes (towards reference) than either of the pulse treatments achieved. Indeed, if the restoration goal was to increase the abundance of MIS, the response to ORV disturbance in the study site could be perceived as positive.

The abundance of the desert gerbil *G. pyramidum* was significantly higher in ORV disturbed dunes, compared with both semi-fixed control and mobile dunes. Thus, the overall rodent assemblage was almost entirely composed of this species. Plausibly, the 'hypermobile' conditions created are only inhabitable by highly adapted sand specialist species. Both gerbils are found in sandy habitats, but *G. pyramidum* is better adapted to exposed, harsh conditions, while *G. a. allenbyi* needs some shrub cover for protection [97,125]. Population levels of the beetle *S. striatus* were also significantly higher in ORV dunes than in semi-fixed dunes, although still somewhat lower than numbers found in mobile dunes. This species is also highly adapted to sandy desert conditions. Compared to untreated semi-fixed dunes, ORV usage resulted in an increase in highly adapted MIS, but overall alpha diversity was reduced.

In a review of coastal restoration projects [19] only two studies did *not* find a significant impact of recreation disturbance, and both were conducted in Israel's Mediterranean dunes in other coastal reserves [114,115]. Despite the differences from mobile dunes, the overall composition in ORV-disturbed dunes in Nizzanim remained the most similar to the reference mobile dunes, compared with the two press disturbances for all three faunal taxa. Thus, we found evidence supporting Hypothesis *H5* that press disturbance is more effective than a single pulse. That said, richness in these dunes was generally lower. Therefore, widespread, pervasive ORV disturbance in NDNR will likely result in an overall decrease in γ-diversity. The LTER experimental design was not balanced across dune types or years, therefore these comparisons should be taken as inferred rather than proven.

#### *4.4. Implications for Conservation Management*

Disturbance is a fundamental part of natural coastal dune processes, and encroachment by grasses, shrubs and woody plants has most likely occurred due to the removal of this process [4,33]. Geomorphologists increasingly recommend that remobilization of inner dunes should consider the dune–beach system as a whole, arguing for an integrated dynamic approach that recreates the natural disturbances such as wind erosion and sand burial [15,24,40,126]. Our results provide support to the growing body of evidence that restoration attempts should promote, not diminish, the degree of disturbance in coastal dune systems. e.g., [31,38,49].

When comparing the pulse treatments, vegetation removal of the wind-facing slopes of semi-fixed dunes appeared to be more effective than grid formation removal of fixed dunes, but since the experiments are not equivalent this can only be inferred. An initial pulse disturbance is likely needed to remove perennial vegetation, but continued press disturbance appeared to more closely match conditions found on mobile dunes, at least for faunal taxa. While the unregulated use of ORVs as a method for maintaining dune mobility is certainly not recommended, previous evidence has demonstrated the benefit of press disturbances to plant and arthropod species [15,31,49]. ORV disturbance has severe impacts on annual plant species, so other forms of press disturbance such as grazing could be more beneficial.

Grazing has been successful in increasing both plant and arthropod diversity on coastal dunes [49,127]. However, domestic grazers may affect the composition (of vegetation) differently than wild grazers [128]. A review of the effects of grazing by domestic livestock showed both positive or negative impact on diversity, often increasing diversity at the landscape level but reducing alpha diversity at the very local scale [19].

Depending on the reference species, grazing can have a positive or negative impact on richness and composition [34,127,129,130]. Ground beetles and reptiles are often negatively affected by management practices such as grazing and trampling [129,131,132]. That said, populations of *A. scutellatus* or *S. striatus* on disturbed dunes were not as low as in semi-fixed control dunes, suggesting that even recreational driving is less damaging than passive management (allowing shrub encroachment) for these species. Attum [133] showed that desert reptile specialists in the Egyptian desert did not respond favorably to heavy grazing in inland sand dunes, despite predictions that they may fare better than generalist species. Another study found arthropod assemblages were more diverse in response to disturbance by trampling, blowouts and burning in Danish dune grasslands [44]. Grazing management applied in coastal dune grasslands of France and Belgium was also able to control highly competitive plant species, and may be beneficial in maintaining mobile assemblages [134]. Grazing could be introduced either through re-wilding of the ecosystem or in the form of managed domestic grazers [49,135].

Mountain gazelles (*Gazella gazelle gazella)*, the only mammalian grazers in NDNR, are unlikely to provide an adequate grazing impact in the near future [136]. Attempts to introduce camel grazing to Nizzanim on semi-fixed dunes had little effect on the perennial vegetation [137]. Other forms of grazing or press disturbances should be investigated. Reintroduction of small livestock herds (sheep and goats) are currently being trialed with Bedouin shepherds, and should be monitored carefully before further recommendations can be made with regarding to grazing.

#### **5. Conclusions**

Conservation managers must carefully consider their priorities when implementing coastal restoration actions. Is the conservation goal a specific species that is regionally or locally endangered, or is the goal to preserve mobile dune communities, or increase alpha or gamma diversity? These decisions affect the desired actions and restoration strategy. The scale and extent of the introduced disturbances for management purposes must be carefully planned to prevent erosion and detrimental effects and should be tailored the specific goals [2].

The methods used to identify IS (*IV* and *wPPC*) appeared to be a good choice for selecting species, and several of these species appeared to respond favorably to treatments, in particular the MIS. The MIS are of particular interest for restoration references, since they are by definition found only on mobile dunes and are therefore most likely to decline across the reserve as a result of widespread dune fixation. Therefore, the positive albeit weak response of these species is encouraging. This research highlights the potential use of IS in assessing the impact of restoration projects, and emphasizes the importance of long-term studies in restoration projects.

Having longer time-series for pre-restored sites would be beneficial [138]. We strongly encourage LTER sites monitoring biodiversity trends to adopt a regularized temporal sampling regime that enables powerful and simpler statistical analyses [73,139]. That said, empirical ecological research projects in natural conditions are often unable to achieve perfect experimental design due to historical and external factors affecting accessibility. Balanced sampling across multiple taxa in a long-term study is invariably difficult for empirical data [140]. The approach to concurrently monitor treated, control and reference dunes allowed us to make effective interpretations, despite the unbalanced design. We hope the statistical approaches we adopted for unbalanced data will be useful for other researchers and practitioners.

Shrub removal experiments allowed us to investigate the structural component of dune states rather than the biological impact of the primary productivity. It is apparent that the shrubs themselves determine biodiversity trends rather than just providing biotic resources. Each taxon responded differently to removal, and each type of treatment had different effects. Restoration manipulations were only partially successful in Nizzanim LTER, but rodents seemed to respond most favorably. With respect to ecosystem functioning, "long-term" can mean anything from decades to centuries, so realistic expectations about restoration of sand dunes could take many years [141]. Over a longer timeframe, the response at the composition level may become more pronounced. Overall, it was apparent that in the absence of any disturbance, gamma diversity of coastal dunes communities would be lost as dune fixation increases.

It is likely that more 'natural' disturbances such as controlled grazing could be more effective and less damaging than mechanical ORV disturbances. Indicator species can represent assemblage level responses, although their response was often more pronounced than assemblages. Finally, it should be remembered that long-term assessment and monitoring of restored coastal dunes should include multiple taxa, and should examine changes in composition rather than abundance and richness alone.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3417/10/7/2310/s1, Table S1: Experimental design for the long-term ecological research project in Nizzanim, showing the sampling time (years) for each taxon, Table S2: Poisson (GLMM\_PQL) and Linear (LME) Mixed effect Models of community abundance (Figure 2) for each taxon, Table S3: Poisson (GLMM\_PQL) and Linear (LME) Mixed effect Models of species richness (Figure 3) for each taxon, Table S4: Pairwise Permanova of species composition among dune states. Each treatment type was tested separately, with some data overlap (e.g., mobile and Untreated Semi-Fixed dunes are used for both treated Semi-Fixed and Disturbed dune comparisons). Bold rows are of interest in terms of meeting the expectations for responses to treatment, Table S5: Poisson (GLMM\_PQL) and Linear (LME) Mixed effect Models of abundance for each Mobile-dune Indicator Species (MIS). Results were used for Figure 11a–d, Table S6: Poisson (GLMM\_PQL) and Linear (LME) Mixed effect Models of abundance for each Fixed dune Indicator Species Fixed-dune Indicator Species. Results were used for Figure 11e–h. Table S7: Permuted Empirical Cumulative Distribution Function (ecdf) models for Mobile dune Indicator Species (MIS) in each treatment type (Figure 12a–d), Table S8: Permuted Empirical Cumulative Distribution Function (ecdf) models for Fixed dune Indicator Species (FIS) in each treatment type (Figure 12e–h).

**Author Contributions:** Conceptualization, A.B., E.G. and P.B.K.; data curation, T.L.F.B., A.B., E.G. and P.B.K.; formal analysis, T.L.F.B.; funding acquisition, A.B., E.G. and P.B.K.; investigation, T.L.F.B., A.B., E.G. and P.B.K.; methodology, T.L.F.B., A.B., E.G. and P.B.K.; project administration, T.L.F.B. and P.B.K.; resources, A.B., E.G. and P.B.K.; supervision, A.B., E.G. and P.B.K.; visualization, A.B., E.G. and P.B.K.; writing—original draft, T.L.F.B.; writing—review and editing, A.B., E.G. and P.B.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by The Israel Nature & Parks Authority.

**Acknowledgments:** Immense that's are due to Yael Zilka, Arnon Tsairi, Boaz Shacham, Zehava Siegal, Ittai Renan, Adi Ramot and the numerous team heads, graduate and undergraduate students from Ben Gurion University for their assistance in collecting monitoring data across the years. Special thanks to Michael Dorman, who provided substantial time and advice on the initial analyses and use of R coding. We thank the entomologists at the Steinhardt Museum of Natural History for their help with identify of beetles, in particular Prof. Chicatunov & Laibale Friedman. Thanks to the Shikmim Field Study Center (Society for the Protection of Nature in Israel), for their hospitality during our fieldwork sessions over the years. We also acknowledge contribution from the

Ministry of Science and Technology (MOST) and the International Arid Land Consortium (IALC). Finally, we thank Yehoshua Shkedi, Yariv Malihi and Israel Nature & Parks Authority (INPA) rangers for continuous support and assistance.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Characteristics of Aeolian Dune, Wind Regime and Sand Transport in Hobq Desert, China**

**Hui Yang 1, Jiansheng Cao 1,\* and Xianglong Hou 1,2,\***


Received: 30 October 2019; Accepted: 12 December 2019; Published: 16 December 2019

**Abstract:** A systematic study of the wind regime characteristics in a region can not only accurately grasp the dynamic factors of the development of aeolian geomorphology, but also provide a scientific basis for the prevention and treatment of regional sand disasters. Taking the Hobq Desert as the study area, the basic characteristics of dune are analyzed by using remote sensing images. Based on the annual meteorological data of six meteorological stations from 2009 to 2018, the spatial and temporal distribution characteristics of wind speed were obtained. With the daily wind data of three stations from 2009 to 2018, we have figured out the wind regime and sand transport characteristics of the Hobq Desert. The results show that the sand dune height of the Hobq Desert ranges large, the highest height is 5010 m and the lowest is 10 m. It decreases gradually from the west to the east. The height of dune mainly distributed below 1500 m, followed by 1500–2000 m. Migratory sand dunes in Hobq Desert accounts for 51.8% and is mainly distributed in the west of the desert. The distribution area of fixation sand dunes in Hobq Desert is the least, accounting for 8.3%. The migratory dune pattern is trellis dune, semimigrated dune and semifixed dune patterns include honeycomb dune, parabolic duneand brush dune, and fixation dune pattern is grass dune. Annual wind speed was greatest in the southeast and decreased moving to the northwest. The dominant wind direction was W and SW from 2009 to 2018 in the Hobq Desert, the average wind speed of the prevailing winds mainly distributed at 4–8 m/s. The frequency of wind speed exceeding 10 m/s is very low, with a maximum value of 10% or below. There is a low energy wind environment in the Hobq Desert, with intermediate annual directional variability and obtuse or acute bimodal wind regime. The resultant drift direction (RDD)at Dongsheng station was relatively constant from 2009 to 2018, it was about 350◦. RDD differed significantly at Baotou and Linhestations were 181 ± 169◦ and 231 ± 121◦, respectively.The relationship between drift potential (DP) and the average and maximum wind speed was expressed as a power function. DP was strongly correlated with them. There is no significant correlated between the temporal changes in DPandprecipitation and temperature from 2009 to 2018 in the Hobq Desert.

**Keywords:** wind erosion; prevailing wind; sand transport; Hobq Desert; Yellow River

#### **1. Introduction**

Wind is one of the fundamental forces shaping geomorphology [1], especially in arid climate, it provides the main driving power to determine the desert surface morphology [2,3]. The distribution of aeolian sand, dune morphology, movement and even the development of grain sand are all related to the wind of different scales [4]. Wind is also the direct power condition causing aeolian sand disasters [5]. Wind energy is a key variable in indices designed to examine potential sand drift [6]. Wind regime (i.e., wind speed and direction) hastherefore been widely studied to support the classification of aeolian geomorphology [3]. Other regional factors, such as grain size and sand surface moisture [7], surface crusting [8], vegetation cover [9] and topography [10] also influence the aeolian sand transport. Scholars have previously made numerous efforts to study the estimation of sand transport through both field measurement and theoretical calculation with models [11–14]. Resulting from regional variability in these factors, the patterns, size and distribution of aeolian dunes usually exhibit distinct regional characteristics [15].

The Hobq Desert (Hobq is in Mongolian) is the seventh largest desert in China. At the beginning of 1959, the sand control team of the Chinese Academy of Sciences carried out a comprehensive investigation of the natural geographical features and natural resources of the Hobq Desert. Previous studies were mainly concentrated on the problem of land desertification in Hobq region [16], variation in grain size of dune [17], the environmental changes [18,19] (Li et al., 2017; Xu et al., 2018), and how the Hobq Desert affected the channel evolution of the Yellow River [20,21]. Ta et al. [22] suggested that channel aggradation of the Hobq Desert reach of the Yellow River was induced by perennial sediment input from the surrounding deserts. There has less sufficient detailed research on basic dune characteristics and temporal changes in wind regime and sand transport of the Hobq Desert.

The aims of this study were to analyze the dune height, dune type and spatial and temporal distribution of wind speed and perform an integrated analysis of recent changes from 2009 to 2018 in the wind regime of the Hobq Desert, and to discuss the influence of climate change on sand transport. The conclusions provide a scientific basis for the prevention and treatment of regional sand disasters and aeolian sand blown into the Yellow River, and the research application is the procedure for long-term monitoring of sand dunes.

#### **2. Study Area**

The Hobq Desert is located in the south bank of the Yellow River and the northern part of the Ordos Plateau in Inner Mongolia, with an area of 1.56 <sup>×</sup> 104 km2. The desert is narrow in the east and wide in the west, with a length of about 300km from the east to the west. Ten tributaries (cross-desert ephemeral tributaries) flow through the east and middle of the Hobq Desert to the Yellow River on its right bank (Figure 1). Set the most western tributary of the ten tributaries as the boundary, the eastern desert is with a width of 8–28km and that of the western desert is 40–90 km. Influenced by the east Asian summer monsoon, the average annual precipitation of the Hobq Desert decreased from 400 mm in the east to 150 mm in the west [23]. The rainfall appears mainly in the form of rainstorms and concentrates from July to September. Strong winds and sandstorms often happen in winter and spring with an annual average frequency of 34.5 and 19.2 days at Dongsheng meteorological station. At Baotou meteorological station, they are 46.8 and 21.6 days respectively [24]. The average day of sandstorms from 1981 to 2013 in Linhe weather station is 20 days [25].

**Figure 1.** Location of study area.

#### **3. Data and Methods**

#### *3.1. Data*

Dune height data are clipped from the map of desert distribution in 1:2,000,000 in China and dune type data are extracted from 1:100,000 scale desert distribution map set of China. These two data setsare both provided by Environmental and Ecological Science Data Center for West China, National Natural Science Foundation of China(http://westdc.westgis.ac.cn).

Remote sensing data used are the satellite images of the Google Earth. They are available and suitable for this study taken at some times in 2019.

Annual meteorological data (average wind speed, precipitation, temperature, relative humidity and annual maximum wind speed) of six meteorological stations around the Hobq Desert: Hangjinhouqi (HJHQ), Linhe (LH) and Baotou (BT) stations in the north and E'tuokeqi (ETKQ), Yijinhuoluo (YJHL) and Dongsheng (DS) stations in the south (Figure 1), and daily meteorological data of Linhe, Baotou and Dongsheng stations from 2009 to 2018 were downloaded from China's National Meteorological Information Center (http://data.cma.cn). Wind data were measured at 10 m above the ground, following World Meteorological Organization standards for anemometer heights.

#### *3.2. Sand Drift Potential*

Drift potential (DP) represents yearly total wind power and therefore is used for describing the potential maximum amount of sand transport by winds with a velocity above the threshold velocity [26]. It is calculated as follows:

$$DP = V^2(V - V\_1)t,\tag{1}$$

where *V* is the measured 10-min wind velocity, in m/s; *Vt* is the threshold wind velocity for winds to entrain sediment, in m/s, it was estimated to be 6 m/s in Hobq Desert, according to the previous study [24] and *t* is the proportion of time during which wind velocity is greater than *Vt*.

Modified by Fryberger and Dean [26], different DP values represent different wind energy environments. If DP is less than 200, the wind energy environment is low. If DP is between 200 and 400, the wind energy environment is intermediate. If DP is greater than 400, that is high.

Resultant drift potential (RDP) represents the net DP or the vectorial sum of the DP values in each compass direction. It is calculated as follows [6]:

$$RDP = \left(C^2 + D^2\right)^{0.5},\tag{2}$$

$$C = \sum \left( V \mathcal{U} \right) \sin(\theta),\tag{3}$$

$$D = \sum \left( V \mathcal{U} \right) \cos(\theta), \tag{4}$$

where *VU* represents the *DP* in each wind direction (in this paper, we grouped winds into 16 directions), in vector units, and θ is the midpoint of each wind orientation class measured clockwise from 0◦ (north) [3].

The resultant drift direction (RDD) represents the direction of net trend of sand drift. It is calculated as follows [6]:

$$RDD = \arctan(\text{C/D}).\tag{5}$$

Directional variability (RDP/DP) is the ratio of the resultant drift potential (RDP) to the drift potential (DP). RDP/DP values close to 1 indicate a narrowly unidirectional wind regime, with a single dominant drift direction, whereas values close to 0 indicate a multidirectional wind regime with multiple significant drift directions [3,6]. Specifically, if RDP/DP is less than 0.3, the directional category is complex or obtuse bimodal. If RDP/DP is between 0.3 and 0.8, that is obtuse bimodal or acute bimodal. If RDP/DP is greater than 0.8, that is wide or narrow unimodal [26].

#### *3.3. Trend and Breakpoint Detection*

The Mann–Kendall (M–K) trend test [27,28] is based on the correlation between the ranks of a time seriesand their time order. For a time series *X* = {*x*1, *x*2, ... , *x*n}, the test statistic is given by

$$S = \sum\_{i$$

where

$$a\_{ij} = \text{sign}(\mathbf{x}\_j - \mathbf{x}\_i) = \text{sign}(\mathbf{R}\_j - \mathbf{R}\_i) = \begin{cases} 1, \mathbf{x}\_i < \mathbf{x}\_j \\ 0, \mathbf{x}\_i = \mathbf{x}\_j \\ -1, \mathbf{x}\_i > \mathbf{x}\_j \end{cases} \tag{7}$$

and *Ri* and *Rj* are the ranks of observations *xi* and *xj* of thetime series, respectively. The mean and variance of the *S* statistic in Equation (6) above are given by [28]

$$E(S) = 0,\tag{8}$$

$$V0(S) = n(n-1)(2n+5)/18,\tag{9}$$

where *n* is the number of observations. If *n* is greater than 8, the *S* values will follow an approximately normal distribution [29]. The Mann–Kendall *Z* is given by [30]

$$Z = \begin{cases} (S-1) / \sqrt{V0(S)}, S > 0 \\ 0, S = 0 \\ (S+1) / \sqrt{V0(S)}, S < 0 \end{cases} \tag{10}$$

A positive value of *Z* indicates an increasing trend, and vice versa. The null hypothesis of no trend is rejected if |*Z*| > 1.96 at the 5% significance level [29].

After trend analysis, we use the Pettitt test [31], a non-parametric method that is widely applied to detect abrupt changes, to find out the breakpoint.

#### *3.4. Other Methods*

To calculate the spatial distribution of dune height inthe study area, we used ArcGIS10.5 to interpolate values between adjacent stations using the ordinary kriging function provided by the Spatial Analyst tool.

Spatial distribution of the average and maximum annual wind speed was calculated through Surfer 13.

We obtained mathematicsstatistical summary of the wind velocity parameters using OriginPro20, including standard deviation, minimum value and maximum value. We performed regressions for the relationships between *DP* and climatic parameters with SPSS 20.

#### **4. Results and Analysis**

#### *4.1. Basic Characteristics of Dune*

The sand dune height of the Hobq Desert ranged largely, the highest value was 5010 m and the lowest was 10 m. It decreased gradually from the west to the east (Figure 2A). The height of the dune mainly distributed below 1500 m, which accounted for 20.6%, followed by 1500–2000 m, which was 17.3%. The distribution area of other heights was basically the same, which was about 15%.

**Figure 2.** Spatial distribution of dune height (**A**) and dune type (**B**) in Hobq Desert.

According to different fractional vegetation cover (FVC), the Hobq Desert includes four dune types, migratory dune (FVC < 5%), semimigrated dune (FVC 5%–20%), semifixed dune (FVC 21%–50%) and fixation dune (FVC > 50%). Migratory, semimigrated, semifixed and fixation sand dunes in Hobq Desert accounted for 51.8%, 25%, 14.9% and 8.3% of the dunes, respectively. Migratory and semimigrated dunes mainly were distributed in the west and southwest of the desert (Figure 2B).

Twenty locations were selected, namely, sixmigratory dunes (M1–M6), five semimigrated dunes (SM1–SM5), seven semifixed dunes (SF1–SF7) and two fixation dunes (F1–F2). Location SF4–SF7, M5–M6 and SM5 lay in the eastern part of the desert, and the rest of the locationslay in the western part of the desert (Figure 2B). Downloading part of high-resolution images from Google Earth (Figure 4), following the Livingstone and Warren (1996) dune classification system, the migratory dune pattern wastrellis dune (M2, M4 and M6), semimigrated dune and semifixed dune patterns including the honeycomb dune (SM1 and SF1), parabolic dune (SM2 and SF3) and brush dune (SM3 and SF6), and the fixation dune pattern was grass dune (F1 and F2).

**Figure 3.** *Cont.*

**Figure 4.** *Cont.*

**Figure 4.** *Cont.*

**Figure 4.** *Cont.*

**Figure 4.** *Cont.*

**Figure 4.** Google Earth high-resolution imagery of different dune field patterns at 20 locations in the study area.

#### *4.2. Annual Wind Characteristics*

For the rest of the stations, the average and maximum annual wind speed fluctuated between about 10 m/s and 15 m/s from 2009 to 2018 (Figure 5).

**Figure 5.** Maximum (**a**,**b**) average annual wind speed at the six meteorological stations from 2009 to 2018 (stations are shown by abbreviated names: HJHQ—Hangjinhouqi station, BT—Baotou station, LH—Linhe station, ETKQ—E'tuokestation, DS—Dongshengstation, YJHL—Yijinhuoluostation).

With the M–K test, there was a significant downward trend of the maximum and average annual wind speed at Dongsheng station during 2009–2018 and had a break in 2013. At Hangjinhouqi and E'tuoke stations, there wasa significant upward trend of the average annual wind speed and had a break in 2017 and 2012, respectively. There was no significant trend at the rest of the stations (Table 1).


**Table 1.** Trend and breakpoint detection of maximum and average annual wind speed from 2009 to 2018 (\* referring to Figure 5 for the abbreviated names of stations).

The greatest maximum and average annual wind speed were both at Yijinhuoluo (16.7 <sup>±</sup> 3.1 m·s−<sup>1</sup> and 3.0 <sup>±</sup> 0.5 m·s<sup>−</sup>1). The smallest of those two parameters were at Dongsheng (10.0 <sup>±</sup> 1.4 m·s<sup>−</sup>1) and Hangjinhouqi (2.2 <sup>±</sup> 0.2 m·s<sup>−</sup>1), respectively (Table 2).


**Table 2.** Statistical summary of wind speed (\* referring to Figure 5 for the abbreviated names of stations).

**\*** SD, standard deviation. Min., the minimum value. Max., the maximum value.

Figure 6 shows the distribution of the average and maximum annual wind speed in the desert and adjacent areas from 2009 to 2018. Annual wind speed was greatest in the southeast and decreased moving to the northwest.

#### *4.3. Dominant Wind Characteristics*

In order to find out the dominant wind direction, we obtained the rose diagrams of the wind direction at Dongsheng, Baotou and Linhe stations during the years from 2009 to 2018 with daily wind data (Figure 7). The dominant wind direction of Dongsheng station was W from 2009 to 2016, and WNW from 2017 to 2018, but the frequency of W wind was less than 3% lower than that of WNW wind.Apart from the ENE wind in 2009, all the other years of Linhe station were SW wind, but after 2010 and 2014, it was mainly WSW. Different from Dongsheng and Linhe stations, the dominant wind direction of Baotou station varied greatly in the decade. To be specific, the prevailing wind was NW from 2009 to 2012, and ESE from 2013 to 2018.

**Figure 6.** Distribution of the average (**a**) and maximum (**b**) annual wind speed (m·s<sup>−</sup>1) in the Hobq Desert and adjacent areas.

According to the prevailing winds of Dongsheng, Baotou and Linhe stations in each year, the frequency of occurrence of different speed intervals was plotted (Figure 8). From the overall wind speed distribution, the average wind speed of the prevailing winds of each station from 2009 to 2018 was mainly distributed at 4–8 m/s. In contrast, the frequency percentage of Baotou station with wind speed greater than 6m/s was the lowest, the mean value was less than 30% and the maximum value wasonly 54%. As for the other two stations, they had little difference, the frequency percentage of the maximum wind speed greater than 6 m/s was about 85% and that of the mean value was about 45%. However, the frequency percentage of wind speed exceeding 10 m/s at each station was very low, with the maximum value of 10% or below.

**Figure 7.** Variations of prevailing winds at Dongsheng, Baotou and Linhe stations from 2009 to 2018 in the Hobq Desert (referring to Figure 5 for the abbreviated names of stations).

**Figure 8.** Distribution of the wind speed of dominant wind direction at Dongsheng (**a**), Baotou (**b**) and Linhe (**c**) stations from 2009 to 2018 in the Hobq Desert (referring to Figure 5 for the abbreviated names of stations).

#### *4.4. Sand Drift Potential*

Sand drift potential (DP) was greatest at Baotou station, where the annual DP was 76.79 ± 68.57 (mean ± SD), followed by the Linhe station, with an annual DP of 45.3 ± 28.55. That was 17.43 ± 11.41 at Dongsheng station (Figure 9a). All those DP values indicated a low energy wind environment in the Hobq Desert.

**Figure 9.** Calculated results of drift potential (DP; (**a**)), resultant DP (RDP)/DP (**b**) and resultant drift direction (RDD; (**c**)) at selected stations Dongsheng, Baotou and Linhe stations from 2009 to 2018 in the Hobq Desert (referring to Figure 5 for the abbreviated names of stations).

The annual directional variability was mostly intermediate, with values between 0.3 and 0.8 (Figure 9b), with one exception: Baotou had low directional variability in 2018 and Linhe had low directional variability in 2013 and 2015. The wind regime in the Hobq Desert was main obtuse or acute bimodal.

RDD at Dongsheng station was relatively constant from 2009 to 2018, it was about 350◦, and the standard deviation was 12◦. RDD differed greatly during 2009 to 2018 at Baotou and Linhe, were 181◦ ± 169◦ and 231◦ ± 121◦, respectively (Figure 9c).

#### *4.5. Sand Transport and Wind Speed*

Using the annual wind speed data from threeweather stations in and around the Hobq Desert, we obtained 10 years of data that provided a better idea of the annual maximum and average wind speed and established a relationship between DP and the average (Figure 10a) and maximum wind speed (Figure 10b). It could be concluded that this relationship was expressed as a power function through comparing the values of *R*<sup>2</sup> and *p* among 11 kinds of math models. DP was strongly correlated with the annual average wind speed (*R*<sup>2</sup> = 0.81, *p* < 0.01) and the annual maximum wind speed (*R*<sup>2</sup> = 0.90, *p* < 0.01). With the increase of wind speed, DP increased, especially the average wind speed was greater than 2.5 m/s and the maximum wind speed was greater than 12 m/s, DP increased rapidly.

**Figure 10.** The relationship between the calculated drift potential (DP) and (**a**) theannual average and (**b**) the maximum annual wind speed from 2009 to 2018 (*vave* and *vmax*, respectively).

#### *4.6. Climate Change Characteristics*

We used precipitation and temperature to present the climate change characteristics. Based on data from three weather stations (Dongsheng, Baotou and Linhe) from 2009 to 2018, we analyzed annual precipitation and annual average temperature characteristics of those three stations. The annual precipitation at Linhe station was 200 mm and below. While, the annual precipitation at Dongsheng and Baotou stations was between 200 and 500 mm, for some year, it was as high as about 700 mm at Dongsheng station (Figure 11a).

**Figure 11.** Characteristics of annual precipitation (**a**) and annual average temperature (**b**) from 2009 to 2018 at three representative weather stations in the Hobq Desert (referring to Figure 5 for the abbreviated names of stations).

For the annual average temperature, it concentrated from 8 to 8.5 ◦C at Baotou station (Figure 11b). The annual average temperature ranged greater at Dongsheng and Linhe stations, were 6.5–8.5◦C and 7.5–9.5◦C, respectively.

#### **5. Discussion**

Wind is the power responsible for movement of aeolian sediments, thus transport should increase with increasing wind speed. However, other factors can disrupt this simple relationship. For example, precipitation can increase vegetation cover, sand moisture content and the air's relative humidity, and all of these factors can decrease sand transport [32,33]. Increased temperatures have the opposite effect: they increase evaporation, and thereby decrease soil cohesion. Thus, climate change will affect the formation and subsequent evolution of sand seas and sand dunes [3]. To figure out the influence of long-term climate change on sand transport, we analyzed the relationship between DP and annual precipitation and annual average temperature. The annual precipitation at Linhe station was positive correlated with DP. While, it was negative correlated with DP at the eastern part of the Hobq Desert (Figure 12a).

**Figure 12.** The relationships between the calculated DP and (**a**) annualprecipitation (Pre) and (**b**) annual average temperature (Tem) from 2009 to 2018 at three representative weather stations in the Hobq Desert.

For the annual average temperature, it had a positive correlation with DP at Baotou station (Figure 12b) and the annual average temperature at Dongsheng and Linhe stations had a negative correlation with DP. However, the relationships between DP and the annual precipitation and temperature at all the stations from 2009 to 2018 were not significant (*p* > 0.05). Although precipitation and temperature directly restricted the species diversity, growth and development of plants, usually, the vegetation developed in dry years is difficult to produce coverage protection effect on the surface [34]. For the Hobq Desert with a large area of migratory dune, there wasno significant correlationbetween the temporal changes in DPand climate change. The modern environmental monitoring of the Maowusu Sandy Land, which is located in the south of the Hobq Desert shows that the vegetation restoration capacity of the east and west of the Maowusu Sandy Land was out of sync with the decreasing wind, decreasing anthropology activities and slightly increasing precipitation. The eastern region gives a priority response to the climate changes, while the western region depends on the continuous improvement of climatic environment [35]. This suggests that the response to climate change in the eastern and western regions of the Hobq Desert may also be out of sync. Under the condition of a better climate, the eastern region may give a priority response to the changes, solidifying the dune and gradually developing the paleosols, while the sandstorm activity in the western region may still continue, at least it did not form the environmental conditions for soil development. Meanwhile, through study on the sedimentary sequence of the strata in the Hobq Desert, it is found that the formation of aeolian sand in some sections is mainly related to the increase of available loose sand sources rather than the arid climate conditions [23].

#### **6. Conclusions**

Wind regimes play an important role in the formation and development of aeolian dune. The calculated DP was strongly related to the annual average and maximum wind speed, and the relationships were described with a power function. With the increase of wind speed, DP increased, especially the average wind speed was greater than 2.5 m/s and the maximum wind speed was greater than 12 m/s, DP increased rapidly. While, considering the other two key climate factors, precipitation and temperature, there was no significant correlation between them and calculated DP from 2009 to 2018.

In the Hobq Desert, there was a low energy wind environment, with intermediate annual directional variability and obtuse or acute bimodal wind regime. Annual wind speed was greatest in the southeast and decreased moving to the northwest. Aeolian dune height mainly distributed below 1500 m, it decreased gradually from the west to the east, and dune patterns includedtrellis dune, honeycomb dune, parabolic duneand brush dune and grass dune.The dominant wind direction was main W and SW from 2009 to 2018 in the Hobq Desert, the average wind speed of the prevailing winds mainly distributed at 4–8 m/s. The frequency of wind speed exceeding 10 m/s is very low, with the maximum value of 10% or below.The complexity of the local and region wind regime controlled the variety of the sand dune patterns.

**Author Contributions:** H.Y. performed the data processing and wrote the manuscript. X.H. provided the ideas and helped process the data. J.C. revised the language and provided financial support. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Plan Project of China grant numbers [2018YFC0406501–02] and National Natural Science Foundation of China grant numbers [41877170]. And the APC was founded by [2018YFC0406501–02].

**Acknowledgments:** This work is supported by CFERN and Beijing Techno Solutions Award Funds on excellent academic achievements. We are grateful for the inputs of the laboratory staff and anonymous reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Wind Tunnel Measurements of Surface Shear Stress on an Isolated Dune Downwind a Bridge**

**Wenbo Wang 1,2, Hongchao Dun 1,\*, Wei He <sup>1</sup> and Ning Huang 1,\***


**\*** Correspondence: dunhc@lzu.edu.cn (H.D.); huangn@lzu.edu.cn (N.H.)

Received: 2 May 2020; Accepted: 8 June 2020; Published: 10 June 2020

**Abstract:** As part of a comprehensive environmental assessment of the Dun-Gel railway project located in Dunhuang city, Gansu Province, China, a wind tunnel experiment was proposed to predict surface shear stress changes on a sand dune when a bridge was built upstream it. The results show that the length of the wall shear stress shelter region of a bridge is about 10 times of the bridge height (H). In the cases that the interval of the bridge and sand dune (S) is less than 5 H, normalized wall shear stress on the windward crest is decreased from 1.75 (isolated dune) to 1.0 (S = 5.0 H, measured downwind bridge pier) and 1.5 (S = 5.0 H, measured in the middle line of two adjacent bridge piers). In addition, the mean surface shear stress in the downstream zone of the sand dune model is reduced by the bridge pier and is increased by the bridge desk. As for the fluctuation of surface shear stress (ζ) on the windward crest, ζ decreases from 1.3 (in the isolated dune case) to 1.2 (in the case S = 5.0 H, measured just downwind the pier) and increases from 1.3 (in the isolated dune case) to 1.6 (in the cases S = 5.0 H, in the middle of two adjacent piers). Taking the mean and fluctuation of surface shear stress into consideration together, we introduce a parameter ψ ranging from 0 to 1. A low value indicates deposition and a high value indicates erosion. On the windward slope, the value of ψ increases with height (from 0 at toe to 0.98 at crest). However, in the cases of S = 1.5 H, ψ is decreased by the bridge in the lower part of the sand dune at y = 0 and is increased at y = L/2 compared with the isolated dune case. In other cases, the change of ψ on the windward slope is not as prominent as in the case of S = 1.5 H. Downstream the sand dune, erosion starts in a point that exists between x = 10 H and 15 H in all cases.

**Keywords:** wind tunnel; surface shear stress; isolated dune; bridge

#### **1. Introduction**

Railways in arid and semiarid regions in Western China suffer from sand hazards [1–3]. On one hand, wind-blown sand movement could bury rail tracks, attack railway power systems and roll over the trains. On the other hand, railway structures such as high subgrades, bridges and windbreak walls are huge enough to perturb airflow over the nearby sand dunes, which are stable and moving regularly under initial airflow conditions [4–6]. Therefore, it is important to study how the railway structures change the shape and moving patterns of sand dunes in the neighborhood. In previous studies, particular attention has been paid to the effects of wind-blown sands on the railways and subgrades [5–8], however, the effects of railways on the wind-blown sand, especially on the sand dunes movement, are not well explored. Based on the bridge engineering of the Dun-Gel railway in a sand valley, we observed and characterized the surface shear stress on the sand dunes downwind of a bridge in the wind tunnel.

Sand particle deposition and erosion processes are closely related to the threshold friction velocity of sand entrainment [9–14]. Many studies have shown the relationships between surface shear stress acting on the ground and sand entrainment [15–17]. To understand the formation and migration of dunes, one first needs to know the stationary wind stress exerted on a given sand topography [18]. Using the visualization techniques to investigate the flow structures around the roughness elements, researchers have shown that vortex strongly governs the pattern and magnitude of surface shear stress at the presence of roughness elements [13]. Hence, it is important to study the fluctuation of surface shear stress.

This study is a part of a comprehensive environmental assessment of the Dun-Gel railway project located in Dunhuang city, Gansu Province, China (Figure 1). Located in the northwest of China, Gansu lies between latitude 32◦11 –42◦57 north and longitude 92◦13 –108◦46 east, with a total area of 425,800 km2. The terrain of Gansu province is long and narrow, and the landform is complex and diverse, including mountains, plateaus, flat rivers, river valleys, deserts and gobi. The climate type of Gansu is mainly temperate continental arid climate. Our purpose is to study how the railway structures affect the shape and moving patterns of sand dunes in the neighborhood by investigating the shear stress on the sand dune surface, which determines the deposition or erosion of the sand dune surface. In this study, spatial variation in the surface shear stress was measured directly by mounting Irwin sensors at numerous discrete points on the sand dune model [19]. The relationship between the surface shear stress on the sand dune and the position of the bridge was analyzed. By analyzing spatial variation changes of surface shear stress on and downstream of the sand dune model, the erosion and deposition zones were identified. Moreover, future changes in the shape and moving patterns of the sand dunes can be predicted.

**Figure 1.** Dun-Gel railway project in Dunhuang city, Gansu Province, China.

#### **2. Methods**

#### *2.1. Experimental Setup*

The experiment was carried out in a multi-functional environment wind tunnel of Lanzhou University. This open-return blow-down low-speed wind tunnel was 22 m long (only for work section) with a cross-section of 1.45 m × 1.3 m [20]. The roughness element in front of the wind tunnel was used to accelerate the development of the boundary layer, which was about 0.2 m thick in the measurement section. The topographic map was measured with an unmanned aerial vehicle, and the dispersion between the crest and the windward slope foot was 10 m, the height of the bridge in front of the sand dune was 12 m. The crests of sand dunes in Figure 1 were parallel to the bridge, suggesting that the prevailing wind direction and dune movement direction were perpendicular to the bridge. In these

types of studies, the Reynolds similarity was always unsatisfied because the roughness height in the wind tunnel was 1–2 orders smaller than that in the atmosphere boundary layer. As a result, researches have to neglect the Reynolds similarity [21–23].

Figure 2 shows the schematic of the models and the coordinate system used in this study. Roughness elements were placed 6 m upstream of the working section to generate a turbulent boundary layer. The plywood model of the sand dune was 10 cm in height (H) and 100 cm in width. The steel model of the bridge was 12 cm in height (H') and 36 cm between two adjacent piers (L). The scale between models and real dimensions of dune and bridge is 1:100. To observe the changes of surface stress on the sand dune when the bridge was installed and moved upstream, the distance between the bridge and the sand dune (S) was set as 1.5 H, 3.0 H, 5.0 H, 10.0 H and 20.0 H, respectively. A pitot tube was used to measure the profile of flow speed. Irwin sensors were mounted on the ground and sand dune models, which were to measure surface shear stress.

**Figure 2.** Schematic diagram of wind tunnel experiment.

#### *2.2. Surface Shear Stress Sensor*

Irwin sensor is a simple omnidirectional pressure meter that determines the friction velocity without the requirement of alignment by measuring the near-surface vertical pressure differential, which can be used to estimate the surface shear stress in the complex flow [24,25]. It has been used successfully in aeolian sand researches in both lab [13,21,22,26] and field [27,28] to explore the relationship between shear stress and sand transport potential in complex, nonuniform airflow. Many studies have confirmed that the sensor can be simply calibrated [24,29] and is applicable to measure the surface shear stress on a sand dune downwind a bridge in the wind tunnel experiment.

Based on the Irwin empirical function, the friction velocity prepared to calibrate the 34 Irwin sensors was determined using the velocity profile technique. As shown in Equation (1), the calibration of one sensor is [29]:

$$\frac{\mu\_\* h}{v} = |a + b\left(\frac{\Delta p h^2}{\rho v^2}\right)^n|\tag{1}$$

where *h* is the height of the obstacle, *v* stands for kinematic viscosity. The calibration coefficients we used is an averaged result of 34 Irwin sensors. The parameter "*a*" equals 3.86, "*b*" equals 0.11 and "*n*" equals 0.52. The calibration result is accurate enough for the wind velocity within 2 m s−<sup>1</sup> and 18 m <sup>s</sup>−<sup>1</sup> [29]. The friction velocity *<sup>u</sup>*<sup>∗</sup> can be converted to shear stress using <sup>τ</sup> = <sup>ρ</sup>*u*<sup>2</sup> ∗.

Figure 3 is the schematic of installation of sensors over the sand dune. In brief, 34 Irwin sensors were used stream-wise. Among them, 8 were set on the windward slope, 5 were set on the leeward slope, 3 were mounted upwind the dune, and 18 were mounted downwind the dune. The measuring position changed three times along the *y*-axis at y = 0 cm, 6 cm (L/6) and 18 cm (L/2), respectively. One sensor was set at 17.5 H upwind the model to detect the shear stress on a bare surface (τ0) for a reference value to normalize surface shear stress (τ = (τ*<sup>x</sup>* − τ0)/τ0). The bridge model, sand dune model and the stream-wise wind velocity were symmetric to the centerline of the wind tunnel. Therefore, the distribution of shear stress on the sand dune was also symmetrically distributed with the centerline.

**Figure 3.** (**a**) Irwin sensor and the installation sites on dune model during the wind tunnel test, (**b**) Installation dimension of Irwin sensors.

#### *2.3. Procedure of Wind Tunnel Experiment*

The wind tunnel experiment was carried out according to the following steps:


**Figure 4.** The data of wind field in wind tunnel without models: (**a**) vertical distribution of time-averaged wind velocity measured by Pitot-static tubes with different *u*<sup>0</sup> and (**b**) fluctuations of friction velocity measured by Irwin sensor with different *u*0.


**Table 1.** The fitting results of wind tunnel experiments.

*<sup>u</sup>*<sup>0</sup> (m s−1) is free-stream wind velocity, *<sup>u</sup>*∗*<sup>a</sup>* (m s−1) is the friction wind velocity computed in profile technique, <sup>τ</sup>*<sup>a</sup>* (N m2) is the shear stress computed in profile technique, FF. *<sup>u</sup>*∗*<sup>b</sup>* (m s−1) is the friction velocity measured by an Irwin sensor, τ*<sup>b</sup>* (N m2) is the average shear stress measured by an Irwin sensor and στ*<sup>b</sup>* (N m2) is the standard deviation of the shear stress measured by an Irwin sensor.

#### **3. Results**

#### *3.1. Changes in Averaged Surface Shear Stress*

Figure 5 shows the spatial variability of the mean value of normalized surface shear stress (τ ) on the sand dune and downstream it. The windward slope, leeward slope and the crest of the sand dune can be easily recognized in the graphs. It should be noted that for saving the space, the distance between the bridge and the sand dune, S, in the graphs was not in actual scale. Figure 5a shows that in the isolated dune case, wall shear stress increased on the windward slope compared with that on the flat surface. Near the crest of the windward slope, τ was the highest and equals 1.5. In addition, a zone with negative τ was formed between x = 2.5 H and 10 H, corresponding to the backward velocity of the airflow (reattachment region was marked in Figure 5a with green lines). When x reaches 20 H, the magnitude of τ approaches 0, indicating that wall shear stress was no more affected by the sand dune. In Figure 5b–f, although bridge reduces wall shear stress on the windward slope, the magnitude of τ on the windward crest was still high (ranges from 1.25 to 1.5). Moreover, the bridge did not affect the shear stress on the leeward slope dramatically, that is, the change of shear stress on the leeward slope was no more than 5% compared with the isolated case at the same positions. In the cases of S = 1.5 H, 3.0 H and 5.0 H, wall shear stress between x = 5 H and 10 H is affected by the bridge pier and desk. That is, the wall shear stress measured at 5 H < x < 10 H, y = L/2 is higher than that measured at 5 H < x < 10 H, y = 0. Moreover, wall shear stress downwind the sand dune is not affected by the bridge in the case of S = 10.0 H and 20.0 H.

To analyze the change of normalized surface shear stress quantitatively, Figure 6 shows normalized surface shear stress (τ ) as a function of the interval distance (*S*) and stream-wise position. Figure 6a shows the variation of the normalized surface shear stress (τ ) along the centerline (*y* = 0). Compared with the isolated dune case, the surface shear stress on the windward slope was reduced due to shelter of the pier (except at the two points at the beginning of the windward slope). When the free-stream wind velocity *u*<sup>0</sup> was 15 m s−1, the τ on the windward crest decreased by 20.4%, 28.2%, 31.7%, 6.4% and 0.5%, corresponding to the cases of S = 1.5 H, 3.0 H, 5.0 H, 10.0 H and 20.0 H respectively. In the case that S equals 20.0 H, the shelter effect of the wall shear stress almost disappears. The shelter region length of the wall shear stress is between 10 and 20 times the bridge height, which is longer than that of a cylinder. We attribute this to the coeffects from the bridge pier and desk.

Figure 6b shows the variation of the normalized surface shear stress (τ ) along y = L/2 with the free-stream wind velocity *u*<sup>0</sup> = 15 m/s. Compared with the isolated dune case, shear stress on the crest decreases by 11.4%, 9.1%, 12.0%, 0.5% and −2.3%, corresponding to the cases of S = 1.5 H, 3.0 H, 5.0 H, 10.0 H and 20.0 H. In the downwind area of the sand dune, the shear stress shows some increase from x = 2.5 H to 15 H compared with the isolated dune case. At the downstream position of x = 20 H, the shear stress became stable and almost equaled the free-steam value.

**Figure 5.** Normalized value of shear stress τ with normalized distance x/H. H is the height of sand dune, L is the distance of piers of bridge, S is the distance between bridge and sand dune and y is the distance between measuring position and the y-axis.

**Figure 6.** Normalized shear stress τ with normalized distance x/H. H is the height of sand dune, L is the distance of piers of bridge, S is the distance between bridge and sand dune and y is the distance between measuring position and the y-axis.

#### *3.2. Changes in the Fluctuation of Surface Shear Stress*

To study the fluctuation of the shear stress measured on the sand dune, a normalized parameter ζ for the standard deviation of shear stress was defined as follows [13]:

$$\zeta = \frac{\sqrt{\overline{\tau\_z'^2}}}{\sqrt{\overline{\tau\_0'^2}}} \tag{2}$$

where τ *z* <sup>2</sup> <sup>=</sup> στ*<sup>z</sup>* is the standard deviation of the surface shear stress when the model was installed, and τ 0 <sup>2</sup> <sup>=</sup> στ<sup>0</sup> is the standard deviation of the surface shear stress measured on the flat floor.

Figure 7 shows the variation in the normalized standard deviation of the surface shear stress ζ with free-stream wind velocity equaled 15 m s<sup>−</sup>1. On the windward crest, magnitude of ζ equals 1.35 in the isolated dune case. However, in the cases that S = 1.5 H, 3.0 H and 5.0 H, parameter ζ measured at the position of y = 0 and L/2 shows some different trends. That is, on the windward crest, wall shear fluctuation at y = 0 is decreased compared with the isolated dune case. Quantitatively, equals 1.2 in the case of S = 1.5 H while ζ equals 1.35 in the isolated dune case. Downstream the sand dune, wall shear fluctuation measured both at y = 0 and L/2 is lower than the isolated case between x = 5 H to 10 H. However, wall shear fluctuation at y = L/2 is higher than at y = 0. That is, the bridge shows some restraint on the wall shear fluctuation downwind the sand dune and the restraint effect of the bridge pier is greater than that of the bridge desk. In addition, the influence from bridge and sand dune almost disappears at x = 15 H. In the cases that S equals 10.0 H and 20.0 H, ζ calculated on and downwind the sand dune shows little divergence at y = 0 and L/2 and the trend is very closed to the isolated case and this feature is very prominent in the case that S = 20.0 H.

**Figure 7.** Normalized fluctuation of shear stress ζ with normalized distance x/H. H is the height of sand dune, L is the distance of piers of bridge, S is the distance between bridge and sand dune and y is the distance between measuring position and the y-axis.

#### **4. Discussion**

The threshold friction velocity is a key parameter to estimate the deposition and erosion potential of the ground. Some experiments show a higher threshold friction velocity than the calculation results using the Bagnold's empirical equation [30]. We calculated the fluid threshold skin friction velocity as:

$$
u\_{\tau\tau} = \left. u\_{\ast t} + 2\sigma(u\_{\tau}) \right. \tag{3}$$

where *<sup>u</sup>*∗*<sup>t</sup>* was 0.29 m s−<sup>1</sup> in our experiment [31]. The standard deviation <sup>σ</sup>(*u*τ) <sup>=</sup> 0.048 m s−<sup>1</sup> at 15 m/<sup>s</sup> was determined from skin friction velocity variations measured with Irwin sensors on the smooth wooden floor. As a result, *u*ττ in our study was 0.387 m s−<sup>1</sup> and the threshold shear stress τ*<sup>t</sup>* equaled 0.179 Pa.

To assess the local dominance of erosion and deposition mechanisms, a threshold parameter of fraction time ψ was proposed [14], which represents the fraction of time when friction velocity exceeds the threshold parameter.

$$\psi(\mathbf{x}, y) = \frac{\Delta t[\pi(t, \mathbf{x}, y) > \tau\_t]}{T} \tag{4}$$

where Δ*t* is the time period during which the shear stress τ is larger than the threshold shear stress τ*t*, and *T* represents the total time. As defined ψ = 1 indicates an erosion zone, ψ = 0 indicates a deposition zone, and 0 <ψ< 1 indicates both erosion and deposition are both possible. Specifically, for ψ < 0.5, deposition dominates the local net and for ψ > 0.5, erosion dominates.

Figure 8 shows the spatial distribution of ψ at the free-stream wind velocity *u*<sup>0</sup> = 15 m s−<sup>1</sup> in our experiment. In the isolated dune case, the value of ψ increased from 0 at the windward toe to 0.98 at the windward crest, indicating that the erosion rate increases with height on the windward slope. On the leeward and downwind the sand dune (0 < x < 10 H), the values of ψ are less than 0.5, indicating that deposition happened in this region. In the case of S = 1.5 H, the value of ψ is decreased by the bridge pier and increased by the bridge desk in the lower part of the windward slope. That is, the value of ψ between x = -5 and -2.5 decreases at y = 0 and increases at y = L/2 compared with the isolated dune case. For instance, the value of ψ is 0.42 at x = −2.5, y = 0 and the value of ψ is 0.65 at x = −2.5, y = L/2, decreased and increased compared with the value of 0.5 in the isolated dune case. In the cases of S = 3.0 H, 5.0 H, 10.0 H and 20.0 H, the change of ψ on the windward slope is not as prominent as in the case of S = 1.5 H. On the leeward slope, ψ in this region always kept low magnitudes that are lower than 0.3, indicating that the incoming sand particles deposit no matter where the bridge was set. Downstream the sand dune, in the cases of S = 1.5 H, 3.0 H and 5.0 H, ψ is decreased at y = 0 and increased at y = L/2 compared with the isolated dune case. However, the value of ψ exceed 0.5 at x = 15 H almost in all cases. That is, between x = 0 and x = 15 H, deposition domains the sand movement and downwind x = 15 H, erosion domains the sand movement. In the cases of S = 10.0 H and 20.0 H, the trend of ψ is similar to the isolated dune case, especially when S = 20.0 H.

Many wind tunnel experiments for shear stress measurement on the sand dune and around buildings have been done for engineering applications [32]. However, results on the influence of a building on sand terrain are still lacking. To the authors' knowledge, the study of wind erosion patterns downstream a bridge has not been proposed yet. The results in the present study not only provides predictions on sand dune moving after a bridge which is built upstream, but also confirms the viewpoints that the size of the sand drifts is very sensitive to the frontage of the upwind collecting area, which is accordant with the comment proposed [9]. Our research on surface shear stress can play an important role in railway construction items, and can give some guidance on future research about the moving pattern of sand dunes affected by building upstream.

**Figure 8.** Fraction time threshold time parameter ψ with normalized distance x/H. H is the height of sand dune, L is the distance of piers of bridge, S is the distance between bridge and sand dune and y is the distance between measuring position and the y-axis.

#### **5. Conclusions**

We used 34 Irwin sensors to measure the surface shear stress on a sand dune model that is downstream a bridge model immersed in a fully developed turbulent boundary layer in a wind tunnel. The results showed that the averaged wall shear stress on the windward crest of the sand dune is decreased by the bridge in the cases of S = 1.5 H, 3.0 H and 5.0 H. The decrease ratio obtained at y = 0 is greater than that obtained at y = L/2.

The bridge also affected the fluctuation of surface shear stress on the sand dune. In the case of S = 1.5 H, 3.0 H and 5.0 H, the fluctuation of surface shear stress at the windward crests decreased at y = 0 and is increased at y = L/2 compared with the isolated dune case. Downstream the sand dune, the restraint effect of the bridge pier on the wall shear fluctuation is greater than that of the bridge desk.

The fraction time parameter ψ used in the study to assess the local dominance of erosion and deposition mechanisms indicated strong erosion on the windward slope. The value of ψ increases with height on the windward slope of an isolated sand dune. However, in the cases of S = 1.5 H, ψ is decreased by the bridge in the lower part of the sand dune at y = 0 and is increased at y = L/2 compared with the isolated dune case. In other cases, the change of ψ on the windward slope is not as prominent as in the case of S = 1.5 H. Downstream the sand dune, erosion starts in a point between x = 10 H and 15 H in all cases. However, differences in the value of ψ downwind the sand dune indicates the erosion rate can be different.

Last to be acknowledged, wind tunnel experiments of geometrically similar models cannot fulfill the Reynolds similarity. To get more precise results, the next study we concentrate on the field observation of shear stress and sand flux so that erosion patterns on the sand dune downstream a bridge. The simple assumption in this paper will be improved.

**Author Contributions:** Conceptualization, H.D. and N.H.; methodology, W.H.; investigation, W.W. and H.D.; resources, N.H.; data curation, W.W. and W.H.; writing—original draft preparation, W.W. and W.H.; writing—review and editing, H.D.; funding acquisition, N.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National key research and development projects (2016YFC0500901), National natural science foundation of China (11772143), Scientific and Technological Services Network Planning Project of Cold and Arid Regions Environmental and Engineering Research Institute, CAS (HHS-TSS-STS-1504), Central university fund (lzujbky-2020-cd06).

**Acknowledgments:** The authors would like to thank Jie Zhang, Lanzhou University, for the constructive suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Using Post-Harvest Waste to Improve Shearing Behaviour of Loess and Its Validation by Multiscale Direct Shear Tests**

#### **Wen-Chieh Cheng 1,2,\*, Zhong-Fei Xue 1,2, Lin Wang 1,2 and Jian Xu 1,2**


Received: 19 October 2019; Accepted: 27 November 2019; Published: 29 November 2019 -

#### **Featured Application: The research work highlights the great potential of application of loess-PHW mixtures to loess erosion control and slope stabilization.**

**Abstract:** Loess and PHW (post-harvest waste) are easily accessible in the Chinese Loess Plateau and have been widely applied to construction of residential houses that have been inhabited for decades under the effect of freeze-thaw cycles. Although many researchers have recognised that the addition of fibers to loess soil is effective in preventing soil erosion and stabilising slopes, a consensus on this claim has not been reached yet. This study investigates the shearing behaviour of the loess-PHW mixture using small-scale and large-scale direct shear (SSDS and LSDS) tests. Four typical shear stress versus horizontal displacement curves from the multiscale direct shear tests are recognised where one is featured with strain-softening shape and the other three with a strain-hardening shape. Two out of the three curves with strain-hardening shape show a gradual increase in the shear stress at additional and larger displacements, respectively, in which some factor starts to have an influence on the shearing behaviour. Comparisons of the shear strength measured in SSDS and LSDS are made, indicating that there are differences between SSDS and LSDS. The effect of PHW addition on shear strength is assessed in order to determine the optimal dosage. The improvement of shear strength is attributed to the effect of particle inter-locking, resulting from the addition of PHW to loess specimens, and takes effect as the water content surpassed a threshold, i.e., >14%, that facilitates particle rearrangement. Particle-box interaction behaviour is assessed at the same time, and the findings satisfactorily address the main cause of the gradual increase in shear stress following the curve inflection point. The improved shearing behaviour proves the ability of the loess-PHW mixture to resist the seepage force and consequently stratum erosion.

**Keywords:** direct shear test; post-harvest waste; Chinese Loess Plateau; optimal dosage

#### **1. Introduction**

Since rapid urbanisation has significantly impacted surrounding environments [1–7], environment-friendly construction material may ease the impacts despite many advanced technologies available for environmental protection [8–12]. The strength and microstructure properties of spent coffee grounds (CG) stabilised with rice husk ash (RHA) were investigated, which are organic wastes derived from agricultural products [13–15]. The results have shown that elevated temperature curing of up to 90 days was deemed the key to secure the strength development of CG-RHA geopolymers. Belhadj et al. [16] performed an experimental work to assess the influence of the addition of barley straw

on the physico-mechanical properties and the microstructure of a concrete consisting only of sand as the main aggregate. The addition of barley straw to the sand concrete greatly improved its thermo-physical properties. On the other hand, significant improvements have also been measured in other properties, such as flexural strength, lightness, deformability, ductility, and toughness, notwithstanding that a decrease in the mechanical strength and an increase in the dimensional variations have been recorded. To solve the problems raised, several treatments were tested for improving the properties of the optimal composition of the studied lightweight concrete, and the barley straws treated with hot water showed good improvements in the flexural strength of the composite. Additionally, the hot water treatment led to acceptable results in the thermal characteristics although the density of the composite did not increase much compared to the concrete containing the untreated barley straws [17].

Loess soils are aeolian deposits containing primarily silt-sized soil particles. The wind-blown depositional process of loess formation, however, promotes development of a relatively loose soil structure prone to changes in hydro-mechanical load conditions [18]. The silt particles are thus easily eroded by seepage or wind. Tabarsa et al. (2018) [19] investigated the effectiveness of the loess stabilisation using nanoclay both in the laboratory and in the field at the Gonbad dam site, considering various fractions of nanoclay ranging from 0.2% to 3% by mass. The field test section with 2% nanoclay showed the highest erosion resistance, while the laboratory specimens exhibited the same general trends in behaviour. The cost for the use of nanoclay, however, may be an issue as it is subjected to a limited budget, which also indicates restricted applicability. More and more natural hazards were associated with the loess soil in northwest of China because of its metastable structure and wetting-induced collapse deformation [20–25]. Most natural hazards were initiated by soil erosion caused by expansion of agricultural activities resulting from rapid population growth in the region [26,27]. This widespread engineering problem led to various scales of catastrophic slope sliding including shallow and deep landslides [28,29]. Recently, plant roots were used for preventing soil erosion [30–33] and stabilising slopes [34–36]. A similar idea was also applied to residential houses built with straw bale containing post-harvest waste (termed PHW hereafter), such as wheat straw and corn cob, and loess soil on the Chinese Loess Plateau. Field investigations have shown that the residential houses were utilised over decades under the effect of freeze-thaw cycles (Figure 1) [37]. The phenomena motivate this study to investigate how the addition of PHW to the slightly cemented loess improves its shear strength properties. The direct shear (DS) method has been deemed to be a quick and economic manner for estimating soil shear strength. Previous studies on direct shear have shown that shear boxes of different size may not lead to similar shear stress-horizontal displacement curve, but shear strength [38–40]. In spite of the fact that recently the effects of straw and biochar amendments on chemical activities in the loess plateau of China have been studied [41,42], to the authors' knowledge no studies have been performed to assess the significance of test size and PHW dosage to the shear strength of the slightly cemented loess soils. One study objective is to compare the shearing behaviour of the loess-PHW mixture specimens measured in large-scale and small-scale direct shear (LSDS and SSDS) tests with the loess specimens. The other objective is to determine the optimal dosage of PHW that matters most to the ability of the loess-PHW mixture to resist the seepage force and consequently stratum erosion.

**Figure 1.** Residential house built from loess-PHW mixture: (**a**) location of Chinese Loess Plateau, (**b**) house wall and (**c**) enlargement of house wall.

#### **2. Materials and Methods**

#### *2.1. Loess-Post-Harvest Waste (PHW) Mixture Specimens*

Approximately 1 m<sup>3</sup> of loess soil material from one sampling spot, Lantian, located in the Chinese Loess Plateau was retrieved for specimen preparation. The retrieved material passing the No. 4 sieve (4.75-mm opening) was used ensuring that the ratio of box length to maximum particle size was at least 10 and that the ratio of box thickness to maximum particle size was at least 6 [43]. The dry unit weight γ<sup>d</sup> of 13.72 kN/m<sup>3</sup> from the Standard Proctor test result was used for controlling compaction of specimens. The physical properties are summarised in Table 1. The particle-size distribution curves are shown in Figure 2. The material classifies as low plasticity silt (ML). The grading characteristics vary greatly with the uniformity coefficient Cu ranging from 9.26 to 15.26 although the majority of the material shares a common USCS (Unified Soil Classification System) designation. There were three block samples prepared at the designed water contents ω, i.e., 14%, 18% and 22%, respectively. Subsequently each sample was blended with the PHW treated with hot water (Figure 3) [44] and placed in a sealed container. The %PHWs by weight, while preparing the loess-PHW specimens, were equal to 0.3, 0.45, 0.6 and 0.75 respectively in SSDS and LSDS.


**Table 1.** Physical properties of tested loess soils.

Note: USCS represents the abbreviation of the Unified Soil Classification System.

**Figure 2.** Particle-size distribution curves for material used in direct shear tests.

**Figure 3.** Loess sampling and hot water treatment: (**a**) sampling location, (**b**) water bath apparatus, (**c**) water bath tank, (**d**) four sub-tanks and (**e**) post-harvest waste (PHW) before (right-hand side) and after (left-hand side) treatment.

#### *2.2. Small-Scale Direct Shear Tests*

We sheared 45 loess-PHW mixture specimens in the SSDS tests using a circular shear box of 61.8 mm in diameter and 20 mm in height. The SSDS tests were performed under normal pressures 100 kPa, 200 kPa and 300 kPa, respectively. The loess-PHW mixtures were compacted in the shear box in one single lift by tamping the top with a steel tamper. All the tests were conducted at a constant shearing rate of 0.8 mm/min to a maximum horizontal displacement (HD) of 7 mm. Measurements of HD and shear force were recorded by an industrial computer.

#### *2.3. Large-Scale Direct Shear Tests*

We sheared 27 loess-PHW mixture specimens in the LSDS tests using a square shear box (Figure 4). The shear box of 300 mm in width contains a 100-mm thick specimen. The loess-PHW mixtures were compacted following the same procedure as described for SSDS. Each mixture was sheared under normal pressures of 100 kPa, 200 kPa and 300 kPa, respectively. The applied normal pressures were manipulated by introducing a feedback-controlled pressure regulator. The prepared specimens were sheared at a rate of 0.8 mm/min to a maximum HD of 50 mm. Since the lower half of the shear box was bolted to the external box, a stepper motor that displaces the external box was used to control the HD. Measurements were also recorded for all the LSDS tests.

**Figure 4.** Schematic illustration of direct shear method: (**a**) large-scale direct shear (LSDS) apparatus and (**b**) enlargement of shear box.

#### *2.4. Repeatability of Direct Shear Methods*

Test repeatability has been deemed as the key to verify the effectiveness of the testing method adopted. Repeatability of the direct shear test was verified by performing six replicate tests. The six replicate tests on loess with %PHW equal to 0 were performed using the SSDS apparatus. The failure envelopes shown in Figure 5 from the replicate tests are similar. Shear strength parameters and coefficients of determination are summarised in Table 2. The φ' value varies from 16.3◦ to 22.3◦, with an average of 19.5◦ and a standard deviation of 2.2◦. The c value varies from 51.61 to 83.80 kPa, with an average of 67.94 kPa and a standard deviation of 11.93.

**Figure 5.** Results of repeatability tests for the small-scale direct shear (SSDS) test: (**a**) specimen at the water content of 14% and (**b**) specimen at the water content of 18%.


**Table 2.** Shear strength parameters for six replicate tests in SSDS.


**Table 2.** *Cont*.

#### **3. Results and Discussion**

#### *3.1. Definition of Failure*

Due to limited space, only some of the SSDS and LSDS results are presented. The HD in DS tests was normalized to specimen diameter, termed the relative horizontal displacement (RHD) hereafter, for the sake of easy comparison. Figure 6 shows the relationships between the shear stress and the RHD for the SSDS tests, while Figure 7 shows the relationships between the shear stress and the RHD for the LSDS tests. Generally, the shear stress versus the RHD relationships for the SSDS tests can be categorised into strain-softening curve (Type 1) and strain-hardening curve (Types 2 and 3) (Figure 8). Shear stress increasing to a peak value and then gradually decreasing with increasing horizontal displacement was categorised into Type 1 (Figure 6a). In the case shear stress increasing to an ultimate value and then remaining essentially constant with increasing horizontal displacement was categorised into Type 2 (Figure 6b). Type 3 comprised two sub-types. Type 3a can be described as shear stress increasing to the curve inflection point where some factor starts to have influence on the shearing behaviour and then gradually increasing at a constant rate with additional horizontal displacement. Type 3b was similar to Type 3a, but followed by a gradual increase in shear stress at another constant rate with larger horizontal displacement. The gradual increases in shear stress observed in Type 3 were most likely attributed to the effect of particle-box interaction, as discussed subsequently. It is evident that Type 1, Type 2 and Type 3a were characterised by two distinguished curve slopes, initial slope and additional-displacement slope and that Type 3b was characterised by three distinguished slopes; that are, initial slope, additional-displacement slope and larger-displacement slope. In Type 3a and Type 3b for which no distinct peak stress was available, a suggested criterion of 4 mm HD (i.e., 6.5% RHD) [44] to determine the failure stress still had some distance from the onset of the tangent line, which also indicated an inappropriate criterion for determining the failure stress. Thus, the failure stress for Type 1 and Type 2 was defined using the peak stress and the horizontal tangent line, respectively, and from the onset of the tangent line, resulting from the additional-displacement slope, for Type 3a and Type 3b.

**Figure 6.** *Cont*.

**Figure 6.** Shear stress-relative horizontal displacement relationships from SSDS tests on the loess-PHW specimens at ω = 14%: (**a**) %PHW = 0 and (**b**) %PHW = 0.3.

**Figure 7.** Shear stress-relative horizontal displacement relationships from LSDS tests on the loess-PHW specimens at ω=14%: (**a**) %PHW = 0 and (**b**) %PHW = 0.45.

**Figure 8.** Typical shear stress-horizontal displacement curves.

The shear stress versus the RHD curves from the LSDS tests exhibited two out of the four relationships observed from the SSDS tests, which are Type 2 (not shown) and Type 3a. The gradual increase in shear stress observed in Type 3a was believed to be related to the effect of particle-box interaction. Similarly, the failure stress for Type 2 and Type 3a was defined using the horizontal tangent line and the onset of the tangent line, respectively. Additionally, the suggested criterion of 4 mm HD was redefined by regarding the specimen size and applied to the shear stress versus the RHD curves for evaluation of its applicability. It was also found that the redefined criterion of 19.4 mm HD (i.e., 6.5% × 300 mm) differed considerably from the onset of the tangent line, indicating that the 19.4 mm horizontal displacement might not be appropriate for determining the failure stress.

#### *3.2. Comparison of Shear Strength between Small-Scale Direct Shear and Large-Scale Direct Shear*

A failure envelope of Mohr–Coulomb failure criterion was established through three DS tests that involved three different normal pressures using the linear least squares regression. The coefficient of determination R<sup>2</sup> for the failure envelopes in SSDS varied from 0.914 to 0.999, while in LSDS it varied from 0.906 to 0.999. Shear strength parameters for which the box friction and the possible non-linear nature of the failure envelopes near the origin were omitted are summarised in Table 3. Except the loess specimens with the %PHW equal to 0, the φ' values were typically 16.3–30.4◦ in SSDS and 15.3–29.3◦ in LSDS. In most cases the φ' value was the highest as the %PHW was equal to 0.60 and the effect of the PHW on the improvement of the φ' value was more significant in SSDS than in LSDS. A comparison of the φ' values from the SSDS tests and those from the LSDS tests is shown in Figure 9. It was noted that the difference between SSDS and LSDS was typically 0.7–8.8◦ except the difference of 17.8◦ for which the specimen prepared at ω = 22% with the %PHW being equal to 0. The PHW was deemed to be effective in achieving the effect of particle inter-locking and thus impeded the development of shear bands near the shear plane as subjected to shear force. The more PHW added, the greater the effect of particle inter-locking and the lesser the significance of the scale effect.

In addition to the φ' value, the RHD at failure was compared between SSDS and LSDS as well. The RHD varied within the range of 1.9–5.4% in SSDS, while it varied within the range of 0.2%–6.7% in LSDS (Figure 10). The difference in the RHD at failure between SSDS and LSDS averaged 1.6%, with a maximum of 4.8%. It is evident that scatter existed in the data although in most cases the difference was generally limited to less than 2.0%. On the other hand, the difference in the RHD in fact was found to be governed by the added water, not by the added amount of PHW. The effect of particle inter-locking restrained the development of shear bands in the vicinity of the shear plane

and became distinct when the water content was >14%. The lesser the amount of water added to loess-PHW mixture specimen, the larger the difference in the RHD and the greater the significance of the scale effect.


**Table 3.** Shear strength parameters of loess-PHW mixture specimens measured in direct shear tests.

**Figure 9.** Comparison of friction angle between SSDS and LSDS.

**Figure 10.** Comparison of relative horizontal displacement (RHD) at failure between SSDS and LSDS for three different water contents.

#### *3.3. Determination of Optimal PHW Dosage*

Additional LSDS tests were conducted for the sake of determining the optimal %PHW to be added, in which the %PHW of 1.0 was considered while preparing the loess-PHW mixture specimens at ω equal to 14%, 18% and 22%, respectively. The results of the LSDS tests are summarised in Table 4. The φ' value was the highest as the %PHW was 0.6. Higher φ' value was found as well at the %PHW of 1.0. It is noted that the φ' value for the specimens at ω = 14% increased from 22◦ to 29.3◦ as the %PHW increased from 0.45 to 0.6. Then the φ' value decreased to 20.3◦ as the %PHW further increased to 0.75. The φ' value, however, increased to 25.1◦ as the %PHW was ultimately increased to 1.0. Similar tendency could also be observed for the specimens at ω = 18% and 22%, respectively. It is evident that the optimal PHW dosage derived from this study was within a range of 0.55–0.65 considering specimen homogenisation and testing reliability. The main cause to as well lead to the increase in the φ' value at the %PHW > 0.8 was investigated by means of a simple test where the height of PHW (hPHW) in the shear box was measured under different normal pressures, hence allowing a back-analysis of the height of loess (hloess) in the shear box. Figure 11 shows the variation of hloess, γ<sup>d</sup> and φ' as the function of %PHWs added (i.e., 0.45, 0.6, 0.75 and 1.0) for the loess-PHW mixture specimens prepared at ω = 14%, 18% and 22%, respectively. The γ<sup>d</sup> values of loess were very close to their maximum of 15.6 kN/m3 from the Proctor curve (not shown) at %PHW of 0.8–1.0. This phenomenon was most likely ascribed to the densified loess resulting from a high %PHW added. It seemed that the shearing behaviour was governed by the amount of added PHW at %PHW <0.8 and that it, in turn, was governed by the densification effect as the %PHW was >0.8. To sum it up, the densification effect led to higher φ' value at %PHW of 0.8–1.0, and this study identified the optimal PHW dosage to be 0.55–0.65 taking into account specimen homogenisation and testing reliability. The impact of the amount of water added would be discussed later in this paper.


**Table 4.** Effect of PHW dosage on the improvement of friction angle φ' measured in LSDS tests.

**Figure 11.** *Cont*.

**Figure 11.** Variation of loess height, unit weight and friction angle from LSDS at %PHW of 0.45, 0.6, 0.75 and 1.0: (**a**) specimen at ω = 14%, (**b**) specimen at ω = 18% and (**c**) specimen at ω = 22%.

#### *3.4. E*ff*ect of Particle–Box Interaction*

As discussed, the gradual increase in shear stress following the curve inflection point was observed both in Type 3a and Type 3b and its main cause was still not clear. This phenomenon was further studied and elaborated upon. In this regard, two more displacement transducers attached to front and rear of the shear box were introduced for assessing particle movement while shearing in LSDS mode. Figure 12 shows the vertical displacement at front and rear of the shear box versus the RHD for the loess-PHW mixture specimen at ω = 18% and 22%, respectively, as subjected to the normal pressures of 100, 200 and 300 kPa. The vertical upwards displacement at front of the shear box was at most 5.3 mm, while the vertical downwards displacement at rear was typically 6–12 mm. The magnitude of the vertical displacements was more distinct at the rear of the shear box than at the front. Additionally, the greater the normal pressure applied, the larger the vertical displacement. The above results indicate that the particles, during shearing, tended to move to front and rear of the shear box and that the greater normal pressure largely promoted the movement of particles. The observed movement of particles was deemed as the main cause to lead to the further gradual increase in the shear stress observed both in Type 3a and Type 3b.

**Figure 12.** Vertical displacement at front and rear of the shear box versus relative horizontal displacement (RHD) considering three normal pressures of 100, 200 and 300 kPa: (**a**) specimen at ω = 18% and (**b**) specimen at ω = 22%.

#### *3.5. Mechanism Leading to Improved Shearing Behaviour*

Addition of PHW hampered particle dislocation due to particle inter-locking. This effect thus impeded the development of shear bands in the vicinity of the shearing plane as the HD increased. The difference in the φ' value between SSDS and LSDS was decreased for this reason. It was also found that the RHD rather was a function of the water content than of the amount of PHW added. In spite that addition of PHW hampered particle dislocation, this could not take effect with the small amount of water contained in the specimen. A more accurate way to describe the mechanism leading to the improved shearing behaviour of the loess-PHW mixture is that the improvement was attributed to the effect of particle inter-locking resulting from the added PHW and this effect took effect as the water content reached a threshold, i.e., >14%, that facilitates particle rearrangement. The findings prove not only the reduced environmental impacts while introducing the loess-PHW mixture as compared with other ordinary materials, but also a great potential of its application to sustainable development of urban areas [45–47].

#### **4. Conclusions**

This study investigated the shearing behaviour of the loess-PHW mixture specimens with various water contents and different PHW dosages in SSDS and LSDS. Some main conclusions can be drawn as follows:


**Author Contributions:** This paper represents a result of collaborative teamwork. Z.-F.X. and L.W. developed the concept and prepared the manuscript; W.-C.C. and J.X. provided constructive suggestions and revised the manuscript. The four authors contributed equally to this work.

**Funding:** This research received no external funding.

**Acknowledgments:** This work would not have been possible without technical supports from Dr. Zhao Duan at Xi'an University of Science and Technology.

**Data Availability:** The laboratory data used to support the findings of this research work are included within the article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Biostimulation in Desert Soils for Microbial-Induced Calcite Precipitation**

#### **Hadas Raveh-Amit <sup>1</sup> and Michael Tsesarsky 2,3,\***


Received: 4 March 2020; Accepted: 19 April 2020; Published: 23 April 2020

**Abstract:** Microbial-induced calcite precipitation (MICP) is a soil amelioration technique aiming to mitigate different environmental and engineering concerns, including desertification, soil erosion, and soil liquefaction, among others. The hydrolysis of urea, catalyzed by the microbial enzyme urease, is considered the most efficient microbial pathway for MICP. Biostimulated MICP relies on the enhancement of indigenous urea-hydrolyzing bacteria by providing an appropriate enrichment and precipitation medium, as opposed to bioaugmentation, which requires introducing large volumes of exogenous bacterial cultures into the treated soil along with a growth and precipitation medium. Biostimulated MICP in desert soils is challenging as the total carbon content and the bacterial abundance are considerably low. In this study, we examined the biostimulation potential in soils from the Negev Desert, Israel, for the purpose of mitigation of topsoil erosion in arid environments. Incubating soil samples in urea and enrichment media demonstrated effective urea hydrolysis leading to pH increase, which is necessary for calcite precipitation. Biostimulation rates were found to increase with concentrations of energy (carbon) source in the stimulation media, reaching its maximal levels within 3 to 6 days. Following stimulation, calcium carbonate precipitation was induced by spiking stimulated bacteria in precipitation (CaCl2 enriched) media. The results of our research demonstrate that biostimulated MICP is feasible in the low-carbon, mineral soils of the northern Negev Desert in Israel.

**Keywords:** microbial-induced calcite precipitation; desert soil; biostimulation; erosion mitigation

#### **1. Introduction**

The ubiquity of bacteria and the diverse roles they play in natural environments have led to growing interest in harnessing bacterial activities for various anthropogenic purposes. From a physical point of view, soil is regarded as an inorganic multiphase system comprising solids, fluids, and gases. However, soil is also a living system, being one of the largest terrestrial carbon pools, constituting about 33% of the total terrestrial carbon [1]. The organic carbon in the top 1 m constitutes more than 50% of the total soil carbon. Prokaryotes comprise up to 17% of the soil organic carbon [2]. These unicellular organisms, mostly bacteria 0.5–5.0 <sup>×</sup> 10−<sup>6</sup> m in size, are about three orders of magnitude smaller than the pore throat size of sand and about the D10 size of kaolinite [3]. Soil bacteria, either motile or fixed to mineral surfaces (grains), may change the chemical and physical properties of their surroundings depending on their metabolism. Microbial biomass and biodiversity often exhibit exponential decreases with depth [4]; nevertheless, there are still active cells in deeper soil horizons [5].

Many bacteria are capable of inducing mineral precipitation through various metabolic paths in both oxic and anoxic environments. Boquet [6] showed that most heterotrophic bacteria can induce precipitation of calcium carbonate (CaCO3), a common natural cementing agent, by various metabolic pathways.

Microbial-induced calcite precipitation (MICP) is an emerging technique aiming to mitigate different environmental and engineering challenges, including soil erosion, soil liquefaction, fracture sealing, restoration of stone monuments, and others [7–13]. Hydrolysis of urea, catalyzed by the microbial enzyme urease (urea amidohydrolase, EC 3.5.1.5), is considered the most efficient microbial pathway for MICP [14,15]. The hydrolysis of urea produces ammonium and carbonate, seen in Equation (1), thus increasing the saturation for calcium carbonate, which could lead to its precipitation, usually as calcite, in Equation (2).

$$2\text{CO}(\text{NH}\_2)\_2 + 2\text{H}\_2\text{O} \overset{\text{urease}}{\rightarrow} 2\text{NH}\_4^+ + \text{CO}\_3^{2-} \tag{1}$$

$$\rm{CaCO\_3^{2-}} + \rm{Ca^{2+}} \leftrightarrow \rm{CaCO\_3} \downarrow \tag{2}$$

Two alternative approaches to facilitate in situ MICP in soils are biostimulation and bioaugmentation. Biostimulation encourages indigenous urea-hydrolyzing bacteria by providing appropriate enrichment and precipitation media; it relies on the natural ubiquity of ureolytic soil bacteria and bacterial spatial distribution [15,16]. Bioaugmentation introduces large volumes of bacterial cultures into the treated soil along with a growth and precipitation medium; therefore, it requires large volumes of pure-cultured specializing ureolytic bacteria, e.g., *Sporosarcina pasteurii*. Producing and transporting large volumes of these cultures is an expensive and delicate procedure; their injection and homogeneous distribution throughout the treated site are difficult to achieve and might encounter regulatory hindrances [17]. Moreover, the introduced bacteria are likely to decline in numbers due to low compatibility with the environment as well as competition and predation by indigenous bacteria [18]. As the two methods involve the introduction of large quantities of urea along with a source of organic carbon, both biostimulation and bioaugmentation are likely to affect and change the indigenous microbial population.

MICP soil improvement was successfully demonstrated at a variety of scales [10,19–22]. However, reliance on bioaugmentation has restricted the technology from becoming a cost-competitive alternative to more traditional ground amelioration techniques. Biostimulation, the use of selective substrates and environmental factors to stimulate the growth of native microorganisms with desirable metabolic capabilities, has been researched extensively in the field of bioremediation [23,24] with success in several notable field-scale applications [25,26]. Despite the frequent use of biostimulation in the field of bioremediation, the use of the biostimulation technique for enabling MICP is more limited [5,27–32].

Biostimulated MICP in desert soils is more challenging than in most soils as the total carbon content is low [33], and the bacterial population is considerably smaller. For example, in coastal sands from southern Israel (31.61◦N, 34.50◦E), the in situ bacterial population was found to be in the order of 10<sup>4</sup> cells/g [31], three to four orders of magnitude lower than in semiarid soils [34]. For this reason, MICP in low-carbon soils were typically attempted via bioaugmentation [35,36].

In this paper, we present the results of a research aimed to study the biostimulation potential in desert soil from the northern Negev Desert, Israel. The results reported here are part of an ongoing research project aimed to develop a MICP-based technique for the mitigation of topsoil erosion in arid environments and is planned to be followed by laboratory and field experiments. We envision different applications for MICP-based topsoil stabilization, such as the suppression of dust emission from vehicles in construction sites and open-pit mines, improving the efficiency of solar energy facilities and the mitigation of erosion by flush flood, among others.

Biostimulation for ureolytic MICP was confirmed at laboratory conditions by incubating soil samples in urea and enrichment media and monitoring urea depletion and pH evolution within a few days. Precipitation experiments using the indigenous, biostimulated bacteria demonstrated that calcite precipitation could be readily induced. Altogether, we report that urea-hydrolyzing bacteria are naturally present in the topsoil sampled from the northern Negev Desert and that these bacteria can be effectively stimulated to induce calcite precipitation.

#### **2. Materials and Methods**

#### *2.1. Soil Sampling and Chemical–Physical Characterization*

Soils were sampled from a depth of 0.4 m from two sites (soil 1 and soil 2 separated by 500 m) of the Rotem Plateau (31.04◦N, 35.08◦E) at the northern Negev Desert in Israel, an arid region with an average annual rainfall of 70 mm. Samples were stored refrigerated at 4 ◦C until the biostimulation experiments began. Soils were analyzed for elemental composition by X-ray fluorescence (XRF) using an EX-Calibur spectrometer (Xenemetrix, Migdal HaEmek, Israel). Mineralogical phase identification was performed by X-ray diffraction (XRD) using a Bruker D8 Advance system (Bruker, Billerica, MA, USA). Particle size distribution (PSD) was performed by laser diffraction using a Mastersizer 3000 system (Malvern Panalytical, Malvern, UK).

#### *2.2. Biostimulation of Indigenous Ureolytic Microbes*

Biostimulation was performed by incubating 10.0 g of soil samples collected from each of the two sites from the Rotem Plateau in 100 mL of different stimulation media at ambient temperature with gentle shaking at 100 rpm for ten days. The media composition is described in Section 2.3. Samples were taken periodically for chemical analysis, as described in Section 2.4.

#### *2.3. Solutions and Stimulation Media*

Four different stimulation media used in this study were prepared in artificial groundwater (AGW) consisting of MgCl2 (1 mM), MgSO4 (1 mM), NaHCO3 (2.56 mM), NaCl (14.35 mM), CaCl2 (2.43 mM), and KCl (0.32 mM); ionic strength was at 31.5 mM and pH at 7.7 [31]. Urea control media consisted of 330 mM urea and served as a control for bacterial stimulation without an energy source. Two stimulation media contained a 330 mM urea supplemented with either a low dose (0.1 g/L) or high dose (1 g/L) of yeast extract. AGW only media was used as a negative control.

#### *2.4. Chemical Analysis of Stimulated Samples*

During the bacterial stimulation experiments, samples were taken periodically for pH and urea concentration measurements. pH was measured using a Metrohm pH-meter (Metrohm, Herisau, Switzerland). Urea concentrations were measured according to the Knorst [37] colorimetric method, with minor modifications, on an 8453 Agilent spectrophotometer (Agilent, Santa Clara, CA, USA). During the stimulation experiment, urea concentration values of the urea control media were elevated by a 2.5% increase per day due to media evaporation. Urea concentrations of all urea-containing samples were therefore normalized accordingly.

#### *2.5. Calcite Precipitation*

Biostimulated bacteria (100 μL aliquots taken from stimulation supernatants) were spiked into 100 mL of yeast extract and CaCl2-containing enrichment media at ambient temperature with gentle shaking to allow for calcite precipitation. Urea concentrations were systematically measured, as described in Section 2.4. At the end of the experiment (12–14 days), filtered and dried precipitates were analyzed by XRD.

#### **3. Results**

#### *3.1. Chemical and Physical Characterization of Soils*

Soils were collected from two sites of the Rotem Plateau at the northern Negev Desert in Israel, a broad plateau covered with clastic sediments of the Miocene Hazeva Fm. Soils are mainly composed of quartz and calcite, as shown by XRD analyses (Figure 1). XRF analyses showed that the two soils are similar in their elemental composition; however, significantly higher levels of Si and lower levels of Ca were identified in soil 1 than soil 2 (Table 1). Particle size distribution analysis showed that the soils are classified as sand (Figure 2 and Table 2), with soil 2 showing a higher coarse sand fraction 94%, compared to 64% in soil 1, and a lower fine sand fraction, 6% and 36% in soil 2 and soil 1, respectively. μ μ

**Figure 1.** X-ray diffraction patterns of soil 1 sampled from the Rotem Plateau. Counts in arbitrary units.



**Figure 2.** Particle size distribution of soils sampled from the Rotem Plateau in Israel.

**Table 2.** Particle size distribution (by volume) properties of soils from the Rotem Plateau.


#### *3.2. Biostimulation of Indigenous Urea-Hydrolyzing Bacteria*

Biostimulation experiments were conducted on soils sampled from the Rotem Plateau to test whether indigenous urea-hydrolyzing bacteria are naturally present in this region and can be effectively stimulated. Soils were incubated in different media for ten days, and urea hydrolysis was tested by pH and urea concentration measurements (Figure 3). Based on our previous experience [31], it was expected that if indigenous urea-hydrolyzing bacteria would be present in the soil, then their stimulation in the presence of urea would result in a pH rise due to the production and accumulation of ammonium ions. Stimulation of urea-hydrolyzing bacteria was evident in soils incubated in the presence of urea and a carbon source (yeast extract), as shown by pH elevation and urea degradation (low yeast extract concentrations (YE\_L) and high yeast extract concentrations (YE\_H)). The addition of a carbon source was required to achieve the stimulation, as no significant ureolysis was observed in the absence of yeast extract (Urea). Low yeast extract concentrations (YE\_L) resulted in inefficient stimulation when compared to high yeast extract concentrations (YE\_H) in the two soils studied. As expected, the pH of the negative control samples from the two sites remained unchanged, and urea concentrations were negligible throughout the biostimulation experiment (Figure 3, AGW).

**Figure 3.** Biostimulation of soils from the stimulated soils from the Rotem Plateau. (**a**) pH of soil 1; (**b**) pH of soil 2; (**c**) urea concentration of soil 1 and (**d**) urea concentration of soil 2. AGW is artificial ground water (o); urea is urea control (∇); YE\_L is low (0.1 g/L) yeast extract enriched medium (Δ); YE\_H is high (1.0 g/L) yeast extract enriched medium (♦).

To provide direct evidence that MICP can be induced in these soils, precipitation experiments were conducted using bacteria stimulated with urea and 1.0 g/L yeast extract concentrations. Aliquots of bacteria stimulated from soils for seven days were spiked into fresh urea, and yeast extract-containing media were supplemented with CaCl2. Within two days, the media became turbid, and white precipitates appeared at the bottom of the Erlenmeyer flasks (data not shown). Urea hydrolysis was confirmed in samples provided with stimulated bacteria from either soil, albeit higher rates were observed in soil 2 than in soil 1 (Figure 4). This effect probably resulted from differences in biostimulation rates between the two soils (Figure 3). After two weeks of incubation, precipitates were collected, dried, and identified by XRD as calcite (Figure 5). Neither precipitation nor urea hydrolysis took place in negative control samples (AGW) incubated in the same conditions.

**Figure 4.** Urea concentrations measured in precipitation experiments following bacterial stimulation in soils from the Rotem Plateau.

**Figure 5.** XRD patterns of precipitates produced by biostimulated soil (1) bacteria grown in CaCl2 precipitation media. Counts in arbitrary units. CaCO3 is calcite reference RRUFF R040070.

#### **4. Discussion**

Soil erosion, carbon sequestration, infrastructure rehabilitation, hazardous waste disposal, and water resources protection are all 21st-century global challenges. Many of these challenges occur within or are supported by soil. The conventional perspective views soil as an infinite resource, composed of discrete functions, e.g., hydraulic, mechanical, and ecological, among others. However, soil is also a diverse ecosystem. The biotic potential in soil offers prospects for innovative and sustainable solutions for some of these challenges. Harnessing natural biogeochemical processes to improve the environmental conditions and/or engineering properties of geological deposits has received significant attention in the scientific community [5,22,36,38,39].

Microbial-induced calcite precipitation (MICP) is one of the most promising biogeochemical treatments that effectively changes the hydromechanical and environmental properties of materials. Potential applications range from groundwater remediation and sequestration of radionuclides [40] to self-healing concrete [41]. For engineered materials, such as concrete, MICP can be achieved via bioaugmentation only, i.e., the introduction of monoclonal cultivated bacteria. In natural environments, such as soils, effective bioaugmentation by exogenous bacteria is not warranted, as the introduced bacteria are likely to decline in numbers due to low compatibility to the environment, as well as competition and predation by indigenous bacteria [18]. Moreover, producing and transporting large volumes of these cultures is an expensive and delicate procedure; their injection and homogeneous distribution throughout the treated site are difficult to achieve and might encounter regulatory hindrances.

The ability to hydrolyze urea is widely distributed among indigenous bacteria in soils [15]. Hence, the biostimulation approach of MICP in soils is tangible. The main challenge with biostimulated MICP is inducing effective urea hydrolysis, which requires a sufficiently large ureolytic population. This challenge is even greater in nutrient-poor soils, where the initial bacterial abundance is low.

Gat [31] showed that for MICP, effective biostimulation of an indigenous, ureolytic population in soil from a semiarid region (BSh (hot semiarid climate) in Köppen climate classification) requires an energy (carbon) source, in addition to urea. No ureolysis was observed in the absence of a carbon source. The authors showed significant changes to the indigenous bacterial population following biostimulation, where ureolytic bacteria population increased from 5% in the native sand up to 99% in carbon high dosage treatments.

The primary goal of our research was to evaluate the potential of biostimulating urea-hydrolyzing bacteria in desert soils of the Rotem Plateau of the northern Negev Desert, Israel (BWh (hot desert climate) in Köppen climate classification), for soil erosion mitigation via MICP biocementation of the topsoil. We found that indigenous urea-hydrolyzing bacteria are naturally present in the desert soils and can be readily stimulated to achieve effective urea hydrolysis and calcite precipitation within days. These results are consistent with the results of [31]. It was also found that bacteria stimulation was not affected by differences in soil mineralogy observed between the two soils. Gomez [5] showed that urea hydrolysis could be effectively stimulated even at greater depths (10 m) using low doses of yeast extract and alkaline pH adjustment of the treatment media. Biostimulation rates in the desert soils were accelerated by providing a yeast extract dose of 1.0 g/L, whereas a lower dose of 0.1 g/L was found to be ineffective. Other more cost-effective energy sources can be utilized to promote biostimulation, such as off-the-shelf molasses [31] with similar concentrations of 1.0 g/L.

To conclude, the results presented above demonstrate that desert soils from the Rotem Plateau of the northern Negev Desert (Israel) are susceptible to the biostimulation of ureolytic native bacteria. Effective urea degradation is a necessary requirement for in situ MICP, which is designed to achieve effective biocementation to mitigate topsoil erosion.

**Author Contributions:** H.R.-A. and M.T. designed the research methodology; H.R.-A. performed the physiochemical analyses; H.R.-A. and M.T. analyzed the data and wrote the original draft; H.R.-A. and M.T. secured the research funding. The authors have read and agreed to the published version of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the Pazy Foundation, grant number 18-2018.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Lal, R. Sequestration of atmospheric CO2 in global carbon pools. *Energy Environ. Sci.* **2008**, *1*, 86–100. [CrossRef]


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Comparison of Diverse Dust Control Products in Wind-Induced Dust Emission from Unpaved Roads**

#### **Itzhak Katra**

Department of Geography and Environmental Development, Ben Gurion University, Beersheba 8410501, Israel; katra@bgu.ac.il

Received: 15 November 2019; Accepted: 27 November 2019; Published: 29 November 2019 -

**Abstract:** Surfaces of unpaved roads are subjected to dust PM10 (particulate matter < 10 μm) emission by wind process regardless of vehicles (wheels) transport. However, there is little quantitative information on wind-induced dust emission from unpaved roads and the efficiency of diverse dust control products. The study aimed to fill this clear applied scientific gap using wind-tunnel experiments under laboratory and field conditions. The wind-tunnel complies with aerodynamics requirements and is adjusted to dynamic similitude by appropriately scaling all variables that affect dust transport. The results of the control sample (no-treatment) clearly show that dust emission by wind from unpaved road could be a substantial contribution to mass transfer and air pollution, and thus should be considered. Diverse dust control products of synthetic and organic polymers (Lignin, Resin, Bitumen, PVA, Brine) were tested. In the first stage, the products were tested under controlled-laboratory conditions. The results enabled quantitative assessment of the product efficacy in wind erosion without the impact of vehicle transport. In the second stage, the products were tested in field experiment in an active quarry, in which the products were applied on plots along the road. The field experiment was conducted after transportation of the quarry-haul trucks in two time points: several days after the application, and several weeks after the application. The results show that in most of the plots the dust emission increases with the wind velocity. PM10 fluxes from the road surface in each plot were calculated to determine the effectiveness of the dust control products. Some products significantly reduced dust emission from quarry roads, especially when using the Hydrous magnesium chloride (Brine). Additional experiments revealed that such Brine can be applied with reduced amounts and still keeping on low emission.

**Keywords:** wind erosion; dust; suppressants; PM10; wind tunnel; lignin; resin; bitumen; PVA; brine

#### **1. Introduction**

Dust emission has significant implications on the earth's systems including atmosphere [1], marine ecosystems [2], and terrestrial ecosystems [3]. Moreover, dust emission is significant for air pollution and human health. The atmospheric PM pollution can increase high above the standard values of air quality [4,5]. During dust storms, average daily PM10 concentrations can be about 40 times higher than the guideline of WHO (World Health Organization) with PM10 concentration of >2000 μg m<sup>−</sup>3. The impact of dust-PM10 on public health is related to asthma in children [6], cardiovascular morbidity [7], and more. In recent years, there is an increasing interest in quantifying PM emission from dust sources due to anthropogenic activities from environmental quality and human health perspectives.

Anthropogenic activities in soils such as quarrying, mining, and off-road vehicles traveling, are major concern for dust emission by wind. The U.S. Environmental Protection Agency (EPA) published a Compilation of Air Pollutant Emission Factors (AP-42) with reference to dust emission in unpaved roads because of vehicle travels [8]. In unpaved roads, dust is released by direct uplift because of the contact between the wheels and the surface, and indirect uplift is caused by turbulence generated by the vehicle [9]. An experimental study showed a linear relationship between unpaved road dust PM10 emissions and vehicle speed by testing re-entrained aerosol kinetic emissions from roads [10]. The PM10 entrained into the atmosphere and transported over relevant distances should be calculated based on the fluxes of dust emissions [11]. A numerical model was formulated using dust emission parametrizations to present the atmospheric PM dispersion from unpaved roads [12].

Regardless of the dust emission by vehicle travels, unpaved roads are subjected to dust emission by wind process, and overall these could be a substantial contribution to mass transfer and air pollution. Dust emission by wind depends on the surface characteristics to determine the threshold of wind velocity that will initiate the transport of soil particles [13,14]. The fine-dust particles are usually held in soil aggregates and thereby are not available for direct lifting at low wind velocities. In this case the emission of the fine particles may be enabled by higher shear velocity or by saltation of sand particles and/or aggregates [15,16]. The soil aggregates can breakdown to release dust particles either by self-saltating or by the impact of other saltators [13,14]. Vehicles on unpaved roads alter the surface properties and may amplify dust emission by wind. Pulverization of the top layer by the tires facilitates the emission as it breaks soil aggregates [17]. The disintegration of the soil increases the amount of fine particles in the surface (dust-sized particles) that are subjected to direct lifting at lower wind velocities [14] with expected higher dust fluxes to the atmosphere.

Dust emission can be reduced by using dust control products. Such products are based mainly on synthetic or natural polymers, which applied by wetting the soil surfaces with the products. Note that this application refers to "dust control" in this study rather than "soil stabilization" that is typically associated with application of mixing (and packing) the soil layers with products to improve engineering properties. A wide range of dust control products have been tested for dust emission by vehicles traveling in unpaved roads [18,19]. There is lack of fundamental research examining the efficiency of diverse products in suppression of dust emission by wind. This study aimed to fill this gap by targeted wind tunnel experiments under laboratory and field conditions.

#### **2. Materials and Methods**

#### *2.1. Laboratory Experiment*

A soil from a typical calcareous quarry (Hartuv, Beit Shemesh, Israel), which produces calcareous materials for construction, was used for the experiment. This soil-material is a product of rock crushing and screening. It uses for bedding of the unpaved roads in the quarry. A soil sample was transported to the laboratory for experiment. The soil sample characterized by fine-earth material and rock fragments up to size of 50,000 μm in diameter (Figure 1a). The fine-earth material in the soil sample (after extraction of the rock fragments in the sample), which is potential of dust emission, was analyzed for particle size distribution (Figure 1b). It was done by the ANALYSETTE 22 MicroTec Plus (Fritsch, Idar-Oberstein, Germany) laser diffraction, which measures particles in the size range of 0.08–2000 μm [15,20]. Replicas (100 mg) were dispersed in Na-hexamvetaphosphate solution (0.5%) by sonication (38 kHz). The data was calculated using the Fraunhofer diffraction model with a size resolution of 1 μm using MasControl software. Mineralogical and elemental compositions show high content of Ca-based components (Figure 1c,d), which obtained by the X-ray power diffraction (XRPD) and X-ray fluorescence (XRF) methods, respectively [21].

**Figure 1.** Characteristics of the soil sample used for the experiments: (**a**) size distribution of the soil product in the quarry obtained by mechanical sieving. (**b**) Particle size distribution of the fine-earth material (after extraction of the rock fragments in the sample) obtained by the laser diffractometer technique. (**c**) Mineralogical composition of the fine-earth material by the X-ray power diffraction (XRPD). (**d**) Elemental composition of the fine-earth material by the X-ray fluorescence (XRF).

The soil that was brought to the lab from the quarry (Figure 1a) was sieved into size material size <8000 μm in diameter, including fine-earth materials (<2000 μm) and small rock fragments (2000–8000 μm). The soil was placed in specific trays that fit the wind tunnel dimensions (see next section). The area of each tray is 0.5 × 1.0 m with height of 0.02 m. Three trays (replicas) were prepared for each soil treatment/product (Figure 2). Table 1 presents the various products and their application on the soil used for the laboratory experiment. The composition of the solution of each substance and the concentration of the solution to be applied by spraying on the soil surface were provided by the suppliers of the products. Note that the brine was supplied as a solution already, in which the concentration applied on the soil (1.5 L m<sup>−</sup>2) is typical in dust control application of Brine. Following this process, the amount of the solution (water content) was vastly different between the soil samples/applications. However, the soil samples left outside for a drying process of several days before the experiment to avoid the impact of soil water content on wind erosion and dust emission.


**Table 1.** List of products for dust control used in the laboratory and field experiments. Note that the composition of the solution of each substance and the concentrations of the solution to be applied by spraying on the soil surface were provided by the suppliers.

The soils were tested in a wind tunnel for dust emission—overall 108 runs were performed by the wind tunnel (9 samples, 3 replicates, 4 wind velocities). These experiments were performed in the portable boundary-layer wind tunnel of the Aeolian Simulation Laboratory at BGU [3,22]. Boundary-layer wind tunnels enable simulations under standardized quasi-natural wind conditions and provide quantitative information on dust emission rates from the soils. The wind tunnel has a cross sectional area in the order of 0.5 × 0.5 m with open-floored working sections of up to 10 m length. Air push or air suction flow in the tunnel is generated by an axial fan up to maximum velocity of 18 m s<sup>−</sup>1. Instruments installed in the test section of the tunnel enable quantification of wind characteristics and dust (PM) concentrations.

The tunnel operated under four wind velocities, 4.0, 6.5, 8.0, and 9.5 m s<sup>−</sup>1, representing a range of common natural winds that are associated with dust emission [13]. Note that the frequency of a specific wind velocity should be considered for calculation of the total dust emission in a period of time (month, year). Each run of specific wind velocity lasted 30 s to record the trend of dust emission. Dust concentrations (μg m−3) PM10 were recorded by a light-scattering device, DustTrak DRX 8534 (ww.tsi.com), in the range of 0.001–150 <sup>μ</sup>g m−<sup>3</sup> (±0.1% of reading) at 1-s intervals, was placed at 25 cm above the tunnel bed. The recorded PM10 concentrations were converted into mass flux emitted from the soil surface (mg m−<sup>2</sup> s−1) based on the wind tunnel dimensions, the area of the soil bed, and the wind velocity (for in depth description refer to Katra et al. 2016 [23]). The PM10 flux from the soil surface can be used in parameterization of dust emission and transport models [12,23].

**Figure 2.** The soil samples in trays after the application of the dust control product used for the laboratory wind tunnel experiment. The size of each tray is 0.5 × 1.0 m with height of 0.02 m. each soil sample was prepared in three replicas.

#### *2.2. Field Experiment*

The field experiments aimed to test the efficiency of the dust control products under quarry conditions over time. The field experiments were conducted in an active quarry near the city of Beit Shemesh, Israel (Figure 3). A typical un unpaved haul road (~200 m long) of the quarry was used for the experiment of dust control. The road is used for the Euclid Trucks, which ride the road dozen times during the day, with empty (~100 ton) or full (~170 ton) container. The dust control products were applied on equal plots (20 m each) along the road. Each plot used for different product (Table 1). The product list for the field experiment included Brine. Lignin, Resin, Bitumen, PVA. The products were prepared and applied in the same procedure to the laboratory experiment (Section 2.1, Table 1). The solutions were sprayed on the soil using water-truck with taps at the back (Figure 4).

**Figure 3.** The active quarry near the city of Beit Shemesh, Israel (31◦45 54.62" N, 35◦1 5.15" E) is shown at the top right side of the aerial photo. The field experiment was conducted on an unpaved haul road (~200 m long) of the quarry, which is within the yellow rectangular.

The road was closed for two days to allow the soil to dry (the experiment performed during July with average daily temperature of 30 ◦C). The experiment was conducted after transportation of the quarry-haul trucks in two time points: several days after the application of the products, in which 150 rides of Euclid Trucks were recorded, and two weeks after the application, in which additional 296 rides (total of 446) of Euclid Trucks were recorded. The measurements of dust emission were obtained by the portable wind tunnel, which was installed on the road at each soil plot/treatment (in three replicas for each plot). The operation procedure of the wind tunnel was the same to that of the laboratory experiment (see last paragraph in Section 2.1).

**Figure 4.** Application of the dust control product on the unpaved haul road before the dust emission experiment.

#### **3. Results and Discussion**

#### *3.1. Product E*ffi*ciency under Laboratory Conditions*

The wind tunnel experiments were aimed to test dust emission by wind after the applications of the dust control products. The PM10 concentration measured in the wind tunnel during the experiment is a result of dust emission from the soil surface. The results of the laboratory wind tunnel experiment are presented in Figure 5. Each soil sample was tested under four wind velocities. It is clearly shown that no dust emission and PM10 concertation was recorded in the application of Brine. The "straight" line that was similarly recorded for all wind velocities in the test of the brine application is a representation of the background level (~20 μg m−3) measured in the wind tunnel before the test. In all the other applications, the averaged PM10 concentrations were significantly

(*p* < 0.05) higher than the background level, including in the lowest wind velocities, 4 and 6.5 m s<sup>−</sup>1, of the PVA test. At wind velocity of 4 m s<sup>−</sup>1, which represents a wind at which distinct dust emission started from loess soils, the resulted PM10 concentration was the lowest for all samples. However, no significant difference (*p* < 0.05) in PM10 between 4 and 6.5 m s−<sup>1</sup> was received for Resin, Lignin, and PVA applications. The PM10 values of the control sample were significantly higher than those of all the other samples in all wind velocities. It means that all the applications suppressed dust emission but with different efficiency.

As mentioned before, the Brine was the best application, preventing any dust emission. The PVA was the second-best application, in which dust emission was significant only at the highest wind velocities (8 and 9.5 m s<sup>−</sup>1). Basically, the highest value of PM10 concentration in all soil samples was recorded in wind of 9.5 m s−1. However, the PM10 concentrations of the PVA at the highest wind velocity were significantly lower (*p* < 0.05) than those of the Resin, Bitumen, and Lignin. At wind velocity of 8 m s<sup>−</sup>1, no significant difference in PM10 concentration was received between the PVA and the organic polymers (Resin, Lignin). Note that because of the experimental procedure, in which the wind velocities of 6.5, 8, and 9.5 m s−<sup>1</sup> were run promptly after the prior velocity (i.e., 6.5 m s−<sup>1</sup> was run after 4 m s<sup>−</sup>1, and so on), and since the amount of dust particles in the soil samples is limited, the dust emission in those velocities could have been an underestimation for non-limited samples. Nonetheless, dust emission measurements in the laboratory present higher PM10 concentrations, compared to field measurements, under the same experimental conditions (soil type and wind velocity), because the relatively high disturbance of the soil structure during the preparation of the experimental samples in the laboratory [12,14].

**Figure 5.** PM10 concentrations recorded during the laboratory wind tunnel experiment on the diverse dust control products under wind velocities 4.0, 6.5, 8.0, and 9.5 m s<sup>−</sup>1. Note, the "control" is a soil sample with no treatment. The background level measured in the wind tunnel before the test was ~20 μg m<sup>−</sup>3.

The recorded PM10 concentrations were used to calculate the average dust-PM fluxes emitted from the soils (Table 2). The calculation was based on the wind tunnel dimensions and the area of the soil bed [23]. Using the mass flux can provide us quantitative information on soil loss by dust emission from the surface as well as parameterization in models of dust emission and atmospheric PM transport. Since the averaged dust fluxes presented in Table 2 were derived from the PM10 concentrations (Figure 5), they are in full correlation with the averaged PM10 concentrations and have the same

differences (and significance level, *p* < 0.05) discussed above. This is relevant also for the convenience of presentation of the dust suppression efficacy of the diverse products presented in the lower part of Table 2. The dust suppression efficiency was derived as the ratio between the averaged dust flux (per wind velocity) of specific application (Brine, Lignin, Resin, Bitumen, PVA) and the averaged dust flux of the control sample (no-treatment). The values are presented in percentage. It clearly shown that the Brine is the most efficient product in dust suppression, reducing more than 95% of the PM10 dust emission. Nonetheless, the minimum efficiency was 77.7% (Bitumen at 6.5 m s<sup>−</sup>1) among all the products and wind velocities under the controlled conditions of the laboratory experiment.


**Table 2.** The calculated PM10 emission fluxes from the soils tested in the laboratory wind tunnel experiment are presented in the upper part of the table. In the lower part of the table, the dust suppression efficiency (%) calculated as a ratio in dust flux between the treated road section and the control (no-treatment) road.

#### *3.2. Product E*ffi*ciency under Field Conditions*

After the laboratory experiment and first characterization of the product efficiency in dust emission suppression, the products were applied in an unpaved road of the quarry for field experiment on dust emission by the portable wind tunnel in two stages. In the first stage, the experiment was conducted several days after the application of the products, in which 150 rides of Euclid Trucks were recorded. The results of the PM10 dust emission from the soil samples by the field experiment are presented in Figure 6. As recorded in the laboratory experiment, the Brine was the best application, preventing any dust emission also in the field conditions with PM10 concentration values at the level of the background (~20 μg m<sup>−</sup>3) measured in the wind tunnel before the test.

Nonetheless, the averaged PM10 concentrations at the lowest wind velocity (4 m s−1) in all the product applications were significantly lower (*p* < 0.05) than that of the control plot and close to the background value. The averaged PM10 concentrations of the synthetic polymer (PVA) and the Bitumen were significantly lower (*p* < 0.05) than those of the control plot also in the other wind velocities (6.5, 8, and 9.5 m s−1). A significant difference between the Bitumen and the PVA was received only at wind velocity of 8 m s<sup>−</sup>1, in which the averaged PM10 concentration of the Bitumen was the lower one. In the organic products (Lignin and Resin), the averaged PM10 concentrations at wind velocities above 4ms−<sup>1</sup> were close to those of the control plot with no significant differences, excluding the Resin at 6.5 m s−<sup>1</sup> with significant lower value.

**Figure 6.** PM10 concentrations recorded during the first field wind tunnel experiment on the diverse dust control products under wind velocities 4.0, 6.5, 8.0, and 9.5 m s<sup>−</sup>1. The experiment was performed several days after the dust control application on the unpaved road with 150 rides of the quarry Euclid Trucks. Note, the "control" is an unpaved road section with no treatment. The background level measured in the wind tunnel before the test was ~20 μg m<sup>−</sup>3.

The averaged PM10 concentrations recorded in the field experiment for each product (Figure 6) were relatively lower than those of the laboratory experiment (Figure 5) because of the conditions of the soil surface structure before the experiment as explained in Section 3.1. This is less relevant for the soil plots with the application of organic products that present higher values in the field experiment at wind velocities of 8 and 9.5 m s<sup>−</sup>1. However, the PM10 emissions in both organic applications were not higher than that of the control plot. The calculated dust emission fluxes, presented in Table 3, were derived from the PM10 concentrations (Figure 6) and as such they are in full correlation with the averaged PM10 concentrations and have the same differences (and significance level, *p* < 0.05) discussed earlier. The dust suppression efficiency (%) is presented in the lower part of Table 3. It clearly shown that the Brine is the most efficient product in dust suppression, reducing more than 90% of the PM10 dust emission at the higher wind velocities 6.5, 8, and 9.5 m s<sup>−</sup>1. At these wind velocities, it revealed that the applications of the organic products were significantly less efficient in dust suppression under field conditions.


**Table 3.** The calculated PM10 emission fluxes from the unpaved roads in the first field wind tunnel experiment are presented in the upper part of the table. In the lower part of the table, the dust suppression efficiency (%) calculated as a ratio in dust flux between the treated road section and the control (no-treatment) road.

In the second stage, the experiment was conducted two weeks after the application, in which additional 296 rides (total of 446) of Euclid Trucks were recorded. The results of the dust suppression efficiency (%) are presented in Table 4. The organic products (Lignin and Resin) were found to be non-efficient applications for dust suppression by wind. The efficiency of the Bitumen and the PVA reduced dramatically compared to the first stage (Table 3) with only minor dust suppression at the higher wind velocity (14.3% and 3.9%, respectively). It revealed that all these product applications were destructed by the truck wheels after 150 rides, leading to dust emission at the same rate of the control (no treatment) road section. However, the application of magnesium-chloride solution (Brine) on the road surface kept on the same high-level of efficiency as in the first stage.

**Table 4.** The efficiency (%) in dust-PM10 emission recorded two weeks after the dust control application on the unpaved road with 446 rides of the quarry Euclid Trucks. The dust suppression efficiency (%) was calculated as a ratio in dust flux between the treated road section and the control (no-treatment) road.


#### *3.3. Improved Application of Brine*

The high efficiency of the magnesium-chloride solution (Brine) in dust emission suppression because of vehicle transportation was already reported [18,24]. The results of this study indicate clearly on the high efficiency of the magnesium-chloride solution in dust suppression because of wind erosion in calcareous soil bed in arid climate. The magnesium-chloride solution (Table 1) provides soluble salts to the soil that act as a cementing material for the formation of soil aggregates and soil structure [25]. Moreover, its hygroscopic property allows absorbing water from the atmosphere, even in low atmospheric humidity, agglomerating the fines and holding the aggregate matrix together through suction forces [18]. The existence of water enables the cohesive forces between the soil particles, in particular the clay-sized particles [26]. In turn, the higher soil water content, the lower wind erosion and dust emission because of the inter-particles connections and larger aggregates that are more resistance to the wind shear stress [14]. Moreover, the Brine hydrous magnesium chloride can be a byproduct, a waste material, from industrial production. In Israel, this Brine is available by the production process of potassium (K) from the saline Dead Sea. The use of waste materials is important in sustainable environment perspective.

The recommended application of the Brine is wetting the soil surface at the rate of 1.5 L m−2, based on the very high efficiency of the Brine to suppress dust emission in the calcareous soil. It was assumed in this study that the application rate can be reduced and still keep on high efficiency. In order to test this, five samples with different rate of Brine concentration (amount of solution in the soil) were prepared for laboratory experiment: 1.5 L m<sup>−</sup>2, 1.2 L m−<sup>2</sup> (that is 80% of the recommended rate), 0.8 L m−<sup>2</sup> (60%), 0.5 L m−<sup>2</sup> (40%), and 0.3 L m−<sup>2</sup> (20%). After the application process, the soil samples were left outside for drying for several days. The wind tunnel experiment was performed in the same procedure described in Section 2.1 with two wind velocities: 6.5 and 9.5 m s<sup>−</sup>1. The dust-PM emission fluxes were derived from the PM10 concentration measured during the experiment. Then, emission coefficient was calculated as the ratio between the dust fluxes of the Brine sample (different concentration) and the control (no-treatment) sample, in which, the lower value, the higher dust emission suppression. The dependence of PM10 emission in the Brine concentration applied on the soil under two wind velocities, 6.5 and 9.5 m s−1, is presented in Figure 7. The results indicate that the Brine concentration can be reduced to 60% (use of 0.8 L m−<sup>2</sup> instead 1.5 L m−2) and 40% (use of 0.5 L m−<sup>2</sup> instead 1.5 L m<sup>−</sup>2) for winds of 9.5 m s−<sup>1</sup> and 6.5 m s<sup>−</sup>1, respectively, and still keep on high efficiency of >90% in dust emission suppression. Although the high efficiency of the Brine to prevent dust emission by wind in both laboratory (Table 2) and field (Tables 3 and 4) conditions, the exact level of reduced Brine concentration that is recommended to be applied in the field should be determined by a targeted experiment.

**Figure 7.** The dependence of PM10 emission in the Brine concentration applied on the soil under two wind 6.5 and 9.5 m s<sup>−</sup>1. The efficiency is represented by the emission coefficient (*Y* axis) that is the ratio between the dust fluxes of the Brine sample (different concentration) and the control (no-treatment) sample, in which, the lower value, the higher dust emission suppression.

#### **4. Conclusions**

Previous works on dust control in unpaved roads have mainly focused on dust emission by vehicle (wheels) traveling. The results of the control soils (no-treatment) in this study clearly show that dust emission by wind from unpaved road can substantially contribute to the mass transfer, and thus should be considered. The experiment on dust control products in laboratory conditions with no vehicle transport on the soil revealed high efficiency in dust emission suppression for all the products (Brine, Lignin, Resin, Bitumen, PVA), reducing at least 77% of the dust emission by wind. However, the efficiency of most of the products was reduced significantly when applied and tested under field

conditions in quarry-haul truck road. Two weeks later and after 446 rides of the quarry trucks, the products were not efficient to prevent dust emission by wind, except of the Brine.

The Brine enables to keep on very high efficiency (>90%) to suppress dust emission by wind over time. Additional experiments revealed that Brine can be applied with reduced amounts and still keeping on low emission. There are only few scientific articles on the application of hydrous magnesium chloride (Brine). Those studies showed the high efficiency of the Brine to reduce dust emission due to vehicle (wheels) traveling. In this study, the Brine was found to be very efficient also in wind erosion and may provide a good solution for product application in disturbed soils especially in the water-limited areas of the arid and semi-arid climates. Nevertheless, further study is needed to investigate the possible environmental impacts of the diverse dust suppression substances, including toxicity of atmospheric particulate matter when dust is emitted from the treated soils and/or soil-groundwater pollution because of vertical fluxes of the applied solutions on the surfaces.

**Funding:** This research was funded by the Ministry of Energy, Israel, grant number 212-17-021.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Technical Note* **A Clay-Based Geopolymer in Loess Soil Stabilization**

#### **Nadav Hanegbi and Itzhak Katra \***

Department of Geography and Environmental Development, Ben Gurion University, 8410501 Beersheba, Israel; nadavhan@post.bgu.ac.il

**\*** Correspondence: katra@bgu.ac.il

Received: 25 March 2020; Accepted: 6 April 2020; Published: 10 April 2020

**Abstract:** Soil erosion has environmental and socioeconomic significances. Loess soils cover about 10% of the global land area. Most of these soils are subjected to increased land uses such as unpaved roads, which increase soil destruction and dust emission to the atmosphere. There is a significant interest in applications for dust control and soil stabilization. Application of geopolymers may significantly reduce environmental impacts. This study examines the use of a metakaolin-based geopolymer for dust control and soil stabilization in a semi-arid loess soil. The application of the geopolymer for dust control in comparison with common products (brine, bitumen, polyvinyl acetate-PVA) resulted in no dust emission. As a soil stabilizer, the geopolymer tested in this study provides remarkably good results in the tensile test. The most successful composition of the geopolymer, which is activation solution of sodium silicate and sodium hydroxide (NaOH) together with an addition of 30% metakaolin, obtained soil strength of 23,900 N after 28 days. The attempt to replace NaOH with lime (CaO) in the activation solution was far inferior to the original composition. There is a strong potential to develop natural soil stabilizers from a mineral base that even surpass their capabilities over existing synthetic stabilizers.

**Keywords:** loess; metakaolin; dust control; geopolymer; soil erosion

#### **1. Introduction**

Soil erosion has environmental and socioeconomic significances. In arid areas, soil erosion by wind leads to dust emission to the atmosphere that increase air pollution and risk exposure to human health [1–3]. Current global estimates of dust-related particulate matter (PM10), which is less than 10 μm in diameter, loading varies widely from ~6 to 30 million tons [4]. There is uncertainty regarding the relative dust contribution of disturbed soils in non-desert areas, which are subjected to human activities. The continuous population growth increases the land uses of agriculture, mining, and vehicles traveling on unpaved roads in these areas. Such activities change the topsoil properties, reduce the soil stability, and thus increase soil erosion and dust emission [5–7].

There is a growing interest in soil stabilization to reduce soil erosion and dust emission from disturbed soil in arid and semi-arid regions. There are various dust-control materials commonly used to reduce dust emission from soils. These products are based on natural and synthetic polymers such as lignin, resin, bitumen, and polyvinyl acetate (PVA). The typical application of these products is by spraying the solution on the soil surface at certain concentration (liter per square meter). Among the common dust-control materials, the hydrous-magnesium chloride (brine) is one of the most efficient application in preventing dust emission [8]. The adverse soil pollution and environmental footprint resulting from the use of the common dust-control materials are still unclear. Moreover, typical applications for dust control do not provide a reinforcement of the topsoil, and thereby stabilization against physical pressures due to vehicle traveling.

One of the most common materials that is used for soil stabilization is PVA. When the PVA solution is absorbed in the soil, the polymers enveloping the soil particles that remain aggregated after the evaporation processes. PVA is especially efficient in soils with granular texture and with small specific surface area, in which it creates a stable crust with strong tensile strength. However, the environmental impacts of this synthetic polymer and its efficiency in fine-grained soils with high content of dust are still unclear. This raises the need to develop an effective material from the perspectives of dust control, soil stabilization, and environmental footprint. Geopolymers are a group of mineral-based materials that commonly used in the building industry [9]. Geopolymer can be synthesized by mixing calcined kaolin and strongly alkaline solutions (such as NaOH or KOH) and then curing at a room temperature [10]. Several studies have examined the use of various geopolymers for loess soil stabilization, including fly ash, lime, and cement, etc., [11–16]. However, the use of the anhydrous calcined form of the clay mineral kaolinite (metakaolin) was given less attention in loess stabilization. Loess is a widespread surface deposit in many parts of the world. There is a strong need for loess stabilization in construction requirements and geohazard mitigation due to its propensity to collapse and subsidence after loading and wetting [11,12]. Moreover, unpaved roads in loess are associated with strong erosion by wind and water, which leads to soil disintegration [7]. This study examined the use of a metakaolin-based geopolymer with different compounds and concentrations for the reinforcement of the topsoil of a semi-arid loess soil.

#### **2. Materials and Methods**

The experiments conducted on a typical loess soil of the northern Negev, Israel. The semi-arid region of the northern Negev is characterized by annual average rainfall of ~200 mm. Most of the loess soils in northern Negev are a secondary material that reworked under pedogenic and geomorphological processes after investment of aeolian and/or pluvial materials [17]. Currently, the northern Negev loess is subjected to intense human activities and related soil erosion and dust emission, e.g., [6,18].

The texture of the soil is silt-loam with a relatively high content of silt and clay particles. Particle size distribution was obtained by the ANALYSETTE 22 MicroTec Plus (Fritsch, Idar-Oberstein, Germany) laser diffraction, which measures particles in the size range of 0.08–2000 μm [19,20]. Replicas (100 mg) were dispersed in Na-hexamvetaphosphate solution (0.5%) by sonication (38 kHz). The data was calculated using the Fraunhofer diffraction model with a size resolution of 1 μm using MasControl software (Figure 1a). Mineralogical analysis, which was obtained by the X-ray power diffraction (XRPD) method, shows typical loess composition with high content of quartz, clays (illite, kaolinite), as well as Ca-based components that are found in arid loess (Figure 1b).

For the dust control experiment, the soil was placed in trays that fit the wind tunnel dimensions (0.5 × 1.0 m, with height of 0.02 m). Three trays (replicas) were prepared for each soil treatment/product. Table 1 presents the products used for the laboratory experiment and their application on the soil. The composition of the solution of each substance and the concentration of the solution to be applied by spraying on the soil surface were provided by the suppliers of the products and are typical in dust control. The soil samples were left outside for a drying process of several days before the experiment to avoid the impact of the amount of the solution (water content) on wind erosion and dust emission. The geopolymer for the dust control experiment was composed of metakaolin (MK) with an activation solution containing NaOH and sodium silicate in a constant ratio. The activation solution was added to the soil mixture in over-saturated until reaching a shed state of the mixture. The MK amount added to the soil sample was 30% of the soil sample weight. The geopolymer was mixed in the soil sample and tested after a drying process of 28 days.

The soils were tested in a wind tunnel for dust emission—overall 30 runs were performed by the wind tunnel (5 samples, 3 replicates, 2 wind velocities). These experiments were performed in the portable boundary-layer wind tunnel of the Aeolian Simulation Laboratory at BGU [8]. Boundary-layer wind tunnels enable simulations under standardized quasi-natural wind conditions and provide quantitative information on dust emission rates from the soils. The wind tunnel has a cross sectional area in the order of 0.5 × 0.5 m with open-floored working sections of up to 10 m length. Air push or air suction flow in the tunnel is generated by an axial fan up to a maximum velocity of 18 m s−1. Instruments installed in the test section of the tunnel enable quantification of wind characteristics and dust (PM) concentrations.

**Figure 1.** Characteristics of the Loess sample used for the experiments: (**A**) particle size distribution obtained by the laser diffractometer technique with the statistical parameters (right side). (**B**) Mineralogical composition of the loess by the X-ray power diffraction (XRPD).

**Table 1.** List of products for dust control in the loess soil used in the laboratory experiment. The composition of the solution of each substance and the concentrations of the solution to be applied by spraying on the soil surface were provided by the suppliers and are typical in dust control. The geopolymer was prepared and applied as described in the text.


The tunnel was operated under two wind velocities, 6.5, and 9.5 m s<sup>−</sup>1, which represent medium and strong natural winds that are associated with dust emission in semi-arid loess soils [13]. Each run of specific wind velocity lasted 30 s to record the trend of dust emission. Dust concentrations (μg m<sup>−</sup>3) PM10 were recorded by a light-scattering device, DustTrak DRX 8534 (TSI, Shoreview, MN, USA), in the range of 0.001–150 <sup>μ</sup>g m−<sup>3</sup> (±0.1% of reading) at 1-s intervals, which was placed at 25 cm above the tunnel bed. The recorded PM10 concentrations were converted into mass flux emitted from the soil surface (mg m−<sup>2</sup> s<sup>−</sup>1) based on the wind tunnel dimensions, the area of the soil bed, and the wind

velocity (for in depth description refer to Katra et al. 2016 [21]). The PM10 flux from the soil surface can be used in parameterization of dust emission and transport models [7,21].

The soil sample mixed with the geopolymer was examined also for soil stabilization by the tensile strength test. A 300 g quantity of soil was mixed with MK (with the activation solution) to form test sample in three replicates. The mixtures were then poured into a 3.5 cm diameter cylinder on a length of 7.5 cm. Each sample series was dried for 28 days until the bonding material stabilized before testing the strength of the stabilized soil sample. The tensile strength of the sample was measured using the Universal Testing Systems for Tensile, Compression, and Flexure Testing of Instron 5900 Series (Instron, Norwood, MA, USA) at a crosshead speed of 10 mm min<sup>−</sup>1. This device is commonly used to measure the strength of samples in cement engineering and other materials. The soil sample was placed at the center of the lower surface of the pressure clamps for pressing small pulses until the graph stabilized after the collapse of the sample. The strength tests were conducted for two groups of geopolymer compositions (Table 2). (1) Geopolymer with different MK percentage (0–30%) in the soil sample. (2) Metakaolin (MK) based geopolymer with an alternative activation solution, replacing the NaOH with lime (CaO). Different concentrations of CaO were tested, where the measured solubility level of lime at 20 ◦C (1.79 g L) was taken as a starting point. The amount of lime suspension (in water) used in the activation solution was the same as the NaOH in the origin activation solution, and up to 20-times. The content of the MK was 30% in this group.

**Table 2.** The compositions of the geopolymer as a loess soil stabilizer used for the tensile strength test. Group 1: different MK content (%) with a constant activation solution of sodium silicate and NaOH. Group 2: different concentrations of CaO as a replacement to the NaOH in the activation solution with MK content of 30%.


#### **3. Results and Discussion**

Figure 2 presents the results of the dust emission experiment. The recorded PM10 concentrations are a result of dust emission from the surface of the loess soil samples under two wind velocities. It is clearly shown that no dust emission and PM10 concertation was recorded in the application of the clay-based geolopyler. The "straight" line that was similarly recorded for the two wind velocities in the test of the geopolymer application is a representation of the background level (~20 μg m−3) measured in the wind tunnel before the test. In all the other applications, bitumen, brine and PVA, the averaged PM10 concentrations were significantly (*p* < 0.05) higher than the background level. In general, the PM10 concentrations received for the brine, PVA, and bitumen, were higher under 9ms−<sup>1</sup> compared with 6.5 m s<sup>−</sup>1. The wind velocity OF 9.5 m s−<sup>1</sup> was run promptly after the prior velocity (6.5 m s−1), and since the amount of dust particles in the soil samples is limited, the dust emission in 9.5 m s−<sup>1</sup> could have been an underestimation for non-limited samples. The PM10 values of the control sample were significantly higher than those of all the other samples in both wind velocities, suggesting that all the applied products are efficient for dust suppression at a certain level. It is interesting to see that the brine, which was found to be very efficient by reducing > 90% of the dust emission in calcareous soils under the same experimental conditions [8], is less effective in loess soils. The results of the brine are at the range of the PVA. However, the brine is not used for the purpose of soil stabilization as the PVA and the geopoylmers.

**Figure 2.** PM10 concentrations recorded during the laboratory wind tunnel experiment on the diverse dust control products under wind velocities 6.5 (**right** side) and 9.5 m s−<sup>1</sup> (**left** side). Note, the "control" is the Loess soil with no treatment. The background level measured in the wind tunnel before the test was ~0.020 mg m<sup>−</sup>3.

Tensile tests were conducted for the geopolymer as a soil stabilizer. The first test was to examine the effect of the MK amount at the range of 0–30% in the geopolymer composition (Group 1, Table 2). The activation solution (sodium silicate and NaOH) remained in a constant ratio. The results in Figure 3a show that the soil sample without MK (0%) resulted in resistant of less than 1000 N. Application of MK content at the range of 2.5–10% improved the soil resistant up to 6000 N. A significant change was obtained for samples with MK content at the range 20–30%. These samples resulted in a consistent increase in resistance to pressure when the peaks ranges from ~11,000 N (20% MK) to ~25,000 N (30% MK). Just for comparison, we tested also the PVA as soil stabilizer, which has been used extensively in soil stabilization [22]. The PVA results after a curing time of 28 days were up to 4200 N. The drying time until significant strength is reached is an important factor when there is a need to stabilize an unpaved road. External factors, such as solar radiation and water regime, can have a great impact on the nature and duration of the soil stabilization. As such, this study examined the hardness of the geopolymer during a short seven-day cure. The geopolymer composition of 30% MK with the activation solution (sodium silicate and NaOH) yielded a result of 15,510 N.

**Figure 3.** (**A**) Resistance to load of stabilized Loess soil with different contents of Metakaolin (MK, 0–30%) in an activation solution containing sodium silicate and NaOH (Group 1, Table 2). (**B**) Resistance to load of stabilized Loess soil with calcium oxide (CaO) as a substitute for NaOH in the activation solution (Group 2, Table 2). The results are CaO at the same amount as the NaOH (X1), which was used in the geopolymer presented Figure 3A, and with enlarged amount of CaO up to 20-times (X5-X20). The content of the Metakaolin was 30% for all samples in Figure 3B.

A second tensile test was conducted for geopolymer with replacement of the NaOH with the base of lime (CaO) that is more "natural" soil material (Group 2, Table 2). The addition of lime is used in applications of soil loess stabilization, e.g., [23]. The lime improves the cation exchange of the soil and absorbs water, which lowers the pH. Adding lime also affects the structure of the soil through the soil ability to absorb water. The results of the test for soil samples stabilized by the geopolymer with different amounts of CaO in the activation solution are presented in Figure 3b. The MK content of 30% was constant for all the samples. The resistance increased as the concentration of the lime was higher. This was true up to lime concentration that is 15 (X15) times that of the original concentration (X1). At the highest concentration (X20), the resistance was lower (peak at 6450 N) than in X15 (peak at 7400 N). This can be explained by a rapid and unsteady solidification of the material, which does not allow us to pour the soil mixture into the molds. However, the results of the NaOH in the activation solution (Figure 3a) are significantly higher than those of the replacement base (CaO) with 30% MK (Figure 3a).

The results in Figures 2 and 3 show the potential of the clay-based geopolymer examined in this study to suppress dust emission and to stabilized loess soil. The application of the geopolymer for dust control in loess soil in comparison with common products (brine, bitumen, PVA) resulted in no dust emission. As a soil stabilizer, the geopolymer tested in this study provides remarkably good results in the tensile test. The most successful composition of the geopolymer, which is activation solution of sodium silicate and NaOH together with an addition of 30% MK, obtained soil strength of 23,900 N after 28 days. The attempt to replace NaOH with CaO in the activation solution was far inferior to the original composition. The factor of relatively rapid firming of the soil after seven days (strength of 15,510 N) is very important in stabilizing unpaved roads.

#### **4. Conclusions**

The clay-based geopolymer tested in this study showed high potential for application as a dust control product as well as soil stabilized in loess soils. There are important factors that should be considered in the development of such a geopolymer. (1) First, some elements of the activation solution are substances that require production by chemical or physical processes such as MK and sodium silicate, which involve certain costs not taken into account in this study. Second, this study did not take into account the effects of the final product on the environment over long time after the application in the soil. Nevertheless, there is a strong potential to develop natural soil stabilizers from a mineral base that even surpass their capabilities over existing synthetic stabilizers. The ability of these materials to generate reactions with the soil components will allow us to reach high levels of strength of the upper crust of the stabilized soil.

**Author Contributions:** Conceptualization, I.K.; methodology, I.K. and N.H.; formal analysis, N.H.; writing —original draft preparation, N.H.; writing—review and editing, I.K. and N.H.; funding acquisition, I.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by KKL-JNF, grant number 03-02-040-14.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18