Distributed soil loss estimation system including ephemeral gully development and tillage erosion

A new modelling system is being developed to provide spatially-distributed runoff and soil erosion predictions for conservation planning that integrates the 2D grid-based variant of the Revised Universal Soil Loss Equation, version 2 model (RUSLER), the Ephemeral Gully Erosion Estimator (EphGEE), and the Tillage Erosion and Landscape Evolution Model (TELEM). Digital representations of the area of interest (field, farm or entire watershed) are created using high-resolution topography and data retrieved from established databases of soil properties, climate, and agricultural operations. The system utilizes a library of processing tools (LibRaster) to deduce surface drainage from topography, determine the location of potential ephemeral gullies, and subdivide the study area into catchments for calculations of runoff and sheet-and-rill erosion using RUSLER. EphGEE computes gully evolution based on local soil erodibility and flow and sediment transport conditions. Annual tillage-induced morphological changes are computed separately by TELEM.


INTRODUCTION
The Revised Universal Soil Loss Equation Version 2 -RUSLE2 (ARS, 2008;Renard et al., 2011) computes sheet and rill erosion along one-dimensional (1D) hillslope profiles, from the top of the hill, where runoff begins, to a location where runoff meets a concentrated flow channel.RUSLE2 is an advancement of the erosion prediction technology of RUSLE (Renard et al., 1997) as it implements sediment transport methods that permit the determination of sediment deposition that occurs in areas of reduced slope steepness frequently found in the concave areas.Application of RUSLE2 involves the selection of representative 1D profiles in the landscape, which requires judgment and expertise to apply the model correctly for conservation planning.Another remaining deficiency of RUSLE2 is its inability to estimate erosion within the channels that end hillslope flow paths, the locations where ephemeral gullies may form.
The combination of terrain analysis algorithms with high-resolution topographic data, such as that available through use of Light Detection and Ranging (LiDAR) and similar technologies, makes it possible define overland flow paths that represent runoff accumulation and convergence created by the field topography, eliminating the need to select 1D profiles for estimation of erosion, and at the same time providing a consistent method to define the location of potential ephemeral gullies.
The primary objective of this paper is to describe the linkage of a distributed version of RUSLE2 (RUSLER, short for RUSLE2-Raster), with a tillage erosion model, TELEM (Tillage Erosion and Landscape Evolution Model; Vieira andDabney, 2009, 2011), and a new ephemeral gully model, EphGEE (Ephemeral Gully Erosion Estimator), to provide high-resolution, spatially distributed estimations of sheet, rill, and ephemeral gully erosion.In a companion paper (Dabney et al., 2014) we illustrate this technology through application to a research watershed located near Treynor, Iowa, USA.

GRID-BASED RUSLE2
RUSLER re-uses RUSLE2's methods for 1D hillslope profiles to compute erosion, sediment transport, and deposition on each cell of a two-dimensional (2D) computational grid representing entire agricultural fields.Using terrain elevation data from high-resolution digital elevation models doi:10.5194/piahs-367-80-2015(DEMs) generated from LiDAR, and soil type and agricultural management spatial distributions defined in matching raster grids, RUSLER utilizes a series of algorithms built into a software library called LibRaster to define drainage-based areas for the calculation of sheet and rill erosion by RUSLE2.First, overland flow paths are defined using flow directions assigned to each grid cell in the DEM according to the direction of steepest slope using the classical D-8 algorithm.The flow directions are then used to determine drainage paths and flow accumulation; grid cells with upslope contributing area exceeding a threshold are considered "channels," locations where RUSLE2 hillslopes end and where ephemeral gullies may form (Dabney et al., 2014).For each channel cell, LibRaster algorithms determine left and right bank "channel cell catchments" that contribute runoff and sediment to the channel system on an event basis.In the RUSLE2 computation, each 2D cell is equivalent to a 1D slope segment, as illustrated in Fig. 1.This approach allows the immediate application of the RUSLE2 computational engine in a 2D-gridded terrain description.To support the application of RUSLE2 in complex 2D landscapes with flow convergence, the way RUSLE2 estimates slope length was modified in the 2014 release (USDA-ARS, 2014).Rather than being calculated based on linear distance from the top of the hillslope, slope length is now calculated based on the ratio of runoff leaving a segment or cell to that of runoff generated within that cell.This approach accounts for spatial variability in runoff generation related to soil type and management effects.Further, RUSLE2 was enhanced to generate a representative series of runoff events (Dabney et al., 2011(Dabney et al., , 2012)).The new methods allow the prediction of daily hillslope runoff and sediment yields suitable for linkage with a physically-based ephemeral gully model.RUSLER's channel configuration, storm event data, and computed runoff and sediment loads are exported as a set of XML files for use by the ephemeral gully calculator EphGEE.

EPHEMERAL GULLY EROSION ESTIMATOR (EPHGEE)
EphGEE simulates ephemeral gully erosion through a procedure that takes into account detachment of soil through the shear force of flowing water, sediment transport capacity, and changing channel dimensions due to erosion and deposition.EphGEE computes ephemeral gully evolution for a series of storm events over a network of potential gully channels, which are either generated through terrain analysis algorithms based on high-resolution elevation data, or using a channel configuration specified by the model user.Gully dimensions (depth and width) change with each runoff event and have soil added by each RUSLE2 tillage operation.

Channel evolution concept
In EphGEE, channels evolve following Foster and Lane's (1983) approach, which assumes gully channels erode in a two-step process, first downcutting to form a rectangular channel whose width is proportional to flow and soil parameters, and once a layer of high resistance is reached, channels widen at decreasing rates, asymptotically approaching the width where shear stress at the toe of the channel bank is equal to the specified critical stress.The equations that describe change in channel dimensions of a reach are similar to those used in channel erosion models of CREAMS and WEPP.Haan et al. (1994) provide a clear conceptual derivation of this channel erosion theory that is based on several assumptions: (1) that Manning's equation applies, (2) that the shear stress distribution around the cross-section of a channel can be represented by a dimensionless distribution, (3) that the soil consists of a uniform erodible layer with characteristic erodibility and critical shear stress values overlying a non-erodible layer at a specified depth, (4) that potential detachment rate is proportional to excess shear stress, (5) that actual detachment is proportional to the unsatisfied transport capacity of a steady-state runoff rate, (6) that transport capacity can be determined by the set of equations proposed by Yalin (1963), and ( 7) that deposition occurs if sediment load exceeds transport capacity.

RUSLER -EphGEE integration
In EphGEE, this conceptual model was implemented as a C++ program, which has as inputs a map of potential gullies, soil properties, sediment size definition, and spatially-distributed runoff and sediment loads for a series of storm events.When linked with the grid-based RUSLE2, all the data are supplied as a set of XML files created by RUSLER.For gully erosion computations, a channel network is created based on the one used by RUSLER, but with node spacing that can be adjusted by the user, and that is typically smaller than the resolution used for the distributed RUSLE2 calculations.Figure 2 illustrates the RUSLER and EphGEE channel networks.At each node, the channel geometry is defined by a series of distance-elevation points which define the crosssectional shape of the gully.This initial shape can be derived from high-resolution terrain data, or a simplified shape such as a triangular cross-section may be assumed.For each storm event, runoff from all RUSLE2 catchments is converted into steady, spatially varying flow discharges distributed among all nodes in the channels.Water surface profiles are computed using the Standard-Step Method (Henderson, 1966) using the flow discharge, hydraulic roughness, and channel geometry and steepness at each node.Channel roughness is determined by RUSLE2 as the result of local grain and form roughness and ground cover conditions, and varies seasonally due to mechanical disturbances and the natural roughness decay with time.In EphGEE, roughness values are also adjusted based on the amount of erosion or deposition occurring at a channel cross-section.Roughness values vary linearly between the RUSLE2 value for a noneroded channel and user-specified roughness values for channels with substantial erosion or deposition.At channel junctions, it is assumed that the thalweg elevation is initially the same at all branches.Flow calculations assume that the water surface level is the same for all nodes sharing the same location (Fig 1b).If however, after erosion takes place, one channel branch is significantly lower than the other, flow in a tributary with a higher thalweg elevation may be controlled by a free-fall condition (critical flow at the brink).Similarly, depending on flow conditions and channel geometries, high flow or significant deposition in a channel may induce backwater in a tributary.

Erosion and sediment transport
With the flow properties known for all nodes in the channel network, erosion is then computed based on erosion rate determined as a function of local shear stress, (Ascough et al., 1997).Shear stress available for sediment transport, τg, is calculated from the total shear stress, the total Manning coefficient, nt, and an assumed grain Manning n value of 0.01: (1) The potential detachment capacity Dp (kg m -2 s -1 ) is given by: where K is the soil erodibility (s m -1 ), τg is the shear stress available for sediment transport (Pa), τc is the soil critical shear stress (Pa).Actual detachment, D, occurs at a rate that will just fill transport capacity, computed by: where g is the total sediment load, including any lateral inflows, and tc is the transport capacity (both in kg m -2 s -1 ).Detachment is assumed to be non-selective, that is, sediments from all size classes are entrained and available for transport.Channel geometry changes for each event are calculated according to the Foster and Lane's (1983) procedure, as explained above.
As equations ( 2) and (3) show, the onset of erosion and its rate of development are controlled by a critical shear stress value and an erodibility coefficient of the erodible soil layer.In the absence of a national database providing these values, soil critical shear stress and erodibility estimates were approximated as a function of soil clay percentage and time since tillage using relationships derived with data in Watson et al. (1986) for seedbed, crop, and no-till conditions: where τc is the soil critical shear stress (Pa), K is the soil erodibility (s m -1 ), C is the soil clay percentage, t is the time since last tillage or sediment deposition date (d), and a and b are coefficients, and b = 0.0152 if t >1000 d.These crude relations produce results consistent with the observations of Franti et al. (1999) for two no-till Iowa soils with about 24% clay (τc = 5.6 Pa, K = 0.004 s m -1 ) but underestimate observed critical shear stress (τc = 3.6 Pa, K = 0.023 s m -1 ) for tilled soils (t < 10 d ).
Grain size distributions used with predictive transport functions are defined as sediment bins based on fall velocities.Because the size and fall velocities of sediment aggregates varies with clay content of the eroding soil, both RUSLE2 and EphGEE place sediment in a number of predefined bins.Sediment distribution from sheet and rill erosion is transferred from RUSLE2 to EphGEE as part of the storm events data.EphGEE is being developed to provide a choice of sediment transport capacity equations.The simplified relationship that approximates the procedures of Yalin (1963) presented by Ferro (1998) has been recently applied successfully (Dabney et al., 2014).Transport capacity may be partitioned among the several size classes as a function of local size distribution of transport sediment (fractional loads), or optionally, based on the finding of Langendoen et al. (2014), fractional transport capacity should be applied only to the two coarsest sediment fractions predicted by RUSLE2 ("sand" and "large aggregate") while the finer "silt," "clay," and "small aggregates" fractions are considered as washload and assumed to be transported without limiting the transport capacity available to coarser sediment fractions.

Channel deposition
If at any point the total sediment load in the channel exceeds transport capacity, deposition is determined independently for each sediment size class, at a rate controlled by the fall velocity of grains and aggregates.Therefore, deposition is selective, with coarser sediment fractions depositing at faster rates, while finer particles are preferentially transported downstream.Because storm events are simulated as steady flows, and the channel geometry is updated only at the end of the event, it is possible that excessive deposition occurs in areas where significant reduction in transport capacity occurs.Sudden decreases in channel slope steepness or backwater regions often cause excessive deposition, especially when flow depths and channel dimensions are small.If such situation is detected, the model limits the amount of deposition at that location so that the resulting cross-sectional deposition area is less than half the flow area for that storm.Flow leaving that location becomes overloaded with sediment, which is distributed and deposited further down the channel, avoiding unrealistic amounts of deposition in transitional regions.This simplified approach is analogous to the assumption of non-equilibrium used in stream sediment transport models (Wu, 2008), where the assumption that sediment transport rate equals transport capacity is relaxed.

TILLAGE EROSION
Tillage erosion is determined separately with TELEM, which computes soil translocation by each tillage operation in the direction of implement travel and perpendicular to that direction as a function of the terrain slopes in those directions (Vieira andDabney 2009, 2011).A version of the model has been recently implemented into LibRaster to facilitate its combined used with RUSLER.Tillage translocation rates are calculated based on implement tillage translocation coefficients, which are available in the RUSLE2 management database.Tillage directions and the presence of tillage boundaries such as vegetative strips, fences, ditches and other barriers may be specified as raster maps.Values of computed soil loss (Mg ha -1 year -1 ) resulting from changes in translocation due to varying terrain slopes or induced by the boundaries are converted to terrain elevation changes using an assumed bulk density value.TELEM is capable of estimating sediment delivered directly to the channels by tillage.When tillage erosion calculations are not performed, gully channels can be assumed to be reset to a predefined cross section shape after each tillage operation.

FIELD SITE FOR VALIDATION
The 2D version of RUSLE2 was applied and compared with observations on Watershed 11 (Rachman et al., 2008) of the USDA-ARS Deep Loess Research Station located near Treynor, Iowa (Karlen et al., 2009).This 6.3 ha watershed, which has an extensive research archive, had daily runoff and sediment yield monitored at the watershed outlet from 1975 to 1991.Throughout the period of record, a grassed waterway was located in the lower portion of the watershed (Dabney et al., 2014).The field was farmed with contour-planted, conventional-tilled (CT) corn (Zea mays, L.).
A detailed account of the field validation tests can be found in Dabney et al. (2014).A DEM was created at 3-m resolution, and RUSLE2 simulations were conducted for the scenario for growing spring-ploughed (CT) corn yielding 7.6 Mg ha -1 with no additional conservation practices.The climate and soils records for Pottawattamie County, Iowa, in the standard RUSLE2 database were used in all simulations (USDA-NRCS, 2014).Event runoff and sediment delivered to the channels by RUSLE2 were used as input to EphGEE.EphGEE simulations were run with and without a grassed waterway in the last 200 m of the watershed's primary channel.Hydraulic roughness and particle size distributions and properties were calculated as described by Dabney et al. (2014).Tillage operations reset channel dimensions prior to runoff events on the tillage dates specified in the RUSLE2 management database records.
The RUSLE2-predicted average annual runoff of 67 mm was higher than the observed 50 mm runoff, possibly because the RUSLE2 simulation did not include a representation of the grassed waterway, which slows down runoff and increases infiltration.The monthly pattern of observed and modelled runoff showed general as did the size of runoff events with varying return periods (Dabney et al., 2012).Average RUSLE2 predicted sediment yield from the hillslopes to the channel system averaged 35 Mg ha -1 year -1 (Dabney et al., 2013), whereas observed sediment yield at the watershed outlet averaged only 14.6 Mg ha -1 year -1 .Ephemeral gully predictions showed gully development in the mid-and upper parts of the watershed where slope steepness is greater, and significant deposition at the lower, flatter part where a grassed waterway was in place.The difference between hillslope erosion predictions and observed sediment delivery suggests that at this site and during this period, waterway sediment deposition exceeded ephemeral gully detachment.

CONCLUSIONS
The approach presented here permits the estimation of distributed runoff and soil loss using highresolution topography, accounting for actual drainage patterns and the resulting flow convergence, considering the spatial variation of soils and vegetation, and considering agricultural operations such as tillage.The computational methods allowed the RUSLE2 engine, used in classical 1D hillslope soil loss estimation, to be also used in 2D, grid-based configuration.
The GIS-inspired, raster-based drainage analysis provides a consistent method for the identification of complex hillslopes for calculation of sheet and rill erosion, and at the same time allows defining of where ephemeral gullies may form.The integrated application of RUSLER and EphGEE provides a mechanism for the estimation of runoff and sediment loads that control the development of ephemeral gullies.Enhancements introduced to the 2014 release of RUSLE2 permit the determination of a set of storm events that represent geographical and seasonal climatic properties, while also accounting for soils and management characteristics in the generation of runoff and soil loss.The availability for spatially-variable storm data permits the prediction of gully evolution using a physically-based model (EphGEE), where erosion or deposition is determined based on flow and sediment transport controlled by local channel geometry and soil conditions.
The relatively close match between predicted and observed watershed sediment yield is encouraging and suggests that the methods, approximations, and simplifications employed may allow distributed assessment of soil degradation and water quality impacts from agricultural management.Further testing and development will be needed to clarify the ability of this technology to approximate observations over a wider range of climate, soil, and management systems.

Fig. 1
Fig. 1 (A) Illustration of 2D "channel cell catchments" (identified in different shades), used for RUSLER sheet-and-rill calculations and to determine loadings to ephemeral gully channels.(B) Linkage between RUSLER 2D cells and the equivalent RUSLE2 1D profile segment, showing how runoff and sediment loads entering each cell are determined.

Fig. 2
Fig. 2 (A) A raster channel network used for RUSLER 2D simulation (shown as square cells with shades identifying different channel links), which is converted into a network for computation of gully evolution with EphGEE, where circles indicate nodes where channel flow, erosion, and sediment transport are calculated.Nodes with a black outline match the locations of the raster channel cells; remaining nodes are added to increase computational accuracy.(B) Illustration of a channel junction, where nodes from three channel links meet at the same location.