Next Article in Journal
Best-Fit Probability Models for Maximum Monthly Rainfall in Bangladesh Using Gaussian Mixture Distributions
Next Article in Special Issue
The Influence of Crude Oil on Mechanistic Detachment Rate Parameters
Previous Article in Journal
Geoheritage, Geotourism and the Cultural Landscape: Enhancing the Visitor Experience and Promoting Geoconservation
Previous Article in Special Issue
Effects of pH-Induced Changes in Soil Physical Characteristics on the Development of Soil Water Erosion
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Geosciences
Volume 8
Issue 4
10.3390/geosciences8040137
geosciences-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modelling Ephemeral Gully Erosion from Unpaved Urban Roads: Equifinality and Implications for Scenario Analysis

by
Napoleon Gudino-Elizondo
1,2,*,
Trent W. Biggs
2,
Ronald L. Bingner
3,
Yongping Yuan
4,
Eddy J. Langendoen
3,
Kristine T. Taniguchi
2,
Thomas Kretzschmar
1,
Encarnacion V. Taguas
5 and
Douglas Liden
6
1
Departamento de Geología, Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE), 22860 Ensenada, Mexico
2
Department of Geography, San Diego State University, San Diego, CA 92182, USA
3
National Sedimentation Laboratory, Agricultural Research Service, USDA, Oxford, MS 38655, USA
4
Office of Research and Development, United States Environmental Protection Agency (USEPA), Research Triangle Park, NC 27711, USA
5
Department of Rural Engineering, University of Córdoba, 14071 Cordoba, Spain
6
San Diego Border Office, United States Environmental Protection Agency (USEPA), San Diego, CA 92101, USA
*
Author to whom correspondence should be addressed.
Geosciences 2018, 8(4), 137; https://doi.org/10.3390/geosciences8040137
Submission received: 26 February 2018 / Revised: 7 April 2018 / Accepted: 7 April 2018 / Published: 17 April 2018
(This article belongs to the Special Issue Soil Hydrology and Erosion)
Browse Figures
Versions Notes

Abstract

:
Modelling gully erosion in urban areas is challenging due to difficulties with equifinality and parameter identification, which complicates quantification of management impacts on runoff and sediment production. We calibrated a model (AnnAGNPS) of an ephemeral gully network that formed on unpaved roads following a storm event in an urban watershed (0.2 km2) in Tijuana, Mexico. Latin hypercube sampling was used to create 500 parameter ensembles. Modelled sediment load was most sensitive to the Soil Conservation Service (SCS) curve number, tillage depth (TD), and critical shear stress (τc). Twenty-one parameter ensembles gave acceptable error (behavioural models), though changes in parameters governing runoff generation (SCS curve number, Manning’s n) were compensated by changes in parameters describing soil properties (TD, τc), resulting in uncertainty in the optimal parameter values. The most suitable parameter combinations or “behavioural models” were used to evaluate uncertainty under management scenarios. Paving the roads increased runoff by 146–227%, increased peak discharge by 178–575%, and decreased sediment load by 90–94% depending on the ensemble. The method can be used in other watersheds to simulate runoff and gully erosion, to quantify the uncertainty of model-estimated impacts of management activities on runoff and erosion, and to suggest critical field measurements to reduce uncertainties in complex urban environments.

1. Introduction

Both rural and urban development can increase erosion and the delivery of land-based sediment into receiving water bodies, including estuaries, coasts, and inland lakes and reservoirs. Unpaved roads, in particular, represent one of the principal landscape features of rural urbanization in developing countries. Ephemeral gullying, including on unpaved roads, is an important soil erosion process reported in many environments [1]. Road drainage impacts erosive processes, increasing flow peaks and total discharge [2], which is also observed in monsoonal climates [3].
Ephemeral gullies are important components of sediment budgets in both natural and human-disturbed environments. The term ephemeral indicates that they are temporary features, commonly removed by tillage operations [4] or filled by road maintenance in urban environments. Ephemeral gully formation is the product of a complex interaction between terrain topography, climate, soil properties, land cover, and management practices [5], and ephemeral gullies can be the primary source of sediment loss in agricultural and urban environments [6,7,8].
Semi-arid watersheds are highly sensitive to soil and stream channel erosion following urbanization [9,10]. Gudino-Elizondo et al. [11] observed that gullies in Tijuana, Mexico, formed almost exclusively on unpaved roads, reflecting the role of roads in routing flow and their vulnerability to gully erosion. Unpaved roads can also be an important component of anthropogenic sediment generation in the study area [11], as has also been reported in other settings [12], including logged forests [13,14] and tropical islands [15].
There are few studies assessing gully erosion in urban settings, as documented in a compilation of gully erosion studies by Castillo & Gomez [16]. Adejiji et al. [8] described the relationship between urban surface characteristics and gully erosion in Nigeria, and found a significant relationship between soil texture, land use, and gully erosion. However, measurements and modelling of ephemeral gully erosion rates on unpaved roads have rarely been conducted in urban watersheds. Control of gully erosion could involve road paving, but at the cost of increasing peak discharge. Other best-management practices (BMPs) include revegetation of hillslopes that produce runoff, which could mitigate both runoff and sediment production, but this strategy has not been quantitatively evaluated. The trade-off between sediment control and runoff production is particularly important, but remains unquantified.
Numerical models have been used to simulate soil and gully erosion rates [17,18,19,20] and to assess the impacts of conservation practices [21]. These models differ in terms of their structure, their assumptions and the input data necessary for model calibration and application [22]. Bull and Kirkby [23] reviewed the conditions for gully formation and noted that gully modelling must be based on the relationship between flow hydraulics and soil properties [24]. Nachtergaele et al. [25] reported a good performance of the Ephemeral Gully Erosion Model (EGEM) predicting gully volumes in agricultural areas of Spain and Portugal.
The Annualized AGricultural Non-Point Source (AnnAGNPS) model was developed to simulate sheet and rill erosion in agricultural environments [19,20], and has been utilized and validated in many studies, including in evaluations of the impact of agricultural BMPs [26,27,28,29,30,31,32]. Gordon et al. [33] improved on the EGEM using more process-based techniques and this revised EGEM has been incorporated in AnnAGNPS [20]. Improvements to the gully widening approach within AnnAGNPS were developed by Bingner et al. [20]. Head-cut migration rates within the model are based on physical approximations of mass, momentum, and energy transfer, described by Alonso et al. [34]. The AnnAGNPS model has not been tested to simulate and monitor ephemeral gully erosion rates in an urban context.
In hydrologic and soil erosion modelling, several parameter sets may adequately simulate the observed behaviour of the system; such models are called “behavioural” [35]. Hornberger and Spears [36] rejected the idea of an optimal model structure or parameter set in favour of multiple parameter combinations, which all provide acceptable fits to observed data, called equifinality by Beven [37]. Equifinality suggests that there are multiple interactions among the parameters within a model to produce simulations that may be equally acceptable. Equifinality is especially important when simulating the impacts of changes in climate, land use, or watershed management, since different parameter ensembles can generate different predictions under change [35]. Field measurements may be taken to constrain model parameters, but those measurements may or may not match the parameters obtained through calibration due to either unsampled heterogeneity, problems with model structure, or to other processes operating at spatial scales larger than that of the field measurements. To our knowledge, no study has addressed equifinality in gully erosion modelling and its impact on scenario analysis, particularly in an urban setting.
This paper aims to generate a set of behavioural gully erosion models in a rapidly urbanizing watershed, and to explore the impact of parameter uncertainty on scenario analysis in a practical management context. We address the following research questions: (a) How well does the AnnAGNPS model predict urban gully erosion? (b) What are the most sensitive AnnAGNPS parameters in urban gully erosion modelling? And (c) What are the implications of parameter uncertainty for evaluation of the impact of road paving and other BMPs on runoff and erosion? The study is novel in terms of evaluating AnnAGNPS’s capabilities in assessing gully erosion in urban watersheds, which included using a high-horizontal-resolution (30 cm cell size) Digital Elevation Model (DEM) generated using a combination of Unmanned Aerial Systems (UAS) and Structure from Motion (SfM) photogrammetric techniques [11] to improve representation of topographic attributes and flow routing to predict ephemeral gully formation. Understanding the process of gully erosion will be critical in describing and quantifying sediment production within urbanized watersheds, and consequent loads of water and sediment to ecosystems downstream.

2. Materials and Methods

2.1. Study Area

The San Bernardo (SB) neighbourhood is located within Los Laureles Canyon Watershed (LLCW), a bi-national watershed that flows from Tijuana, Mexico, into the southwestern arm of the Tijuana River Estuary, Imperial Beach, CA, USA. The LLCW drainage area is 11.58 km2, with 10.8 km2 in Mexico and 0.75 km2 in the United States (Figure 1b). The climate is Mediterranean, with a wet season from November to April and annual precipitation of approximately 240 mm per year. Soils in SB are sandy uplifted marine terraces with steep slopes (mean 15 degrees), resulting in high vulnerability to soil and gully erosion [11]. Based on a soils map of San Diego County and samples of soil texture taken in the watershed, the soils are similar to the Las Flores soil group, which are described as having loamy sand A horizons with greyish brown and light brownish grey colour, and a sandy clay B horizon grading to weakly consolidated siliceous marine sandstone in the C horizon [38]. SB has typical mixed urban-rural land cover (Figure 1a) with high population density (~6500 people·km−2). SB was urbanized in 2002, and has unauthorized housing developments (”invasiones”). The construction of unpaved roads on highly erodible soils enhances gully formation, affecting the quality of life for the residents [39], and is likely a significant contributor to total sediment production at the watershed scale. The gully network in SB is filled in with sediment at specified dates to represent road repair. However, road repair was not important, because gully formation was simulated from a single storm event.
Excessive erosion, transport and deposition of sediment have many detrimental effects on the people living in the watershed (Figure 1c) and have impaired conditions for aquatic life in the Tijuana River Estuary (Figure 1b). The Tijuana River National Estuarine Research Reserve, located in the United States, is listed as “impaired” by the State of California due to excessive sediment loads [40]. Several U.S. government agencies spend approximately $3 million per year to remove sediment produced in Mexico [41], and it is therefore important to quantify soil erosion rates in the upper watershed in order to identify cost-effective solutions to reduce sediment loads into the Estuary.

2.2. Observed Gully Erosion

Both ground- and UAS-based surveys of a gully network that formed in SB following a large storm event on 5–7 January 2016 were conducted on 16 January 2016 (Figure 2). The storm was the largest of the water year (~50 mm of total rainfall) and had a 15 min maximum rainfall intensity of 4.8 mm, which has a 1 year recurrence interval [42]. Other storms occurred during the year, but all were smaller than the threshold precipitation typically required to produce gullies in SB (~25–35 mm), as observed on other field visits following storm events during three hydrological years (2013–2016) [42]. The observed sediment production during this storm event was used to test the performance of the AnnAGNPS model in simulating gully erosion on unpaved roads.
A total of nine sub-watersheds were used to estimate gully erosion rates. Gully perimeters were digitized manually from a UAS-SfM-derived orthophoto, and field measurements were used to assist with visual estimation of the gully depth of each digitized gully in order to calculate Specific Soil Loss (SSL, [11]). We used the orthophoto to interpolate 48 field measurements (Figure 3b) of gully depth. Polygons delineating gully sections with the same depth were created based on the shadows and colours of the section. Gully sections without a nearby field measurement were identified, delineated, and assigned a depth based on the shadows and colour likeness with other gully sections containing field measurements [11]. The volume of gully erosion was calculated as the product of the gully area times the gully depth. The specific soil loss (SSL, which is the average depth of soil loss in the watershed), was then calculated as the total volume of gully erosion (m3) normalized by each drainage area (Ad) (m2, Figure 3). See Gudino-Elizondo et al. [11] for a full description of methods.

2.3. AnnAGNPS Model

The AnnAGNPS model is a distributed-parameter numerical model developed by the Agricultural Research Service (ARS) and the Natural Resources Conservation Service (NRCS) of the US Department of Agriculture (USDA) to simulate water and sediment loads from any source within a watershed on a daily time step [20]. AnnAGNPS has been used to assess watershed response to different conservation practices [20]. The spatial distribution of soils, land use, and terrain attributes is used to discretize the watershed into topographically defined sub-watersheds (AnnAGNPS cells) that are assumed to be homogeneous in land cover and soil type. The homogeneous spatial distribution of soils used in this analysis was based on field observations, visual interpretation of high-resolution satellite imagery in GoogleEarthTM and soil samples taken for texture (N = 4) and jet-erosion tests (N = 8). AnnAGNPS simulates the contribution of different erosion processes, including sheet, rill, gullies, and streambed and bank.
Total runoff is calculated following the SCS curve number method [43]. Peak discharge, time base and the storm type are calculated using methods described in USDA-NRCS Technical Release 55 (TR-55) [44]. A type-II, 24 h rainfall distribution (TR-55) was used, and the type was determined by comparing cumulative rainfall observed at a nearby rain gauge [42] with the cumulative distribution functions from TR-55. Type II is representative of intense rainfall observed during convective events in semi-arid regions of the south-west United States. The model does not distribute the rainfall data over the day (e.g., minute or hourly), but rather uses the storm type distribution (here, type II) to assign regression coefficients that determine the peak discharge as a function of the ratio of initial abstraction to 24 h precipitation that is then used in determining sediment transport. A topography-based method (TopAGNPS) was used to map the location of the most downstream point of the potential ephemeral gullies [5]. This approach automates identification of the location of potential ephemeral gullies based on the comparison of the runoff erosivity estimated from topographic attributes (i.e., local slope and drainage area) with soil properties. The gully erosion model in AnnAGNPS requires a model estimate of the peak discharge at the incision point (head-cut or nickpoint) where gullies form. If the shear stress exerted by the runoff erosivity exceeds the soil critical shear stress, the gully incises. Once the incision reaches a non-erodible soil layer, defined as TD in AnnAGNPS, the nickpoint migrates upslope at a rate dependent on streamflow conditions and soil resistance to erosion [20,21,33]. The gully width was calculated within AnnAGNPS using the Wells’ Equation [45], which was developed in experimental conditions using packed soil beds under similar soil textures as those observed in SB, expressed as:
W = 9.0057 × ( Q p × S ) 0.2963
where W is the gully width (m); Qp the peak discharge at the gully head (m3/s); and S is the average bed slope above the gully head (m/m). Other relationships were investigated for use by AnnAGNPS [46], with the Wells approach providing the best response for gullies that are repaired. Many other empirical relationships have been developed for gullies or channels that were not repaired, but in the watershed for this study, gullies are repaired after precipitation events, and therefore encouraged us to use it for this analysis.
Rainfall intensity and SCS Curve Number (CN) are the most important parameters for the peak discharge and total runoff calculations using the AnnAGNPS model, and both determine the fraction of the rainfall contributing to overland flow. Manning’s roughness coefficient is also an important parameter in runoff production and runoff erosivity.

2.4. Model Setup

The topographic attributes, such as total and individual cell areas, length of channels, and the USLE-LS (Slope Length and Steepness) factors, have been calculated using the TOPAGNPS algorithms [47] from the DEM generated using a combination of UAS-SfM photogrammetric techniques on the data collected in January 2016 [11]. The DEM has a 0.3 m horizontal spatial resolution, with a Root Mean Square Error (RMSE) of 0.07 m in the vertical and 0.03 m in the horizontal dimensions.
AnnAGNPS can utilize input parameters from the NRCS database developed for any location in the USA, including climate, soil, land use, and management properties [48]. For our field site in Mexico, fieldwork and laboratory analyses were necessary to acquire the needed information to apply AnnAGNPS in an ungauged watershed. Geologic maps may be relatively common, but the utility of such maps and their relationship to soil types must be determined with site-specific data. Soil candidates from the SSURGO database [38] were tested to choose the most suitable soil data, and were validated with field and laboratory measurements [48]. The Las Flores soil type was the most suitable SSURGO soil type for soils in SB, which are characterized by gentle to strong sloping on marine terraces, being moderately well-drained, having medium to rapid runoff, and very slow permeability. This description matches field observations in SB, and the corresponding soil samples are representatives of highly erodible soils according to Hanson’s soil classification diagram [49]. Percentage of impervious cover (IC) was calculated for the study area in SB from a vegetation-impervious-soil (VIS) map by Biggs et al. [50], as updated in Taniguchi et al. [51] to support the SCS curve number (CN) values used in this analysis. A land use map was generated by visual interpretation using the GoogleEarthTM imagery (11 November 2012, 2017 DigitalGlobe) into three land use categories: unpaved road (20%), housing (75%), and vacant lots (5%). IC was then calculated for each category and used to determine the default CN values. A composite curve number was calculated as the sum of the product of the fractional area coverage of each land cover category (unpaved road, housing, and vacant lots) multiplied by the CN associated with that category [52]. The same value (82 for soil type B) was used for all three cover categories, because in standard tables [52], residential areas have a lower CN than unpaved roads, but in our study area, vegetation was relatively sparse, and there is high connectivity between the roofs and the unpaved roads. Lacking additional data on runoff production from different surfaces, we left the CN for housing equal to the CN of unpaved roads and assume that the increased runoff from roofs is balanced by increased infiltration in vegetated areas on the lots. The USLE soil erodibility factor (K, 0.006 t·h·MJ−1·mm−1), and the saturated hydraulic conductivity (0.77 mm·h−1), were taken from the NRCS database [38] for the Las Flores soil type, and were assumed to have a uniform spatial distribution over the study area.
Eight soil samples were collected in the study area to estimate the critical shear stress (τc) and soil erodibility using a mini-jet erosion test following Hanson [49]. The submerged jet-test measures depth-of-scour, manually using a point gauge at known increments over time. τc is determined by the logarithmic-hyperbolic method described in Hanson and Cook [53]. Gordon et al. [33] noted that measured values of τc would be more accurate than any calculated values due to the large range and temporal and spatial variation of τc in the landscape. In our model, we use the measurements of τc and soil (head-cut) erodibility to set a default value, and use uncertainty analysis to determine if the final parameter range includes the measured values, as described below. Head-cut erodibility can be predicted as a power function of τc with coefficients a and b. The results from the jet-test suggested no consistent exponent value, so we assumed b = 0, and that erodibility was a constant value a, with the default value determined from the jet-test results.

2.5. Sensitivity Analysis

Sensitivity analysis was performed to explore and quantify the effect of input parameter variability on the output results. Many sensitivity analyses have been conducted on the AnnAGNPS model [26,27,28,31,54,55,56]. This paper focuses on the input parameters used to evaluate the capability of AnnAGNPS to simulate sediment production from gully erosion on unpaved roads. The sensitivity approach included varying the basic input variables that impact gully erosion modelling, with emphasis on runoff and soil erodibility, using Latin Hypercube Sampling (LHS) [57,58] to analyse the effect on gully erosion modelling, and to explore parameter sets that successfully simulate observed gully erosion. LHS was selected over other techniques such as orthogonal grid sampling because it is more efficient in terms of computational resources requirements. Orthogonal sampling requires more computational resources to perform the same analysis (6 parameters with 5 bins = ~15,000 models). LHS also ensures that each sample is collected in a fully stratified manner [59].
The most important parameters for gully erosion modelling [21] selected for LHS were (1) τc, (2) potential maximum soil moisture retention Smax = (1000/CN) − 10, (iii) TD, (iv) saturated hydraulic conductivity, (v) head-cut erodibility coefficient a, and (vi) Manning’s n for overland flow. A feasible parameter range was specified for each parameter (Table 1) based on measured (Jet-erosion test) and literature values [38,60,61]. For runoff generation, ranges were applied to Smax (Smax = (1000/CN) − 10) instead of the CN because the CN assigned to unpaved roads in SB is close to the upper limit value of CN (100) complicating the evaluation of higher values in the LHS.
LHS subdivides the range of each input parameter into N intervals of equal probability [57,58] then one value from each bin is chosen at random for each parameter to fit the desired sampling range. We used 15 bins to generate an initial 15 parameter ensembles. Preliminary tests suggested that these 15 parameter sets were insufficient to generate ensembles with the full range of parameter combinations, so 500 ensembles were generated by randomly selecting one of the 15 LHS-derived parameter values for each of the six parameters (Table 1).
The sensitivity of the simulated sediment load to variation in each input parameter was quantified using correlation analysis [58]. The linear correlation coefficient (LCC) measures the strength of the linear association between two parameters [57]. Partial correlation coefficient (PCC) measures the relationship between two parameters, with the effects of all other parameters constant. The PCC values were calculated using the algorithms within the pcor library of the R statistical software package [57]. Percent bias PBIAS was also used to estimate whether the average tendency of the simulated gully erosion rates was higher or lower than the observed data [56].
P B I A S = ∑ i = 1 n ( o b s e r v e d − s i m u l a t e d ) × 100 ∑ i = 1 n o b s e r v e d
According to Gupta et al. [62], a positive PBIAS value indicates a model underestimation, and a negative PBIAS indicates model overestimation.

2.6. Model Equifinality and Scenario Analysis

The 500 parameter ensembles were used to assess parameter identifiability and model equifinality for gully erosion modelling on unpaved roads using the AnnAGNPS model. A threshold of goodness of fit between observed and simulated gully erosion rates (SSL) was used to identify parameter ensembles that could be considered acceptable for a behavioural model. A threshold value of RMSE less than 1.2 mm (41% of the mean) was used as the threshold for behavioural models in this study, based on the comparison between observed and simulated SSL in nine sub-watersheds. An RMSE larger than 1.2 mm (41% of the mean) resulted in models with large errors for individual sub-watersheds, and were not used to test model equifinality. The threshold selected to divide behavioural from non-behavioural models is always subjective [63], and is based on the objectives of the analysis. Here, we aimed to quantify the impact of parameter uncertainty on scenario analyses, so we selected a threshold that yielded a tractable number of models for analysis (~20).
In order to test for trade-offs and compensation in parameter values, the correlation among parameters for behavioural models was quantified using Pearson’s correlation coefficients.
Behavioural models were used to quantitatively evaluate the runoff and sediment production from gully erosion on unpaved roads under two scenarios: (1) current conditions, and (2) paving all roads. Runoff production under the paved condition was simulated by increasing the CN values to reflect the runoff producing potential of impervious surfaces (CNscenario), which was calculated as:
C N s c e n a r i o = C N c c ( 1 − f A r o a d s ) + ( f A r o a d s × C N p a v e d )
where CNcc is the Curve Number under current conditions; fAroads is the fractional area of roads (0.2), and CNpaved is the Curve Number for paved roads (98).
We assumed that the CN for paved roads would be uniform, with relatively little uncertainty, so we did not perform a sensitivity analysis for the CN of paved roads. Gully sediment production was turned off under paved conditions since gullies formed exclusively on unpaved roads in the SB area. We assume that the drainage channel network is not modified in the paved scenario, since the change in elevations will be relatively minor. Road paving results in micro-topographic changes, such as routing flow from the centre of the street to side channels, but those alterations should not affect drainage areas or flowpaths at the sub-watershed scale.

3. Results

3.1. Sensitivity Analysis

Smax, TD, and τc were the most sensitive parameters in the gully erosion modelling (Table 2). TD correlated positively with gully erosion (Table 2), because higher scour depths erode more sediment during the upstream migration of the head-cut. Conversely, increasing Smax (decreasing the CN) and increasing τc reduced gully erosion, since increasing Smax reduces runoff, and increasing τc increases the resistance of the soil to detachment and erosion. Head-cut erodibility coefficient, Manning’s n, and saturated hydraulic conductivity did not have statistically significant correlations with sediment production from gully erosion (Table 2).

3.2. Behavioural Models and Parameter Identification

A total of 21 behavioural models were identified using the RMSE < 1.2 mm criterion. Simulated values of gully sediment production correlated with the observed values at the event scale, which illustrates the model’s ability to simulate gully erosion on unpaved roads over the study area (Figure 4). The RMSE of the simulated gully erosion rates using the default model was acceptable (2.1 mm, 70% of the mean), but a significant improvement was observed for the behavioural models (Figure 4). The AnnAGNPS behavioural models had relatively low errors (mean percent bias (PBIAS) ranging from −14.2 to 22.7). Model efficiencies were classified by Moriasi et al. [64] and Parajuli et al. [65] as being very good for ±16 ≤ PBIAS ≤ ±30 for SSL [56].
The behavioural models generally underestimated the largest observed sediment production from gully erosion (SSL > 5 mm) and tended to overestimate sediment production from sub-watersheds with less gully erosion (SSL < 4 mm) (Figure 4). Gully erosion contributes between 80% and 90% (87% on average) to the total sediment production among the behavioural models under current conditions.
The parameter ranges of the behavioural models were smaller than the initial ranges (Table 1), suggesting that the LHS method improves parameter identifiability in our watersheds. For example, τc in the behavioural models was 0.05–1.79 N·m-2, compared with the original range of 0.04 to 4 N·m−2. This corresponds to a soil texture of fine silt (0.05 N·m−2) to fine gravel (1.79 N·m−2) [60]. The parameter range for TD in the behavioural models was 0.63–0.95 m, compared with the original range of 0.3 to 2.4 m. Smax was the most sensitive parameter and was also relatively well constrained in the behavioural models between 35 nm and 57 mm (CN 82–88), compared with the original range of 28–84 mm (CN 75–90). Manning’s n, head-cut erodibility coefficient, and saturated hydraulic conductivity were not well constrained, but did not have a large impact on model output.
Some parameters were correlated in the behavioural models, suggesting that their values traded off or compensated for each other, resulting in parameter uncertainty (Table 3). Smax was inversely correlated with τc (p < 0.05), where lower values of τc were compensated by higher values of Smax in the behavioural models (Table 3). Higher values of Smax, which resulted in low runoff, required lower values of τc to maintain the same total sediment production.
The τc from the soil samples (N = 8) ranged from 0.15 to 1.9 N·m−2, and the erodibility ranged from 103 to 879 cm3 N−1·s−1 (Figure 5). τc from the samples spanned the range of the τc from the behavioural models, though there were some combinations of τc and erodibility that were slightly outside of the combinations observed, especially where τc was lower than observed for a given erodibility. Note that two behavioural models had the same critical shear stress value from the LHS, so only 20 open circles are visible in Figure 5.

3.3. Scenario Analysis: Equifinality

Total sediment loads and total and peak runoff for all of the behavioural models under the current conditions and roads paved scenarios are presented in Table 4.
Among the twenty-one behavioural models, sediment reduction from paving all roads varied from 90 to 94%, while total runoff was 1.46 to 2.27 times the unpaved condition, and peak runoff was 1.78 to 5.75 times the unpaved condition. The decrease in discharge under unpaved conditions results from higher potential maximum soil moisture retention for the modelled event. Other events that occur under higher antecedent moisture conditions may show a lower impact of paving.
A total of 3 out of the 21 behavioural models were identified as outliers in the equifinality analysis for scenario implications (Figure 6a,b). These 3 parameter ensembles, which showed the highest impacts on total and peak discharge, were characterized by high values of Smax, which results in lower runoff production under unpaved conditions and a larger percentage increase in overland flow under the road paved scenario. The sediment production ratio showed more robust results on the total sediment reduction (90 to 94%) for all the behavioural models (Figure 6c).

4. Discussion

Simulated sediment production from gullies was similar to the observed gully erosion in SB, suggesting that the AnnAGNPS model is able to estimate sediment production from unpaved roads in the study watersheds. The AnnAGNPS behavioural models had relatively similar errors (mean percent bias (PBIAS) ranged from −14.2 to 22.7, and 9 of 21 models have PBIAS less than 10) compared to previous AnnAGNPS applications simulating annual sediment loads (PBIAS = −7.1, Chahor et al. [56]), which used nine years of observed data for model calibration in a Mediterranean agricultural watershed.
Smax, TD and τc were the most sensitive parameters for the gully erosion model (Table 2). These results were consistent with previous studies that showed the importance of the runoff production in generating ephemeral gullies [21,26,31], and suggests that field measurements that can determine these parameters are useful for decreasing model uncertainty.
Runoff and soil resistance-to-erosion properties play an important role in gully erosion. The influence of these parameters on modelling the erosive process were reflected in the sensitivity analysis, where a trade-off was observed between the parameters related to runoff and soil erodibility (especially Smax and τc) in order to balance their respective influence in the gully erosion modelling, which is consistent with the significant correlation (p < 0.05) between Smax and τc (Table 3). Other parameters (saturated hydraulic conductivity, head-cut erodibility, Manning’s n) were correlated with other parameters in the behavioural models (p < 0.10), but these three parameters did not have significant impacts on sediment production from gullies, and further analysis on those correlations are beyond the scope of this paper.
Smax, TD, and τc were well-constrained in the behavioural models (Table 1). Manning’s n, head-cut erodibility and saturated hydraulic conductivity showed a wider range of values in the behavioural models, suggesting that these parameters are more influenced or compensated by other parameter combinations, complicating parameter identification.
High soil infiltration rates due to low antecedent soil moisture played a critical role in surface runoff generation in SB during the simulated storm event. CN Type II under “normal” soil moisture conditions was 82, while the CN for the modelled storm event, which was adjusted for antecedent soil moisture, was much lower (30), showing the impact of soil moisture on CN and runoff production. Low runoff production on dry soils had very important implications on the scenario analysis, resulting in a large increase in peak discharge (~1.8 to 5.7 times) under paved conditions. Other events that occur under conditions of higher antecedent moisture condition may show less impact of paving.
Field-measured values of τc helped to constrain the initial value for modelling. For example, the parameter range for τc in the behavioural models was 0.05–1.79, compared with the original range of 0.05 to 4. Approximately 80% of the behavioural models spanned the range of τc from of 0.05 to 1.1 N·m−2, which corresponds to a soil texture of fine silt (0.05) to very coarse sand (1.1) [60]. This suggests that field-measured τc is representative of the τc that controls the simulated sediment production in the study area.
Using different parameter ensembles generated by LHS allowed us to identify the range of the parameters and resulted in a better fit between the observed and the simulated gully erosion rates. Observed gully erosion rates were successfully reproduced using the 21 behavioural models (RMSE < 1.2 mm, < 41% of the mean), indicating robust simulated sediment production by gully erosion from unpaved roads.
The 21 behavioural models were consistent in terms of total sediment reduction (90–94%). Conversely, total runoff of the behavioural models increased from approximately 1.5 to 2.3 times and peak runoff increased by 1.8 to 5.7 times under the paved condition compared to the current condition. This could have significant impacts on the receiving earthen stream channels. The large increase in runoff generation under the paved roads scenario could be related to the large range in Smax in the behavioural models, and suggests that field data on infiltration rates, which could be used to generate values of Smax, is most critical for reducing model uncertainty. Soil compaction by car traffic on unpaved roads, which reduces infiltration rates, also has an impact on parameter uncertainty.
Increased runoff and changes in soil erosion rates due to road construction are well-known processes [13,15]. However, to our knowledge, this is the first attempt to simulate and evaluate model equifinality and implications for scenario analysis of ephemeral gully erosion rates in an urban environment. The AnnAGNPS model provides the capability of evaluating the impact of sediment management activities designed to mitigate gully erosion on unpaved roads. Road paving can be an effective sediment conservation practice, but the overall impact at the watershed scale—for example, the effect on receiving stream channels—needs to be assessed.

5. Conclusions

Gully formation and sediment yield were successfully simulated in an urban setting. Simulated Specific Soil Loss (SSL) using a model of gully erosion (AnnAGNPS) was similar to the observed SSL from gully erosion, with RMSE in SSL ranging from 0.96 to 1.2 mm for the twenty-one behavioural models, compared to 2.1 mm for the default parameter (Smax, TD, Manning’s n, saturated hydraulic conductivity, and head-cut erodibility) set. In the study area, gullies formed almost exclusively on unpaved roads, highlighting them as a major sediment source. Gully erosion may contribute significantly to the total sediment production, but other processes in the sediment budget need to be quantified for comparison. Smax (curve number), TD and τc were the most sensitive parameters in gully erosion modelling. The 21 behavioural models were consistent in their estimates of total sediment reduction when paving all roads (decrease of 90–94%). Conversely, total runoff of the behavioural models increased by approximately 1.5 to 2.3 times under the paved condition compared to under the current conditions. Our results suggest urgency in implementing management practices such as pavement or other stabilization measures of unpaved roads to mitigate soil erosion, but that paving may increase peak discharge significantly (by 1.8–5.7 times) at the neighbourhood scale. Our sensitivity analysis also identified the most uncertain parameters requiring further investigation to quantify the impacts of management on runoff and sediment production, especially parameters relating to infiltration capacity and runoff production. Future studies evaluating the effect of different soil types on gully erosion modelling using AnnAGNPS, as well as modelling the effect of other management actions (i.e., revegetation of hillslopes) on soil erosion and sediment loads, are crucial for proper management of sediment in our study area, and potentially in other urban areas in developing countries.

Acknowledgments

This study was funded by Consejo Nacional de Ciencia y Tecnología (CONACyT, Mexico) and the US Environmental Protection Agency (EPA) (Interagency Agreement ID # DW-12-92390601-0) in collaboration with the US Department of Agriculture (USDA, Agreement # 58-6408-4-015), San Diego State University (USA), University of Córdoba (Spain), and the Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE, Mexico). We thank Rajbir Parmar and Megan Mallard, for review and technical editing. Special thanks to the anonymous reviewers, who provided valuable feedback on the manuscript. The views expressed in this paper are those of the authors and do not necessarily represent the views or policies of the U.S. Environmental Protection Agency (EPA). Mention of trade names or products does not convey, and should not be interpreted as conveying, official EPA approval, endorsement or recommendation.

Author Contributions

All authors drafted the manuscript and participated in setting the objectives of this research. Trent W. Biggs made a significant contribution in optimizing LHS, and Ronald L. Bingner with model calibration. All authors contributed to the manuscript writing and revision.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Poesen, J.; Nachtergaele, J.; Verstraeten, G.; Valentin, C. Gully erosion and environmental change: Importance and research needs. Catena 2003, 50, 91–133. [Google Scholar] [CrossRef]
  2. Montgomery, D.R. Road surface drainage, channel initiation, and slope instability. Water Resour. Res. 1994, 30, 1925–1932. [Google Scholar] [CrossRef]
  3. Ziegler, A.D.; Giambelluca, T.W. Importance of rural roads as source areas for runoff in mountainous areas of northern Thailand. J. Hydrol. 1997, 196, 204–229. [Google Scholar] [CrossRef]
  4. Poesen, J.; Govers, G. Gully erosion in the loam belt of Belgium. In Soil Erosion on Agricultural Land, Proceedings of a Workshop Sponsored by the British Geomorphological Research Group, Coventry, Chicester, UK, January 1989; Boardman, J., Foster, I.D.L., Dearing, J., Eds.; Wiley: Chicester, UK, 1989; pp. 513–530. [Google Scholar]
  5. Momm, H.G.; Bingner, R.L.; Wells, R.R.; Wilcox, D. AGNPS GIS-based tool for watershed-scale identification and mapping of cropland potential ephemeral gullies. Appl. Eng. Agric. 2012, 28, 17–29. [Google Scholar] [CrossRef]
  6. Bingner, R.L.; Czajkowski, K.; Palmer, M.; Coss, J.; Davis, S.; Stafford, J.; Widman, N.; Theurer, F.; Koltum, G.; Richards, P.; et al. Upper Auglaize Watershed AGNPS Modelling Project: Final Report; Research Report No. 51; USDA-ARS National Sedimentation Laboratory: Oxford, MS, USA, 2006. [Google Scholar]
  7. Guerra, A.J.T.; Hoffman, H. Urban gully erosion in Brazil. Geography 2006, 19, 26–29. [Google Scholar]
  8. Adediji, A.; Jeje, L.K.; Ibitoye, M.O. Urban development and informal drainage patterns: Gully dynamics in Southwestern Nigeria. Appl. Geogr. 2013, 40, 90–102. [Google Scholar] [CrossRef]
  9. Trimble, S.W. Contribution of stream channel erosion to sediment yield from an urbanizing watershed. Science 1997, 278, 1442–1444. [Google Scholar] [CrossRef] [PubMed]
  10. Taniguchi, K.; Biggs, T.W. Regional impacts of urbanization on stream channel geometry: A case study in semi-arid southern California. Geomorphology 2015, 248, 228–236. [Google Scholar] [CrossRef]
  11. Gudino-Elizondo, N.; Biggs, T.; Castillo, C.; Bingner, R.; Langendoen, E.; Taniguchi, K.; Kretzschmar, T.; Yuan, Y.; Liden, D. Ephemeral gully erosion rates and topographical thresholds in an urban watershed using Unmanned Aerial Systems and structure from motion photogrammetric techniques. Land Degrad. Dev. 2018, in press. [Google Scholar]
  12. Wemple, B.C.; Browning, T.; Ziegler, A.D.; Celi, J.; Chun, K.P.; Jaramillo, F.; Leite, N.; Ramchunder, S.J.; Negishi, J.N.; Palomeque, X.; et al. Ecohydrological disturbances associated with roads: Current knowledge, research needs, and management concerns with reference to the Tropics. Ecohydrology 2017. [Google Scholar] [CrossRef]
  13. Reid, L.M.; Dunne, T. Sediment production from forest road surfaces. Water Resour. Res. 1984, 20, 1753–1761. [Google Scholar] [CrossRef]
  14. Montgomery, D.R.; Dietrich, W.E. Where do channels begin? Nature 1988, 336, 232–234. [Google Scholar] [CrossRef]
  15. Ramos-Scharrón, C.E.; MacDonald, L.H. Runoff and suspended sediment yields from an unpaved road segment, St. John, US Virgin Islands. Hydrol. Process. 2007, 21, 35–50. [Google Scholar] [CrossRef]
  16. Castillo, C.; Gómez, J.A. A century of gully erosion research: Urgency, complexity and study approaches. Earth-Sci. Rev. 2016, 160, 300–319. [Google Scholar] [CrossRef]
  17. Merkel, W.H.; Woodward, D.E.; Clarke, C.D. Ephemeral gully erosion model (EGEM). In Modeling Agricultural, Forest, and Rangeland Hydrology, Proceedings of the International Symposium, Chicago, IL, USA, 12–13 December 1988; American Society of Agricultural Engineers: St. Joseph, MI, USA, 1988; pp. 315–323. [Google Scholar]
  18. Woodward, D.E. Method to predict cropland ephemeral gully erosion. Catena 1999, 37, 393–399. [Google Scholar] [CrossRef]
  19. Bingner, R.L.; Theurer, F.D. Physics of suspended sediment transport in AnnAGNPS. In Proceedings of the 2nd Federal Interagency Hydrologic Modeling Conference, Las Vegas, NV, USA, 28 July–1 August 2002. [Google Scholar]
  20. Bingner, R.L.; Theurer, F.D.; Yuan, Y.P. AnnAGNPS Technical Processes. Washington, D.C. US Department of Agriculture (USDA)—Agricultural Research Service (ARS). 2015. Available online: https://www.wcc.nrcs.usda.gov/ftpref/wntsc/H&H/AGNPS/downloads/AnnAGNPS_Technical_Documentation.pdf (accessed on 5 July 2017).
  21. Taguas, E.V.; Yuan, Y.; Bingner, R.L.; Gomez, J.A. Modeling the contribution of ephemeral gully erosion under different soil managements: A case study in an olive orchard microcatchment using the AnnAGNPS model. Catena 2012, 98, 1–16. [Google Scholar] [CrossRef]
  22. Merritt, W.S.; Letcher, R.A.; Jakeman, A.J. A review of erosion and sediment transport models. Environ. Model. Softw. 2003, 18, 761–799. [Google Scholar] [CrossRef]
  23. Bull, L.J.; Kirkby, M.J. Gully processes and modelling. Progress Phys. Geogr. 1997, 21, 354–374. [Google Scholar] [CrossRef]
  24. Casali, J.; Gimenez, R.; Bennett, S. Gully erosion processes: Monitoring and modelling. Earth Surf. Process. Landf. 2009, 34, 1839–1840. [Google Scholar] [CrossRef]
  25. Nachtergaele, J.; Poesen, J.; Vandekerckhove, L.; Oostwould-Wijdenes, D.J.; Roxo, M. Testing the ephemeral gully erosion model (EGEM) for two Mediterranean environments. Earth Surf. Process. Landf. 2001, 26, 17–30. [Google Scholar] [CrossRef]
  26. Yuan, Y.; Bingner, R.L.; Rebich, R.A. Evaluation of AnnAGNPS on Mississippi Delta MSEA watersheds. Trans. ASAE 2001, 44, 1673–1682. [Google Scholar] [CrossRef]
  27. Yuan, Y.; Bingner, R.L.; Rebich, R.A. Evaluation of AnnAGNPS nitrogen loading in an agricultural watershed. J. Am. Water Resour. Assoc. 2003, 39, 457–466. [Google Scholar] [CrossRef]
  28. Yuan, Y.; Bingner, R.L.; Theurer, F.D.; Rebich, R.A.; Moore, P.A. Phosphorus component in AnnAGNPS. Trans. ASAE 2005, 48, 2145–2154. [Google Scholar] [CrossRef]
  29. Baginska, B.; Milne-Home, W.A. Parameter sensitivity in calibration and validation of an annualized agricultural non-point source model. Water Sci. Appl. 2003, 6, 331–345. [Google Scholar] [CrossRef]
  30. Suttles, J.B.; Vellidis, G.; Bosch, D.; Lowrance, R.; Sheridan, J.M.; Usery, E.L. Watershed-scale simulation of sediment and nutrient loads in Georgia Coastal Plain streams using the Annualized AGNPS model. Trans. ASAE 2003, 46, 1325–1335. [Google Scholar] [CrossRef]
  31. Licciardello, F.; Zema, D.A.; Zimbone, S.M.; Bingner, R.L. Runoff and soil erosion evaluation by the AnnAGNPS Model in a small Mediterranean watershed. Trans. ASAE 2007, 50, 1585–1593. [Google Scholar] [CrossRef]
  32. Shamshad, A.; Leow, C.S.; Ramlah, A.; Wan Hussin, W.M.A.; Mohd Sanusi, S.A. Applications of AnnAGNPS model for soil loss estimation and nutrient loading for Malaysian conditions. Int. J. Appl. Earth Obs. Geoinf. 2008, 10, 239–252. [Google Scholar] [CrossRef]
  33. Gordon, L.M.; Bennett, S.J.; Bingner, R.L.; Theurer, F.D.; Alonso, C.V. Simulating ephemeral gully erosion in AnnAGNPS. Trans. ASAE 2007, 50, 857–866. [Google Scholar] [CrossRef]
  34. Alonso, C.V.; Bennett, S.J.; Stein, O.R. Predicting headcut erosion and migration in concentrated flows typical of upland areas. Water Resour. Res. 2002, 38. [Google Scholar] [CrossRef]
  35. Beven, K.J.; Freer, J. Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology. J. Hydrol. 2001, 249, 11–29. [Google Scholar] [CrossRef]
  36. Hornberger, G.M.; Spear, R.C. An approach to the preliminary analysis of environmental systems. J. Environ. Manag. 1981, 12, 7–18. [Google Scholar]
  37. Beven, K.J. Prophecy, reality and uncertainty in distributed hydrological modelling. Adv. Water Resour. 1993, 16, 41–51. [Google Scholar] [CrossRef]
  38. Natural Resources Conservation Service Soils. 2018. Available online: https://www.nrcs.usda.gov/wps/portal/nrcs/main/soils/survey/ (accessed on 4 February 2018).
  39. Grover, R. Local Perspectives on Environmental Degradation and Community Infrastructure in Los Laureles Canyon, Tijuana. Master’s Thesis, San Diego State University, San Diego, CA, USA, 2011. [Google Scholar]
  40. CalEPA, California Environmental Protection Agency. 2018. Available online: https://www.waterboards.ca.gov/sandiego/water_issues/programs/tmdls/TijuanaRiverValley.shtml (accessed on 4 February 2018).
  41. Liden, T.D.; US Environmental Protection Agency, San Diego, CA, USA. Personal communication, 2016.
  42. Biggs, T.W.; Taniguchi, K.T.; Gudino-Elizondo, N.; Yongping, Y.; Bingner, R.L.; Langendoen, E.J.; Liden, D. Field Measurements to Support Sediment and Hydrological Modelling in Los Laureles Canyon; U.S. Environmental Protection Agency: Washington, DC, USA, 2018.
  43. SCS. Hydrology, National Engineering Handbook; Natural Resources Conservation Service, US Department of Agriculture: Washington, DC, USA, 1972; Chapter 4. Available online: https://www.nrcs.usda.gov/wps/portal/nrcs/detailfull/national/water/?&cid=stelprdb1043063 (accessed on 15 November 2017).
  44. SCS. Technical Release 55: Urban Hydrology for Small Watersheds; Natural Resources Conservation Service, US Department of Agriculture: Washington, DC, USA, 1986. Available online: https://www.nrcs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb1044171.pdf (accessed on 15 November 2017).
  45. Wells, R.R.; Momm, H.G.; Rigby, J.R.; Bennett, S.J.; Bingner, R.L.; Dabney, S.M. An empirical investigation of gully widening rates in upland concentrated flows. Catena 2013, 101, 114–121. [Google Scholar] [CrossRef]
  46. Bingner, R.L.; Wells, R.R.; Momm, H.G.; Rigby, J.R.; Theurer, F.D. Ephemeral gully channel width and erosion simulation technology. Nat. Hazards 2016, 80, 1949–1966. [Google Scholar] [CrossRef]
  47. Garbrecht, J.; Martz, L.W. TOPAGNPS, An Automated Digital Landscape Analysis Tool for Topographic Evaluation, Drainage Identification Watershed Segmentation and Subcatchment Parameterization for AGNPS. Watershed Modelling Technology; Agricultural Research Service: Beltsville, MD, USA, 1999. Available online: https://www.ars.usda.gov/ARSUserFiles/60600505/AGNPS/Dataprep/TOPAGNPS_userman.pdf (accessed on 1 June 2017).
  48. Biggs, T.W.; Taniguchi, K.T.; Gudino-Elizondo, N.; Yongping, Y.; Bingner, R.L.; Langendoen, E.J.; Liden, D. Geology, soil properties and erosion on marine terraces along the US-Mexico Border. Manuscript in preparation. 2018. [Google Scholar]
  49. Hanson, G.J. Surface erodibility of earthen channels at high stresses part II—Developing an in situ testing device. Trans. ASAE 1990, 33, 132–137. [Google Scholar] [CrossRef]
  50. Biggs, T.W.; Atkinson, E.; Powell, R.; Ojeda-Revah, L. Land cover following rapid urbanization on the US-Mexico border: Implications for conceptual models of urban watershed processes. Lands. Urban Plan. 2010, 96, 78–87. [Google Scholar] [CrossRef]
  51. Taniguchi, K.T.; Biggs, T.W.; Langendoen, E.J.; Castillo, C.; Gudino-Elizondo, N.; Yongping, Y.; Liden, D. Stream channel erosion in a rapidly urbanizing region of the US–Mexico border: Documenting the importance of channel hardpoints with Structurefrom-Motion photogrammetry. Earth Surf. Process. Landf. 2018. [Google Scholar] [CrossRef]
  52. Dunne, T.; Leopold, L.B. Water in Environmental Planning; W.H. Freeman and Company: New York, NY, USA, 1987; p. 818. ISBN 0-71670079-4. [Google Scholar]
  53. Hanson, G.J.; Cook, K.R. Procedure to estimate soil erodibility for water management purposes. In Proceeding of the Mini-Conference Advance in Water Quality Modeling, Toronto, ON, Canada, 18–21 July 1999; ASAE 1999, Paper No. 992133; ASAE: St. Joseph, MI, USA, 1999. [Google Scholar]
  54. Das, S.; Rudra, R.P.; Gharabaghi, B.; Gebremeskel, V.; Goel, P.K.; Dickinson, W.T. Applicability of AnnAGNPS for Ontario conditions. Can. Biosyst. Eng. 2008, 50, 1.1–1.11. [Google Scholar]
  55. Zema, D.A.; Bingner, R.L.; Denisi, P.; Govers, G.; Licciardello, F.; Zimbone, S.M. Evaluation of runoff, peak flow and sediment yield for events simulated by the AnnAGNPS model in a Belgian agricultural watershed. Land Degrad. Dev. 2012, 23, 205–215. [Google Scholar] [CrossRef]
  56. Chahor, Y.; Casalí, J.; Giménez, R.; Bingner, R.L.; Campo, M.A.; Goñi, M. Evaluation of the AnnAGNPS model for predicting runoff and sediment yield in a small Mediterranean agricultural watershed in Navarre (Spain). Agric. Water Manag. 2014, 134, 24–37. [Google Scholar] [CrossRef]
  57. Kim, S.H.; Yi, S. Understanding relationship between sequence and functional evolution in yeast proteins. Genetica 2007, 131, 151–156. [Google Scholar] [CrossRef] [PubMed]
  58. Momm, H.G.; Bingner, R.L.; Yuan, Y.; Locke, M.A.; Wells, R.R. Spatial Characterization of Riparian Buffer Effects on Sediment Loads from Watershed Systems. J. Environ. Qual. 2014, 43, 1736–1753. [Google Scholar] [CrossRef] [PubMed]
  59. McKay, M.D.; Beckman, R.J.; Conover, W.J. A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Am. Stat. Assoc. 1979, 21, 239–245. [Google Scholar] [CrossRef]
  60. U.S. GEOLOGICAL SURVEY, Scientific Investigations Report 2008–5093, Table 7. 2018. Available online: https://pubs.usgs.gov/sir/2008/5093/table7.html (accessed on 5 January 2018).
  61. Engman, E.T. Roughness coefficients for routing surface runoff. J. Irrig. Drain. Eng. 1986, 112, 39–53. [Google Scholar] [CrossRef]
  62. Gupta, H.V.; Sorooshian, S.; Yapo, P.O. Status of automatic calibration for hydrologic models: Comparison with multi-level expert calibration. J. Hydrol. Eng. 1999, 4, 135–143. [Google Scholar] [CrossRef]
  63. Beven, K.J.; Binley, A. The future of distributed models: Model calibration and uncertainty prediction. Hydrol. Process. 1992, 6, 279–298. [Google Scholar] [CrossRef]
  64. Moriasi, D.; Arnold, J.; Van Liew, M.; Bingner, R.L.; Harmel, R.; Veith, T. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 2007, 50, 885–900. [Google Scholar] [CrossRef]
  65. Parajuli, P.B.; Nelson, N.O.; Frees, L.D.; Mankin, K.R. Comparison of AnnAGNPS and SWAT model simulation results in USDA-CEAP agricultural watersheds in south-central Kansas. Hydrol. Process. 2009, 23, 748–763. [Google Scholar] [CrossRef]
Geosciences 08 00137 g001 550
Figure 1. (a) UAS-SfM-derived orthophoto for San Bernardo (SB), and the 9 study watersheds with their outlets; (b) Geographic location of the Los Laureles Canyon Watershed (LLCW), SB, and the Tijuana River Estuarine Reserve (TJE); (c) one example of land degradation caused by gully erosion in Tijuana, Mexico.
Figure 1. (a) UAS-SfM-derived orthophoto for San Bernardo (SB), and the 9 study watersheds with their outlets; (b) Geographic location of the Los Laureles Canyon Watershed (LLCW), SB, and the Tijuana River Estuarine Reserve (TJE); (c) one example of land degradation caused by gully erosion in Tijuana, Mexico.
Geosciences 08 00137 g001
Geosciences 08 00137 g002 550
Figure 2. Daily rainfall time series for the 2016 water year. The grey box represents the rainfall threshold (~25–35 mm) for gully formation observed in the study area.
Figure 2. Daily rainfall time series for the 2016 water year. The grey box represents the rainfall threshold (~25–35 mm) for gully formation observed in the study area.
Geosciences 08 00137 g002
Geosciences 08 00137 g003 550
Figure 3. (a) Digitized gullies, watershed boundary, outlet, and locations of field measurements of gully depths; (b) An example of field measurement of gully depth.
Figure 3. (a) Digitized gullies, watershed boundary, outlet, and locations of field measurements of gully depths; (b) An example of field measurement of gully depth.
Geosciences 08 00137 g003
Geosciences 08 00137 g004 550
Figure 4. Relationship between observed and simulated Specific Soil Loss (SSL, the average depth of soil loss in the watershed in mm) from gully erosion in San Bernardo, Tijuana, Mexico, obtained from 21 behavioural models. The blue dots show the results from the default model parameters (Table 1).
Figure 4. Relationship between observed and simulated Specific Soil Loss (SSL, the average depth of soil loss in the watershed in mm) from gully erosion in San Bernardo, Tijuana, Mexico, obtained from 21 behavioural models. The blue dots show the results from the default model parameters (Table 1).
Geosciences 08 00137 g004
Geosciences 08 00137 g005 550
Figure 5. τc and head cut erodibility as measured by the jet-test (black dots) compared with other values from the literature (lines), and with the parameters from the behavioural models (open circles).
Figure 5. τc and head cut erodibility as measured by the jet-test (black dots) compared with other values from the literature (lines), and with the parameters from the behavioural models (open circles).
Geosciences 08 00137 g005
Geosciences 08 00137 g006 550
Figure 6. Impacts on water and sediment load ratios between current conditions and paving-all-roads scenario using the 21 behavioural models.
Figure 6. Impacts on water and sediment load ratios between current conditions and paving-all-roads scenario using the 21 behavioural models.
Geosciences 08 00137 g006
Table
Table 1. Parameter default values, parameter range, and the actual parameter ranges obtained using LHS and for the parameter ensembles that gave acceptable errors (behavioural models).
Table 1. Parameter default values, parameter range, and the actual parameter ranges obtained using LHS and for the parameter ensembles that gave acceptable errors (behavioural models).
ParameterDefault ValuesParameter RangeLHS-Derived Parameter Range, All Models (N = 500)Behavioural Models Parameter Range, (N = 21)
MinMaxMinMaxMinMax
Smax55.75 mm27.8783.6327.9380.8435.1856.85
Saturated conductivity50 mm·d−155005.514385.51438
Critical shear stress1 N·m−20.0440.053.250.051.79
Manning’s n0.150.0150.30.0170.290.0170.22
Tillage depth0.60 m0.32.40.332.310.630.95
Head-cut erodibility1000 g·N−1·s−1150175021317132131562
Table
Table 2. Sensitivity analysis of the effect of variability in potential maximum soil moisture retention, tillage depth, critical shear stress, head-cut erodibility, Manning’s n, and saturated hydraulic conductivity on sediment production by gully erosion using the Linear (LCC) and Partial (PCC) correlations.
Table 2. Sensitivity analysis of the effect of variability in potential maximum soil moisture retention, tillage depth, critical shear stress, head-cut erodibility, Manning’s n, and saturated hydraulic conductivity on sediment production by gully erosion using the Linear (LCC) and Partial (PCC) correlations.
VariableLCCPCC
Smax−0.58 *−0.77 *
Tillage depth0.44 *0.72 *
Critical shear stress−0.48 *−0.71 *
Headcut erodibility−0.10−0.03
Manning’s n0.010.05
Saturated conductivity0.020.01
* p < 0.05.
Table
Table 3. Correlation coefficients for input parameters of the behavioural models.
Table 3. Correlation coefficients for input parameters of the behavioural models.
ParameterSmaxHead Cut ErodibilitySaturated ConductivityCritical Shear StressManning’s nTillage Depth
Smax10.030.05−0.51 *−0.18−0.31
Head cut erodibility 1−0.42 †0.14−0.270.24
Saturated conductivity 10.110.110.10
Critical shear stress 1−0.210.43 †
Manning’s n 1−0.44 †
Tillage depth 1
* indicates p < 0.05; Numbers with the symbol (†) indicate p < 0.10.
Table
Table 4. Modelled peak discharge (L/s), total discharge volume (Q, m3), and sediment load (tons) at the outlet under unpaved and paved conditions for 21 behavioural models.
Table 4. Modelled peak discharge (L/s), total discharge volume (Q, m3), and sediment load (tons) at the outlet under unpaved and paved conditions for 21 behavioural models.
Peak (L/s)Q (m3)Sediment (tons)
Unpaved
min4148513
mean50500787
max1017391048
Paved
min2033749
mean10579959
max181107867
Ratio of Paved: Unpaved
min1.781.460.06
mean2.731.700.08
max5.752.270.10

Share and Cite

MDPI and ACS Style

Gudino-Elizondo, N.; Biggs, T.W.; Bingner, R.L.; Yuan, Y.; Langendoen, E.J.; Taniguchi, K.T.; Kretzschmar, T.; Taguas, E.V.; Liden, D. Modelling Ephemeral Gully Erosion from Unpaved Urban Roads: Equifinality and Implications for Scenario Analysis. Geosciences 2018, 8, 137. https://doi.org/10.3390/geosciences8040137

AMA Style

Gudino-Elizondo N, Biggs TW, Bingner RL, Yuan Y, Langendoen EJ, Taniguchi KT, Kretzschmar T, Taguas EV, Liden D. Modelling Ephemeral Gully Erosion from Unpaved Urban Roads: Equifinality and Implications for Scenario Analysis. Geosciences. 2018; 8(4):137. https://doi.org/10.3390/geosciences8040137

Chicago/Turabian Style

Gudino-Elizondo, Napoleon, Trent W. Biggs, Ronald L. Bingner, Yongping Yuan, Eddy J. Langendoen, Kristine T. Taniguchi, Thomas Kretzschmar, Encarnacion V. Taguas, and Douglas Liden. 2018. "Modelling Ephemeral Gully Erosion from Unpaved Urban Roads: Equifinality and Implications for Scenario Analysis" Geosciences 8, no. 4: 137. https://doi.org/10.3390/geosciences8040137

APA Style

Gudino-Elizondo, N., Biggs, T. W., Bingner, R. L., Yuan, Y., Langendoen, E. J., Taniguchi, K. T., Kretzschmar, T., Taguas, E. V., & Liden, D. (2018). Modelling Ephemeral Gully Erosion from Unpaved Urban Roads: Equifinality and Implications for Scenario Analysis. Geosciences, 8(4), 137. https://doi.org/10.3390/geosciences8040137

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop