Introduction:

a) Soil erosion “is a complex and dynamic process representing continuous removal of [the] top soil layer…” (Choudhury and Sil 2010) that occurs over many different scales.  The research problems are that “accelerated soil erosion is among the most pressing…environmental problems, resulting in degradation of ecosystem function and decreased productivity and sustainability of agriculture…” (Bowker et al. 2008).  The United States has lost “3.06 billion tons of soil in 1982 and 1.73 billion tons of soil in 2007” from erosion on cropland (NRCS 2007) in areas with industrial agriculture. “Ninety percent of U.S. cropland is losing soil…17 times faster than formation…” (Kimbrell 42).  Accurately calculating soil erosion and the resultant impacts “is one of the greatest challenges in natural resources and environmental planning” (Bhuyan et al. 2002) due to dynamic changes in physical and chemical weathering and erosion rates on landscapes from water, wind and ice.

Simulation methods use empirical measurement and models with calibrated parameters using methods from geomorphology or fluvial geomorphology.  The equations are generally related to the elevation, local relief, slope angle, infiltration, soil traits/erodibility, vegetation presence, weathering and climate, transport (e.g., rill flow, overland flow, etc.) and deposition of materials over space and time.  “Measurement techniques differ in accuracy, equipment and personnel cost” (Stroosnijder 2005).  Models, like “the Water Erosion Prediction Project (WEPP)” (Bhuyan et al. 2002), may provide the easiest and most cost-effective way to estimate erosion (based on the change in elevation) over time.  Often, “to address soil and water resource problems…for large basins…hydrologic models are coupled with sediment transport models to increase prediction accuracy” (Choudhury and Sil 2010).  For larger spatial extents, erosion is often estimated or calculated by “sediment collection…[or]…dam [measurement] comparison” (ibid.) or some version of Remote Sensing (RS).  The latter method might involve using bi-temporal change detection to find the delta of elevation.  “The yield of sediment…is a complex process responding to all variations…in precipitation, soils, vegetation, runoff and land use” (Langbien and Schumm 1958).

b) The reason that complex spatial modeling approach is necessary to complete the objective is to create a complex process-response model that emulates “time variation in sediment discharges” (Choudhury and Sil 2010), at a large spatial extent, to “mitigate [potential] risks” (Tramblay et al. 2009) and examine how local interactions lead to system level patterns.  Models are necessary to better “manage…natural and disturbed [fluvial and terrestrial] systems” (Coulthard et al. 2012) “…to understand…ongoing erosion processes” (Stroosnijder 2005) so agencies can prioritize and “allocate conservation resources and develop conservation regulations…[for] erosion control technology…” (ibid.).  Calculating soil erosion is problematic because of “…the large temporal and spatial variation…[and] the paucity of accurate erosion measurements” (Stroosnijder 2005).  Soil erosion is due to a number of spatially random variables (climate and runoff), uncertain variables (flow direction) and feedback patterns (from flow accumulation, soil and geological characteristics) that result in un-deterministic outcomes, making prediction difficult if not impossible.  “Recent investigations have shown that the extent of the channel network in some drainage basins is controlled by a threshold for overland flow erosion” (Tucker and Slingerland 1997).

c) The objective of this project is to look at how changes in environmental (rainfall based-runoff and soil hardness – erodibility) and biological (vegetation) parameters affect soil erosion on terrain and hill slope environments.  The project goals are to identify what and how local interactions between rainfall-based runoff, soil hardness and vegetation lead to higher-level spatial patterns of cumulative erosion (overall lowering of the landscape).  The model calculates soil erosion (A) based on “…properties of soil which influence soil erodibility…” (Bowker et al 2008) using soil hardness (K) and “plant canopy cover” (ibid.) (C).  In this model, erosion is A = ((1 – A) – C).

d) The project will be completed using a complex systems approach to examine the agent interactions.  The erosion model (Dunham et al 2004) from the NetLogo (Wilensky 1999) was used as the base model, modified for this project by adding a vegetation parameter.  The model emulates rainfall, erosion and drainage network to look at how environmental variables affect soil erosion.  A sensitivity analysis was completed for both terrain and hill slope based watersheds with and without vegetation.

Methods Section:
Overview
Purpose:
The purpose of this model is to examine how rainfall-based runoff and vegetation affect erosion in a watershed through space and time.

Entities, state variables and scales:
The entities are rainfall representing the climate, the land representing environmental space (the watershed) and drains that represent low-elevation exit points for the system runoff (e.g., tributaries).

The state-variables are elevation, water and erosion.

The fixed spatial scales represent the environment (151 x 151 units), appropriate for a watershed (kilometers) with a cell size of 4 x 4 units.  The temporal scale is synchronous, each time-step represents a daily amount of time.

Process overview, scheduling:

The agents are placed as land (bumpy terrain or hill) then drains with some elevation.  Rainfall begins at t = 1 accumulating with every time-step, “initially…[creating]…lakes” (Dunham et al.).  Drains are reset to a negative elevation at every time-step where runoff will exit the system.  Water flows from a neighbor to a contiguous neighbor with low elevation and erosion occurs at every time step.  At each time-step, the state variables (elevation, erosion and water) are updated and the process repeats.  Time is modeled as continuous.

Design concepts

Basic principles:

The model appears to be based on “the Revised Universal Soil Loss Equation (RUSLE)” (Bowker et al. 2008) using some of the equation parameters.  Flow accumulation and direction are determined by elevation and rainfall.  Water accumulates in areas of low elevation and flows to the lowest land unit (low cell of Moore neighborhood) simulating gully erosion and drains out at the edges of the environment.  The abstract erosion model takes soil hardness, rainfall and vegetation into account when calculating soil erosion.

Emergence:

The emergence spatial patterns produced are landscape erosion and channelization (drainage network) on the landscape over time.  The erosion and channelization patterns vary based on the parameter settings.

Objectives:
The objective is to emulate rainfall-based erosion and drainage network formation over space and time.  Water flows downhill carrying sediment with it, creating fluvial channels.  Erosion is measured cumulatively as well as the cells with water at the system and agent level.

Sensing:

Sensing occurs (but not displayed) when agents (Moore neighborhood) compare elevation and flow to neighbors to find the flow direction of runoff.  Sensing is an explicit process within the flow direction within the model.

Interaction:

Direct interactions are between the land and water.  Path dependence occurs as river channels where runoff will likely continue in the same direction of areas of lower relief (feedback).  Communications are represented through state comparison.

Stochasticity:
Stochasticity is implemented in this model by rainfall on random areas over the landscape through time.  Elevation is semi-stochastic with a terrain setting, but is uncertain because upper and lower elevation limits are specified within the model.

Observation:

Observation includes the measurement of cumulative erosion and the number of landscape cells where erosion occurs over the landscape.

Details:

Initialization:
The initial state (at t = 0) includes the landscape (hill, bumpy or flat terrain) with a roughness (6) and drains at the edge of the environment.  The environmental parameter settings are rainfall (set at 0.10), vegetation (set at 0.20) and soil hardness (set at 0.30).  The states of the entities are elevation at the maximum height because no erosion or water has occurred.  The values of the state variables are set arbitrarily.  Initialization is the same and includes the landscape parameters  (terrain roughness, hill slope or karst), environmental parameters (rainfall, vegetation and soil hardness), but the landscape (world) is slightly different every time.

b) The sensitivity analysis will be conducted (behavior space tool) for the base erosion model and the modified model (with a vegetation parameter) over ten iterations with incremental changes in the parameter settings for both rough terrain and a hill slope.  The rainfall, soil hardness (K) parameters will be evaluated for values 0.25, 0.5 and 1 for both models and landscape.  The vegetation parameter will be evaluated for values 0, 0.5 and 1.

c) The model will be evaluated by comparing the results of the sensitivity analysis (hill slope and rough terrain for the base model and abstract model with vegetation) against results from scientific literature.

3. Results Section:

a) The basic model represents a landscape on a regular 151 x 151 square lattice.  Rain falls randomly over the landscape flowing downhill. Soil erosion occurs as water flows over the landscape.  Interactions are between the land and the water resulting in soil erosion from neighbor to a contiguous neighbor with the minimum elevation within any neighborhood (Moore neighborhood).  Erosion changes through time following:

A_t+1 = ((1 – K) – C)

Where, t is the current time-step, K is the erodibility resistance (soil hardness) and C is the vegetation.  The environmental inputs apply equally to all land units.  It is the heterogeneity of the landscape (elevation differences), land-rainfall interactions, feedback and path dependence that produce the emergent spatial patterns.  Cumulative erosion changes as follows:

∑ A (t_{1 – ∞)}

At the cellular level, the environment starts as a rough landscape (t = 0).  At each time-step, rain falls, water accumulates and flow towards a neighbor with less elevation.  The flow is half of the delta of water plus elevation, if the amount is greater than zero and erosion occurs (sediment is carried downhill to that neighbor).  The elevation of each land unit is recalculated by subtracting the amount of soil erosion from the starting elevation (z₁) of that cell.  As erosion continues drainage channels develop creating a pathway from erosion (path-dependence) within neighborhoods.  Drains are set to an extreme negative elevation and reset at each time-step.  The drain cells provide an outlet for runoff to exit the system.  This entire process repeats at every time step.

At the system level lakes develop over the landscape because of rainfall.  Erosion occurs, lakes are drained and a drainage network develops. Flow direction is variable at first but as the system becomes more stabilized (channelization), runoff flows in similar directions.  Runoff becomes path dependent from the erosion-based feedback.  A positive feedback results from erosion leading to path-dependence creating a drainage network over time.  Path-dependence also occurs because of the stochasticity of rainfall.  The emergent spatial patterns are denudation and enlarged drainage channels because the system is in disequilibrium.  Thresholds could not be attained because of the variability among iterations.

The image sets and videos below show the parameter settings (left) where the plot shows cells with water (in salmon) and without (in brown) and the monitor shows cumulative erosion (parameters are applied equally).  The emergent patterns that are produced after 730 time-steps (middle image) and the individual land agent (right) where the plot shows elevation (in brown) and water (in blue).  Each image set has a video that is matched to the image set.  The sets are variations of the rainfall parameter for rough terrain (first two sets) and hill slope (second two sets) using the abstract model.

 

 

b) The findings from the sensitivity analysis indicate the hill slope model results in less erosion than the bumpy terrain in the abstract model (not in the baseline model).  The abstract and base model had about the same maximum cumulative erosion.  In the base model both the first quartile and the minimum were equal to zero while in the abstract model both of those were negative and the median was zero.  Erosion was most sensitive to rainfall but was highly dependent on the other parameter settings (vegetation and soil hardness).  Erosion is higher when the soil hardness and vegetation are low and rainfall is high.  The median value shows a difference of 45.6 thousand units of erosion between the abstract and base model.  The third quartile value is double in the base model at (81.7 thousand units) indicating that vegetation minimizes erosion and may build up soil over time. This model emulates a supply limited environment where sediment production is less than erosion rates.

As the rainfall parameter increases the cumulative erosion doubles.  There is a proportional relationship between rainfall, soil hardness and vegetation (the formula that calculates erosion uses vegetation and soil hardness parameters).  If the soil hardness and vegetation parameters equal or add to one, there will be negative erosion, and soil aggradation occurs.  Erosion is less sensitive to terrain smoothness because the parameter creates variability in elevation (not specifically part of the erosion calculation).  When rainfall is high the drainage network develops faster.  The base and abstract sensitivity analysis results are displayed below for comparison.

c) The model shows continuous erosion. However, “sediment yield is at a maximum at about 10 to 14 inches of precipitation, decreasing sharply on both sides…[because] of deficiency of runoff…[or] increased density of vegetation” (Langbein and Schumm 1958).  The results of the model match real-world observations.  In a semi-arid environment, 12 inches of annual precipitation result in an annual sediment yield of ~825 tons per square mile (ibid.).  In a forested watershed, 60 inches of annual precipitation, the annual sediment yield is ~250 tons per square mile (ibid.).

The results indicate that cumulative erosion of 45,000 cubic feet (the median value of the base model or third quartile of the abstract) would be ~2,127 tons per year.  The watershed sediment yield (151 square kilometers) would result be ~36.6 tons of erosion per square mile.  The maximum would be ~175,000 cubic feet with a sediment yield of ~120 tons per square mile matching what Langbien and Schumm measured in a desert watershed. However, this environment would have less than 5 inches of precipitation annually.

Sediment yields from reservoirs are higher, averaging 1,100 – 1,500 tons per square mile in places like California, Wyoming, Utah and Colorado (Langbien and Schumm). “Several studies have shown that sediment yields have decreased with increased drainage area” (ibid.).  Sediment yields are highly dependent on the lithography, vegetation and precipitation.

Comparing the results from this model to other models from research by B. Cheviron et al., this model does not agree well with other erosion models.  There are differences in spatial scale (of the watershed) and methods.  The models, in order from small-scale to large-scale spatial extent applications, are:  Distributed Hydrological Modeling for Agrosystems (MHYDAS), the Semiempirical Sealing and Transfer by Runoff and Erosion Related to Agricultural Management (STREAM), the Pan-European Soil Erosion Risk Assessment (PESERA) and the Regional Modeling of Soil Erosion Risk (MESALES) (Cheviron et al.).  MHYDAS seems to be the closest to the abstract model in terms of the extent and continuous rainfall.  Erosion for the abstract project model, discussed above (for 2,127 tons/year), would have an annual soil loss of 0.14 tons per hectare.  The models from the article resulted in similar values for runs using lower precipitation values.  Their model setup and calculations are not completely clear.

When comparing the abstract model with the models from the research by Bhuyan et al., there is more agreement.  The models included WEPP, the Erosion Productivity Impact Calculator (EPIC) and the Areal Nonpoint Source Watershed Environment Response Simulation (ANSWERS) (Bhuyan et al.).  These models evaluated different cultivation methods (no till to plow) with the WEPP and EPIC models resulting in similar values.

4. Discussion:

This Cellular Automata (CA) model emulates sheetwash, overland flow, rill and gully erosion for either a hill slope or rough terrain landscape system (middle numbered system occurring at the human time-scale).  The model shows how interactions between water and land result in soil erosion, sediment transport and drainage network formation over space and time.  Rainfall is the stochastic component in this model.  Autonomous decisions occur when runoff flows in an unbiased direction based on the elevation delta (Moore neighborhood).

The results of rainfall and erosion on the land produces positive feedback encouraging more erosion at the neighborhood level.  These feedbacks are a result of erosion and channelization that reinforce channel development.  The feedbacks and the stochastic rainfall create path dependency in flow direction making it more likely that the flow direction will continue to the same lower point within a neighborhood as channels become enlarged over time.  At the system level, feedback is negative because the system is moving towards equilibrium. Rainfall, vegetation and soil hardness are all linear components that lead to non-linear erosion results, producing higher-level drainage patterns (drainage and erosion).

Preventing erosion can be achieved by not disturbing sensitive areas or by re-establishing vegetation on hill slopes and terrain.  Preventing gully erosion is particularly important in minimizing accelerated erosion. This model could be enhanced by adding parameters for slope angle, slope length, infiltration, uplift and different climate types (e.g., arid, semi-arid, forested, etc.).  Currently, the model simulates continuous rainfall but there are very few areas where rainfall is so extensive (i.e., it doesn’t rain all year).  Another issue is the temporal scales, a drainage network could take millennia to develop.  In any case, reporting units need to be adjusted accordingly to fit the watershed type, temporal and spatial scales.

5. Conclusion: This project examined how local interactions  lead to soil erosion due to rainfall-based runoff resulting in denudation and drainage network formation over space through time.  It is important to understand how local processes lead to higher-level patterns because soil is key to our food systems.  It is also necessary to prevent catastrophic land use changes (e.g., deforestation, agricultural operations) in unstable landscapes.  This project used an abstract of the erosion model (Dunham et al 2004) from NetLogo (Wilensky 1999) by adding a vegetation parameter.  A sensitivity analysis for cumulative soil erosion was evaluated for terrain and hill slope watershed.  Interactions occur between land and water show variability in flow direction until channels develop when path dependency is created by positive feedbacks and stochasticity.  At the system-level, feedback is negative driving the system towards stability.  The emergent spatial patterns that are produced over time are a drainage network and denudation of the landscape.  Not disturbing the soil, especially in “semi-arid” (Langbien and Schumm 1958) environments, and reestablishing vegetation (where applicable) is the best strategy for building soil tilth and quantity increasing infiltration and water retention, decreases soil exposure and can reduce or prevent soil erosion.

 

Citations

Bowker, Matthew A., Jayne Belnap, V. Bala Chaudhary and Nancy C. Johnson. “Revisiting Classic Water Erosion Models in Drylands: The Strong Impact of Biological Soil Crusts.” Soil Biology and Biochemistry 40 (2008): 2309-316. Soil Biology and Biochemistry. Elsevier. <www.elsevier.com/locate/soilbio>. Web. May 27, 2015.

Bhuyan, Samar J., Prasanta K. Kalita, Keith A. Janssen and Philip L. Barnes. “Soil Loss Predictions with Three Erosion Simulation Models.” Environmental Modeling & Software 17 (2002): 135 – 44. Environmental Modeling & Software. Elsevier. <www.elsevier.com/locate/envsoft>. Web. May 27, 2015.

Cheviron, B., Y. Le Bissonnais, J.F. Desprats, A. Couturier, S.J. Gumiere, O. Cerdan, F. Darboux and D. Raclot. “Comparitive Sensitivity Analysis of Four Distributed Erosion Models.” Water Resour. Res Water Resources Research 47 (2011). Water Resources Research. <http://onlinelibrary.wiley.com.libproxy.uoregon.edu/doi/10.1029/2010WR009158/abstract>. Web. June 5, 2015.

Choudhury, Parthasarathi and Briti Sundar Sil. “Integrated Water and Sediment Flow Simulation and Forecasting Models for River Reaches.” Journal of Hydrology 385 (2010): 313 – 22. Journal of Hydrology. Elsevier. <www.elsevier.com/locate/jhydrol>. Web. May 27, 2015.

Coulthard, Tom J., Greg R. Hancock and John B.C. Lowry. “Modeling Soil Erosion with a Downscaled Landscape Evolution Model.” Earth Surf. Process. Landforms Earth Surface Processes and Landforms. 37. 10 (2012): 1046-055. Web. May 27, 2015.

Dunham, G., Tisue, S. and Wilensky, U.. NetLogo Erosion model. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL. <http://ccl.northwestern.edu/netlogo/models/Erosion>. 2004.

Grimm, V., Polhill, G., & Touza, J. (2013). “Documenting social simulation models: the ODD protocol as a standard.” In Simulating Social Complexity (pp. 117-133). Springer-Verlag Berlin Heidelberg. 2013.

Kimbrell, Andrew. Fatal Harvest: The Tragedy of Industrial Agriculture. Washington: Published by the Foundation for Deep Ecology in Arrangement with Island, 2002. Print.

Langbein, W. B., and S. A. Schumm. “Yield of Sediment in Relation to Mean Annual Precipitation.” Transactions, American Geophysical Union Trans. AGU 39.6 (1958): 1076. AGU. <http://onlinelibrary.wiley.com/doi/10.1029/TR039i006p01076/abstract>. Web. 5 June 2015.

“Natural Resources Conservation Service (NRCS).” Soil Erosion on Cropland 2007. United States Department of Agriculture (USDA). <http://www.nrcs.usda.gov/wps/portal/nrcs/detail/mo/programs/?cid=stelprdb1041887>. Web. May 23, 2015.

Stroosnijder, Leo. “Measurement of Erosion: Is it Possible?” Catena 64 (2005): 162 – 73. Wageningen University, Erosion & Soil and Water Conservation Group. <www.elsevier.com/locate/catena>. Web. May 22, 2015.

Tramblay, yves, Christophe Bouvier, Claude Martin, Jean-Francois Didon-Lescot, Dragana Todorovik and Jean-Marc Domergue. “Assessment of Initial Soil Moisture Conditions for Event-based Rainfall-runoff Modelling.” Journal of Hydrology 387 (2009): 176 – 87. Journal of Hydrology. Elsevier. <www.elsevier.com/locate/jhydrol>. Web. May 27, 2015.

Wilensky, U.. NetLogo. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL. <http://ccl.northwestern.edu/netlogo/>. 1999.

Tucker, Gregory E., and Rudy Slingerland. “Drainage Basin Responses to Climate Change.” Water Resour. Res. Water Resources Research 33.8 (1997): 2031-047. Water Resources Research. Wiley. <http://onlinelibrary.wiley.com.libproxy.uoregon.edu/doi/10.1029/97WR00409/pdf>. Web. June 1, 2015.

Wilkinson, S. N., C. Dougall, A. E. Kinsey-Henderson, R.D. Searle, R. J. Ellis and R. Bartley. “Development of a Tima-stepping Sediment Budget Model for Assessing Land Use Impacts in Large River Basins.” Science of The Total Environment 468 – 469 (2014): 1210 – 224. Elsevier. <www.elsevier.com/locate/scitotenv>. Web. May 31, 2015.