Predicting location and length of ephemeral gullies with a process-based Topographic Index model

Ephemeral gully (EG) erosion can be a major source of sediment in agricultural watersheds and predicting the location and length of EGs is important to assess sediment source areas. Topographic index (TI) models utilize topographic-only factors to locate concentrated flow paths on agricultural fields. In this paper, the TI-model is expanded by incorporating an overland flow equation that describes surface excess rate, roughness, and soil critical shear stress. An un-calibrated Process-Based (PB) topographic index model was applied to an agricultural watershed in Central Kansas and the results were compared with a calibrated Slope-Area (SA) model. The analysis showed that the PB model slightly over-predicted topographic properties of ephemeral gullies compared to the observed data. Accuracy of the SA model depended on the selected threshold value.


INTRODUCTION
The National Resource Inventory on soil erosion from cropland (NRI, 2007) has reported a 43% decrease in soil erosion in the United States since 1982.This reduction has been mainly due to the adoption of various soil conservation practices; however, the impact on ephemeral gully (EG) erosion is still unclear.EGs are concentrated flow channels in agricultural fields that form when surface water runoff energy exceeds critical shear stress of soil.Such channels are defined as ephemeral channels because they can be temporarily removed by tillage or grow into permanent gullies if not mitigated.Studies by Foster (1986), Poesen et al. (1996), Vandaele et al. (1996), Casalí et al. (1999) and Douglas-Mankin et al. (2011) report that soil loss due to EG erosion can contribute from 30% to as high as 83% of total erosion worldwide.
The main factors that affect the formation and development of EGs include rainfall amount and intensity, field topographic features, soil properties, land cover and management practices.EG initiation is known to be a threshold phenomenon (Knapen and Poesen, 2010).A precipitation event that forces EG erosion normally exceeds a threshold precipitation amount, but also depends on soil conditions and the initial moisture content (Capra et al., 2009).For example, Casalí et al. (1999) reported that within a three-year period from October 1994 to September 1997, only three rainfall events were able to promote EGs.Topographic attributes such as upstream drainage area, slope, and plan of curvature are key topographic controls in the formation process (Foster, 1986;Desmet et al., 1999;Knapen and Poesen, 2010;Daggupati et al., 2013).Soil properties play an important role in allowing soil particles to detach and travel with overland flow.The erosion resistance of topsoil is commonly referred to as soil critical shear stress, the value at which shearing forces of concentrated water flow initiate soil erosion.Soil moisture content can affect the soil critical stress, thus influence formation and progression of EG (Nearing et al., 1989).Land cover and management practices greatly influence vegetation cover and the ability of land to slow the overland flow (Poesen et al., 2011).Vandekerckhove et al. (2000) reported that land cover has a greater influence than climatic conditions in explaining topographic thresholds for different areas.
Several empirical and physically-based models have been developed to quantify EG erosion at both field and watershed scales.The physically-based models, for example, the Ephemeral Gully Erosion Model (EGEM) (Woodward, 1999) and the Revised EGEM model (Gordon et al., 2007), require a wide set of input physical parameters that are difficult to validate at large scales.Accordingly, simple empirical-and regression-based models were also developed (Nachtergaele et al., 2001;Capra et al., 2005;Daggupati et al., 2013).Empirical models utilize a limited set of input parameters.They normally consider only topographic indices and neglect physical processes.The goal of this study was to explore a mixed approach by combining a physically-based formulation to develop a Topographic Index (PB-TI) model.This approach uses GIS to identify, locate, and statistically assess EGs in an agricultural watershed, and evaluates performance of the model when compared to the Slope-Area (SA-TI) model.

STUDY AREA
The Goose Creek watershed (12-digit Hydrologic Unit Code [HUC] 110300140204) is a 13 306 ha watershed that drains into North Fork river in Reno and Kingman counties, northwest of the City of Wichita in central Kansas, USA (Fig. 1).Land use in the watershed is dominated by cropland (64%) and rangeland (29%).Soils are predominantly silt-loams (KDHE, 2000) and the topography is generally flat with a median slope of 0.9%.The major crop type is winter wheat followed by grain sorghum, corn, and soybeans.A weather station is located near the city of Pretty Prairie, Kansas (NCDC USC00146573), east of the watershed.Following the procedure outlined in Daggupati et al. (2013), EGs were digitized using aerial photographs (Google Earth, 2010) and Digital Elevation Models (DEM; USDA-NRCS, 2011).The hill-shade raster files were used to enhance surface relief in ArcGIS (ESRI, 2011).The 874 EG fields identified ranged in length from 14 to 819 m.A field survey was conducted to verify more than 65% of the identified EGs.

SLOPE-AREA TOPOGRAPHIC INDEX MODEL
The Slope-Area (SA-TI) TI model uses two primary topographic attributes, slope S (%) and upstream drainage area A (m 2 ) to calculate an index value TSA at each point in a field (Moore et al., 1988): A critical threshold value of TSA determined a priori for the entire field is compared with the index in equation ( 1) at each point in a field, and EGs were identified if the threshold is exceeded.All points within a field where the threshold is exceeded compose a flowpath of an EG.Daggupati et al. (2013) compared the SA-TI model with three other TI models and reported that the SA-TI model statistically outperformed the others with an optimal threshold of 50.

PROCESS-BASED TOPOGRAPHIC INDEX MODEL
Most topographic index models account for topographic factors such as drainage area, local slope, curvature.However, they often fail to include other characteristics related to soil properties and initial soil moisture conditions, land cover and management as well as precipitation, runoff and peak discharge.Montgomery and Dietrich (1994) proposed an approach that assumes steady-state precipitation intensity over a surface with uniform infiltration capacity and estimates a critical drainage area at any point in the field for channel erosion initiated by overland flow using Manning's equation for flow velocity: Acr=n (ρwg) 5/3 τcr -5/3 E S 7/6  (2) where Acr (m 2 /m) is the critical drainage area divided by the unit contour length, n is Manning's roughness coefficient (s/m 1/3 ), S is local surface slope (m/m), τcr (N/m 2 ) is the critical shear stress, ρw is density of water (kg/m 3 ), g is acceleration due to gravity (m/s 2 ), and E (m/s) is a steady-state precipitation excess rate.The unit contour length was taken as a resolution scale of DEM.This process-based TI (PB-TI) model in equation ( 2) is based on the assumption that EGs form by turbulent overland flow when the stream power exceeds critical shear stress.At each point in the field, an EG is assumed to occur if the local drainage area AS is greater than Acr, This approach provides additional insight in EG identification by accounting for topographic features (S and A), land cover condition (n), weather and hydrology (E), and soil property (τcr).Topographic attributes can be determined from a DEM, the roughness coefficient can be found from land use datasets and critical shear stress is available from online soil databases.However, the precipitation excess rate (E) is event-specific and it must be modelled separately.The antecedent soil moisture content preceding a runoff event is not explicitly included in the formulation but assumed to be accounted in the simulation of E.
To calculate E, the Soil and Water Assessment Tool (SWAT) (Arnold et al., 1998) was applied to the studied area.SWAT is a continuous, process-based model that calculates surface runoff for each field based on the SCS Curve Number method (USDA-SCS, 1972).Watershed input geospatial datasets were 10 m × 10 m DEM (USGS, 1999), SSURGO soil dataset (USDA-NRCS, 2005;Sheshukov et al., 2011), field-based cropland NASS dataset and management schedule data (Gali et al., 2012).The SWAT model was run from 1 January 1996 to 31 December 2008 and calibrated for streamflow.Average daily runoff rates E were extracted for each catchment during the storm event on 3 October 2002 that was assumed to cause EG erosion.

RESULTS
Two TI models (PB-TI and SA-TI) were applied to cropland fields of Goose Creek watershed.Catchments with EGs were identified and the length of each EG was calculated.A map of agricultural fields and EG lines are shown in Fig. 1.
Both the occurrence and length of EGs were predicted by two models and analysed using a statistical error matrix method (Gómez-Gutiérrez et al., 2009;Meyer and Martinez-Casosnovas, 1999).The method is based on treating the presence of EG as unity and absence as zero.The false positive is defined as predicted/not observed and the false negative as not predicted but observed.The occurrence statistic kappa also used in the analysis; it reaches a higher value with better model performance.The resulting statistics are presented in Table 1.The PB-TI model was found to have a greater false positive rate (34%), but lower false negative rate (18%) compared to the SA-TI model.Both statistics represent better model performance for lower values, so neither model had a clear advantage.The kappa statistics were similar for the two models, but the purely topographic based SA-TI model exhibited slightly better overall accuracy of EG occurrence prediction.
EG length was calculated for the two models and compared for each catchment with EGs digitized from the aerial image.The length of predicted EGs ranged from 50 m to 1253 m for the PB-TI model and from 50 m to 898 m for the SA-TI model.There were 714 matching EG catchments for the PB-TI model and 610 for the SA-TI model.
Figure 2 presents the number of EG catchments for each 20-m increment of EG length.The maximum number of EGs in the watershed was 83 for EGs with a length between 80 and 100.The number of EGs was low for short gullies, approached a maximum at about 100 m, and then decreased with increase of length, reaching single digits for gullies longer than 600 m.Both models underestimated the number of EGs that were shorter than 200 m, but over predicted the number of longer EGs.The statistics (R 2 , NSE, and PBIAS) that evaluate performance of the two models are shown in Table 1.Probability of exceedence of EG lengths for the observed and modelled data are presented in Fig. 3.Over the entire EG length range, both models exhibited overprediction of EG length, as shown by negative PBIAS values in Table 1 and the curves in Fig. 3.The three calculated statistics were better for the SA-TI model (Table 1), which was also supported by Figs 2 and 3 that show better agreement with the observed curve than the PB-TI model.The observed EG length for 50% probability was 110 m, while it was 130 m for the SA-TI model and 165 m for the PB-TI model.
Better performance of the SA-TI model likely resulted due to model's prior calibration in a small area of the Goose Creek watershed and identification of the optimal threshold of 50 prior to application to the entire watershed.Initial calibration of the SA-TI model in a small area within a watershed can result in reasonable model performance when applied to a topographically similar region or larger watershed.In contrast, the PB-TI model was not calibrated, but parameters in equation ( 2), such as critical shear stress and roughness coefficient, were taken from their respective databases without adjustment.However, the excess rate E was carefully adjusted based on comprehensive continuous hydrologic model.
A comparison of the two TI models showed that performance of the calibrated SA-TI model was only slightly better that the one of un-calibrated process-based PB-TI model.It demonstrates a robustness of the process-based model and its capabilities to produce reasonable EG length estimation without calibration.The PB-TI model included a steady-state surface excess rate that is difficult to find without the use of a hydrologic model.Further model improvements may be possible by reducing uncertainties with runoff estimates.Future work will focus on development of a dynamic TI model that would incorporate time-dependent surface runoff estimates that could provide more realistic EG location, length, and dimension estimations.The ongoing EG monitoring study will provide additional insight on site-specific physical properties and improve model performance in predicting location and length of EGs, which is essential for accurate estimation of EG sediment erosion rates.

Fig. 1
Fig. 1 Map of Goose Creek watershed in central Kansas, USA, showing cultivated cropland fields (dark areas), ephemeral gullies (short solid lines), stream network (black lines), and subwatersheds.The period from 2002 to 2003 was selected for assessment of EG locations and length because of previous knowledge of field management, topography, aerial imagery, precipitation and streamflow data.Crop tillage and planting for winter wheat was completed in September 2002.A Digital Ortho Quarter Quadrangle aerial image of Goose Creek watershed from 22 March 2003 was acquired (Google Earth, 2010) and it showed a network of developed EGs.Analysis of precipitation records revealed that EGs most likely resulted from precipitation events from 2 October to 4 October 2002, because no significant storm event was recorded during winter of 2002 and spring of 2003.Following the procedure outlined inDaggupati et al. (2013), EGs were digitized using aerial photographs (Google Earth, 2010) and Digital Elevation Models (DEM; USDA-NRCS, 2011).The hill-shade raster files were used to enhance surface relief in ArcGIS(ESRI, 2011).The 874 EG fields identified ranged in length from 14 to 819 m.A field survey was conducted to verify more than 65% of the identified EGs.

Fig. 2
Fig. 2 Number of EG catchments for each 20 m range of EG length for observed EGs, PB-TI and SA-TI models.

Table 1
EG occurrence (rows one to three) and length (rows four to six) statistics for observed, PB-TI and SA-TI models.