Journal cover Journal topic
Geoscientific Model Development An interactive open-access journal of the European Geosciences Union
Journal topic
GMD | Articles | Volume 12, issue 7
Geosci. Model Dev., 12, 2837-2854, 2019
https://doi.org/10.5194/gmd-12-2837-2019
© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.
Geosci. Model Dev., 12, 2837-2854, 2019
https://doi.org/10.5194/gmd-12-2837-2019
© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.

Model description paper 11 Jul 2019

Model description paper | 11 Jul 2019

r.sim.terrain 1.0: a landscape evolution model with dynamic hydrology

r.sim.terrain 1.0: a landscape evolution model with dynamic hydrology
Brendan Alexander Harmon1, Helena Mitasova2,3, Anna Petrasova2,3, and Vaclav Petras2,3 Brendan Alexander Harmon et al.
  • 1Robert Reich School of Landscape Architecture, Louisiana State University, Baton Rouge, LA, USA
  • 2Center for Geospatial Analytics, North Carolina State University, Raleigh, NC, USA
  • 3Department of Marine, Earth, and Atmospheric Sciences, North Carolina State University, Raleigh, NC, USA
Abstract
Back to toptop

While there are numerical landscape evolution models that simulate how steady-state flows of water and sediment reshape topography over long periods of time, r.sim.terrain is the first to simulate short-term topographic change for both steady-state and dynamic flow regimes across a range of spatial scales. This free and open-source Geographic Information Systems (GIS)-based topographic evolution model uses empirical models for soil erosion and a physics-based model for shallow overland water flow and soil erosion to compute short-term topographic change. This model uses either a steady-state or unsteady representation of overland flow to simulate how overland sediment mass flows reshape topography for a range of hydrologic soil erosion regimes based on topographic, land cover, soil, and rainfall parameters. As demonstrated by a case study for the Patterson Branch subwatershed on the Fort Bragg military installation in North Carolina, r.sim.terrain simulates the development of fine-scale morphological features including ephemeral gullies, rills, and hillslopes. Applications include land management, erosion control, landscape planning, and landscape restoration.

1 Introduction
Back to toptop

Landscape evolution models represent how the surface of the Earth changes over time in response to physical processes. Most studies of landscape evolution have been descriptive, but a number of numerical landscape evolution models have been developed that simulate elevational change over time (Tucker and Hancock2010; Temme et al.2013). Numerical landscape evolution models such as the Geomorphic – Orogenic Landscape Evolution Model (GOLEM) (Tucker and Slingerland1994), CASCADE (Braun and Sambridge1997), the Channel-Hillslope Integrated Landscape Development (CHILD) model (Tucker et al.2001), CAESAR (Coulthard et al.2002, 2012), SIBERIA (Willgoose2005), LAPSUS (Schoorl et al.2000, 2002) r.landscape.evol (Barton et al.2010), and eSCAPE (Salles2019) simulate landscape evolution driven primarily by steady-state flows over long temporal scales. Landlab (2019) (http://landlab.github.io/, last access: 3 July 2019), a new Python library for numerically modeling Earth surface processes (Hobley et al.2017), has components for simulating landscape evolution such as the Stream Power with Alluvium Conservation and Entrainment (SPACE) model (Shobe et al.2017). While Geographic Information Systems (GIS) support efficient data management, spatial and statistical modeling and analysis, and visualization, there are few GIS-based soil erosion models (see Table 1) or landscape evolution models. Thaxton (2004) developed the model r.terradyn as a Geographic Resources Analysis Support System (GRASS) GIS shell script module to simulate terrain evolution by steady-state net erosion–deposition rates estimated by the Simulation of Water Erosion (SIMWE) model (Mitas and Mitasova1998) and gravitational diffusion. Barton et al. (2010) developed a long-term landscape evolution model in GRASS GIS called r.landscape.evol that integrates the Unit Stream Power Erosion Deposition (USPED) model, fluvial erosion, and gravitational diffusion. r.landscape.evol has been used to simulate the impact of prehistoric settlements on Mediterranean landscapes. In spite of the recent progress in landscape evolution modeling and monitoring, there are still major research questions to address in the theoretical foundations of erosion modeling such as how erosional processes scale over time and space, and how sediment detachment and transport interact (Mitasova et al.2013). While most numerical landscape evolution models simulate erosion processes at steady-state peak flows, short-term erosional processes like gully formation can be driven by unsteady, dynamic flow with significant morphological changes happening before flows reach steady state. A landscape evolution model with dynamic water and sediment flow is needed to study fine-scale spatial and short-term erosional processes such as gully formation and the development of microtopography.

Mitasova et al. (1996)Mitasova et al. (1996)Mitas and Mitasova (1998)Flanagan et al. (2013)Guertin et al. (2015)Ad de Roo et al. (1996)Hobley et al. (2017)

Table 1Examples of geospatial soil erosion models.

Download Print Version | Download XLSX

At the beginning of a rainfall event, overland water flow is unsteady – its depth changes at a variable rate over time and space. If the intensity of rainfall continues to change throughout the event, then the flow regime will remain dynamic. If, however, overland flow reaches a peak rate, then the hydrologic regime is considered to be at steady state. At steady state,

(1)h(x,y,t)t=0,

where (x,y) is the position [m], t is the time [s], and h(x,y,t) is the depth of overland flow [m].

Gullies are eroded, steep-banked channels formed by ephemeral, concentrated flows of water. A gully forms when overland water flow converges in a knickzone – a concave space with steeper slopes than its surroundings (Zahra et al.2017) – during intense rainfall events. When the force of the water flow concentrated in the knickzone is enough to detach and transport large amounts of sediment, an incision begins to form at the apex of the knickzone – the knickpoint or headwall. As erosion continues, the knickpoint begins to migrate upslope and the nascent gully channel widens, forming steep channel banks. Multiple incisions initiated by different knickpoints may merge into a gully channel and multiple channels may merge into a branching gully system (Mitasova et al.2013). This erosive process is dynamic; the morphological changes drive further changes in a positive feedback loop. When the gully initially forms, the soil erosion regime should be detachment capacity limited with the concentrated flow of water in the channel of the gully detaching large amounts of sediment and transporting it to the foot of the gully, potentially forming a depositional fan. If the intensity of rainfall decreases and transport and detachment capacity approach a balance, then the soil erosion regime may switch to a variable erosion–deposition regime, in which soil is eroded and deposited in a spatially variable pattern. Subsequent rainfall events may trigger further knickpoint formation and upslope migration, channel incision and widening, and depositional fan and ridge formation. Between high-intensity rainfall events, lower-intensity events and gravitational diffusion may gradually smooth the shape of the gully. Eventually, if detachment capacity significantly exceeds transport capacity and the regime switches to transport capacity limited, the gully may fill with sediment as soil continues to be eroded but cannot be transported far.

Gully erosion rates and evolution can be monitored in the field or modeled on the computer. Field methods include dendrogeomorphology (Malik2008) and permanent monitoring stakes for recording erosion rates, extensometers for recording mass wasting events, weirs for recording water and suspended sediment discharge rates, and time series of surveys using total station theodolites (Thomas et al.2004), unmanned aerial systems (UASs) (Jeziorska et al.2016; Kasprak et al.2019; Yang et al.2019), airborne lidar (Perroy et al.2010; Starek et al.2011), and terrestrial lidar (Starek et al.2011; Bechet et al.2016; Goodwin et al.2016; Telling et al.2017). With terrestrial lidar, airborne lidar, and UAS photogrammetry, there are now sufficient-resolution topographic data to morphometrically analyze and numerically model fine-scale landscape evolution in GIS, including processes such as gully formation and the development of microtopography. Gully erosion has been simulated with RUSLE2-Raster (RUSLER) in conjunction with the Ephemeral Gully Erosion Estimator (EphGEE) (Dabney et al.2014), while gully evolution has been simulated for detachment-capacity-limited erosion regimes with the Simulation of Water Erosion (SIMWE) model (Koco2011; Mitasova et al.2013). Now numerical landscape evolution models that can simulate steady-state and unsteady flow regimes and can dynamically switch between soil erosion regimes are needed to study fine-scale spatial and short-term erosional processes.

The numerical landscape evolution model r.sim.terrain was developed to simulate the spatiotemporal evolution of landforms caused by shallow overland water and sediment flows at spatial scales ranging from square meters to kilometers and temporal scales ranging from minutes to years. This open-source GIS-based landscape evolution model can simulate either steady-state or unsteady flow regimes, dynamically switch between soil erosion regimes, and simulate the evolution of fine-scale morphological features such as ephemeral gullies (Fig. 1). It was designed as a research tool for studying how erosional processes scale over time and space, comparing empirical and process-based models, comparing steady-state and unsteady flow regimes, and studying the role of unsteady flow regimes in fine-scale morphological change. r.sim.terrain was tested with a subwatershed scale (450 m2) case study and the simulations were compared against a time series of airborne lidar surveys.

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f01

Figure 1The digital elevation model (DEM) (a) before and (b) after simulated landscape evolution with r.sim.terrain for a subwatershed of Patterson Branch, Fort Bragg, NC, USA. The before DEM was generated from an airborne lidar data acquired in 2012. The simulation used the SIMWE model for a 120 min rainfall event with 50 mm h−1 for a variable erosion–deposition regime at steady state. In the evolved DEM, the gully channel has widened with depositional ridges forming along its thalweg.

Download

2 r.sim.terrain
Back to toptop

The process-based, spatially distributed landscape evolution model r.sim.terrain simulates topographic changes caused by shallow, overland water flow across a range of spatiotemporal scales and soil erosion regimes using either the Simulated Water Erosion (SIMWE) model, the 3-Dimensional Revised Universal Soil Loss Equation (RUSLE3D) model, or the USPED model (Fig. 2). The r.sim.terrain model can simulate either steady-state or dynamic flow regimes. SIMWE is a physics-based simulation that uses a Monte Carlo path sampling method to solve the water and sediment flow equations for detachment-limited, transport-limited, and variable erosion–deposition soil erosion regimes (Mitas and Mitasova1998; Mitasova et al.2004). With SIMWE, r.sim.terrain uses the modeled flow of sediment – a function of water flow and soil detachment and transport parameters – to estimate net erosion and deposition rates. RUSLE3D is an empirical equation for estimating soil erosion rates in detachment-capacity-limited soil erosion regimes (Mitasova et al.1996, 2013). With RUSLE3D, r.sim.terrain uses an event-based rainfall erosivity factor, soil erodibility factor, land cover factor, and 3-D topographic factor – a function of slope and flow accumulation – to model soil erosion rates. USPED is a semi-empirical equation for net erosion and deposition in transport-capacity-limited soil erosion regimes (Mitasova et al.1996, 2013). With USPED, r.sim.terrain uses an event-based rainfall erosivity factor, soil erodibility factor, land cover factor, and a topographic sediment transport factor to model net erosion or deposition rates as the divergence of sediment flow. For each of the models, topographic change is derived at each time step from the net erosion–deposition rate and gravitational diffusion. Depending on the input parameters, r.sim.terrain simulations with SIMWE can represent variable soil erosion–deposition regimes, including prevailing detachment-capacity-limited or prevailing transport-capacity-limited regimes.

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f02

Figure 2Conceptual diagram for r.sim.terrain.

Download

The r.sim.terrain model can simulate the evolution of gullies including processes such as knickpoint migration, channel incision, channel widening, aggradation, scour pit formation, depositional ridge formation along the thalweg of the gully, and depositional fan formation at the foot of the gully. Applications include geomorphological research, erosion control, landscape restoration, and scenario development for landscape planning and management. This model can simulate landscape evolution over a wide range of spatial scales from small watersheds less than 10 km2 with SIMWE to regional watersheds of 100 km2 with USPED or RULSE3D, although it does not model fluvial processes. It has been used at resolutions ranging from submeter scale to 30 m. The model has been implemented as a Python add-on module for the free, open-source GRASS GIS (https://grass.osgeo.org/, last access: 3 July 2019) (GRASS Development Team2019). The source code is available at https://github.com/baharmon/landscape_evolution (last access: 3 July 2019) under the GNU General Public License v2 (Harmon2019a). It supports multithreading and parallel processing to efficiently compute simulations using large, high-resolution topographic datasets. The landscape evolution model can be installed in GRASS GIS as an add-on module with the command (Harmon2019f):

g.extension extension=r.sim.terrain

2.1 Landscape evolution

Landscape evolution in r.sim.terrain is driven by change in the elevation surface caused by soil erosion and deposition. During storm events, overland flow erodes soil and transports sediment across landscape, and under favorable conditions deposits the sediment. Gravitational diffusion, applied to the changed elevation surface, simulates the smoothing effects of localized soil transport between events.

2.1.1 Elevation change

Assuming negligible uplift, the change in elevation over time is described by the continuity of mass equation expressed as the divergence of sediment flow (Tucker et al.2001):

(2)zt=-qsρs-1=dsρs-1,

where z is elevation [m], t is time [s], qs is sediment flow per unit width (vector) [kg m−1 s−1], ds is the net erosion–deposition rate [kg m−2 s−1], and ρs is sediment mass density [kg m−3].

In r.sim.terrain, the net erosion–deposition rate ds driven by overland flow is estimated at different levels of complexity based on the simulation mode selected by the user. Gravitational diffusion is then applied to the changed topography to simulate the smoothing effects of localized soil transport between rainfall events. The change in elevation due to gravitational diffusion is a function of the diffusion coefficient and the Laplacian of elevation (Thaxton2004):

(3)zt=εg2z,

where εg is the diffusion coefficient [m2 s−1].

The discrete implementation follows Thaxton (2004):

(4)zt+Δt1=zt+Δzs(5)zt+Δt1+Δt2=zt+Δt1+Δzg,

where Δzs is elevation change [m] caused by net erosion or deposition during time interval Δt1 (Eq. 2), and Δzg is the diffusion-driven elevation change [m] during time interval Δt1 (Eq. 3).

2.1.2 Erosion–deposition regimes

Following experimental observations and qualitative arguments, Foster et al. (1977) proposed that the sum of the ratio of the net erosion–deposition rate ds to the detachment capacity Dc [kg m−2 s−1] and the ratio of the sediment flow rate qs=|qs| to the sediment transport capacity Tc [kg m−1 s−1] is a conserved quantity (unity):

(6)dsDc+qsTc=1.

The net erosion and deposition rate ds can then be expressed as being proportional to the difference between the sediment transport capacity Tc and the actual sediment flow rate qs:

(7)ds=DcTc(Tc-qs).

This principle is used in several erosion models including the Water Erosion Prediction Project (WEPP) (Flanagan et al.2013) and SIMWE (Mitas and Mitasova1998).

Using this concept, it is possible to identify two limiting erosion–deposition regimes. When TcDc leading to Tcqs, the erosion regime is detachment capacity limited and net erosion is equal to the detachment capacity:

(8)ds=Dc.

For this case, the transport capacity of overland flow exceeds the detachment capacity, and thus sediment flow, erosion, and sediment transport are limited by the detachment capacity. Therefore, no deposition occurs. An example of this case is when a strong storm producing intense overland flow over compacted clay soils causes high-capacity flows to transport light clay particles, while the detachment of compacted soils is limited. When DcTc, sediment flow is at sediment transport capacity qs=Tc, leading to a transport-capacity-limited regime with deposition reaching its maximum extent for the given water flow. Net erosion–deposition is computed as the divergence of transport capacity multiplied by a unit vector s0 in the direction of flow:

(9)ds=Tcs0.

This case may occur, for example, during a moderate storm with overland flow over sandy soils with high detachment capacity but low transport capacity. For 0<(Dc/Tc)<, the spatial pattern of net erosion–deposition is variable and depends on the difference between the sediment transport capacity and the actual sediment flow rate at the given location.

The detachment capacity Dc and the sediment transport capacity Tc are estimated using shear stress and stream power equations, respectively, expressed as power functions of water flow properties and slope angle. The relations between the topographic parameters of well-known empirical equations for erosion modeling, such as the Universal Soil Loss Equation (USLE) and stream power, were presented by Moore and Burch (1986) and used to develop simple, GIS-based models for limiting erosion–deposition cases such as RUSLE3D and USPED (Mitasova and Mitas2001). The SIMWE model estimates Tc and Dc using modified equations and parameters developed for the WEPP model (Flanagan et al.2013; Mitasova et al.2013).

The simulation modes in r.sim.terrain include (Fig. 2)

  • the process-based SIMWE model for steady-state and unsteady shallow overland flow in variable erosion–deposition regimes with ds computed by solving the shallow water flow and sediment transport continuity equations,

  • the RUSLE3D model for detachment-capacity-limited cases with ds given by Eq. (8), and

  • the USPED model for transport-capacity-limited regimes with ds given by Eq. (9).

The following sections explain the computation of ds for these three modes in more detail.

2.2 Simulation of Water Erosion (SIMWE)

SIMWE is a physics-based simulation of shallow overland water and sediment flow that uses a path sampling method to solve the continuity equations with a 2-D diffusive wave approximation (Mitas and Mitasova1998; Mitasova et al.2004). SIMWE has been implemented in GRASS GIS as the modules r.sim.water and r.sim.sediment. In SIMWE mode, for each landscape evolution time step, r.sim.terrain

  • computes the first-order partial derivatives of the elevation surface z/x and z/y,

  • simulates shallow water flow depth, sediment flow, and the net erosion–deposition rate, and

  • then evolves the topography based on the erosion–deposition rate and gravitational diffusion.

The first-order partial derivatives of the elevation surface are computed using the GRASS GIS module r.slope.aspect using the equations in Hofierka et al. (2009). r.sim.terrain simulates unsteady-state flow regimes when the landscape evolution time step is less than the travel time for a drop of water or a particle of sediment to cross the landscape, e.g., when the time step is less than the time to concentration for the modeled watershed. With longer landscape evolution time steps, the model simulates a steady-state regime.

2.2.1 Shallow water flow

The SIMWE model simulates shallow overland water flow controlled by spatially variable topographic, soil, land cover, and rainfall parameters using a Green function Monte Carlo path sampling method. The steady-state shallow water flow continuity equation relates the change in water depth across space to source, defined in our case as rainfall excess rate:

(10)q=(hv)=n-1h5/3s1/2s0=ie,

where q is the water flow per unit width (vector) [m2 s−1], h is the depth of overland flow [m], v is the water flow velocity vector [m s−1] whose magnitude is computed with Manning's equation v=n-1h2/3s1/2, n is Manning's coefficient [s m-1/3], s is slope (unitless), and ie is the rainfall excess rate [m s−1] (i.e., rainfall intensity infiltration vegetation intercept).

An approximation of diffusive wave effects is incorporated by adding a diffusion term proportional to 2[h5∕3]:

(11)-εw22h5/3+n-1h5/3s1/2s0=ie,

where εw is a spatially variable diffusion coefficient [m4∕3 s−1].

The path sampling method solves the continuity equation for h5∕3 through the accumulation of the evolving source (Mitasova et al.2004). The solution assumes that water flow velocity is largely controlled by the slope of the terrain and surface roughness and that its change at a given location during the simulated event is negligible. The initial number of particles per grid cell is proportional to the rainfall excess rate ie (source). The water depth h5∕3 at time τ during the simulated rainfall event is computed as a function of particle (walkers) density at each grid cell. Particles are routed across the landscape by finding a new position for each walker at time ττ:

(12)rmnew=rm+Δτv+g,

where r=(x,y) is the mth walker position [m], Δτ is the particle routing time step [s], and g is a random vector with Gaussian components with variance Δτ [m].

The mathematical background of the method, including the computation of the temporal evolution of water depth and incorporation of approximate momentum through an increased diffusion rate in the prevailing direction of flow, is presented by Mitas and Mitasova (1998) and Mitasova et al. (2004).

2.2.2 Sediment flow and net erosion–deposition

The SIMWE model simulates the sediment flow over complex topography with spatially variable overland flow, soil, and land cover properties by solving the sediment flow continuity equation using a Green function Monte Carlo path sampling method. Steady-state sediment flow qs is approximated by the bivariate continuity equation, which relates the change in sediment flow rate to effective sources and sinks:

(13)qs=sources-sinks=ds.

The sediment flow rate qs is a function of water flow and sediment concentration (Mitas and Mitasova1998):

(14)qs=ρscq=ρschv=ϱv,

where ρs is sediment mass density in the water column [kg m−3], c is sediment concentration [particle m−3], and ϱ=ρsch is the mass of sediment transported by water per unit area [kg m−2].

The sediment flow equation (Eq. 13), like the water flow equation, has been rewritten to include a small diffusion term that is proportional to the mass of water-carried sediment per unit area 2ϱ (Mitas and Mitasova1998):

(15)-εs22ϱ+(ϱv)+ϱDcTc|v|=ds,

where εs is the diffusion constant [m2 s−1].

On the left-hand side of Eq. (15), the first term describes local diffusion, the second term is drift driven by water flow, and the third term represents a velocity-dependent “potential” acting on the mass of transported sediment. The initial number of particles per grid cell is proportional to the soil detachment capacity Dc (source). The particles are then routed across the landscape by finding a new position for each walker at time ττ:

(16)rmnew=rm+Δτv+g,

while the updated weight is

(17)wmnew=wmexp-Δτurmnew+urm/2,

where u=Dc/Tc|v|.

Sediment flow is computed as the product of weighted particle densities and the water flow velocity (Eq. 14), and the net erosion–deposition rate ds is computed as the divergence of sediment flow using Eq. (13). See Mitas and Mitasova (1998) and Mitasova et al. (2004) for more details on the Green function Monte Carlo solution and equations for computing Dc and Tc.

This model can simulate erosion regimes from prevailing detachment-limited conditions when TcDc to prevailing transport-capacity-limited conditions when DcTc and the erosion–deposition patterns between these conditions. At each landscape evolution time step, the regime can change based on the ratio between the sediment detachment capacity Dc and the sediment transport capacity Tc and the actual sediment flow rate. If the landscape evolution time step is shorter than the time to concentration (i.e., the time for water to reach steady state), then net erosion–deposition is derived from unsteady flow.

2.3 Revised Universal Soil Loss Equation for Complex Terrain (RUSLE3D)

RUSLE3D is an empirical model for computing erosion in a detachment-capacity-limited soil erosion regime for watersheds with complex topography (Mitasova et al.1996). It is based on the USLE, an empirical equation for estimating the average sheet and rill soil erosion from rainfall and runoff on agricultural fields and rangelands with simple topography (Wischmeier et al.1978). It models erosion-dominated regimes without deposition in which sediment transport capacity is uniformly greater than detachment capacity. In USLE, soil loss per unit area is determined by an erosivity factor R, a soil erodibility factor K, a slope length factor L, a slope steepness factor S, a cover management factor C, and a prevention measures factor P. These factors are empirical constants derived from an extensive collection of measurements on 22.13 m standard plots with an average slope of 9 %. RUSLE3D was designed to account for more complex, 3-D topography with converging and diverging flows. In RUSLE3D, the topographic potential for erosion at any given point is represented by a 3-D topographic factor LS3-D, which is a function of the upslope contributing area and the angle of the slope.

In this spatially and temporally distributed model, RUSLE3D is modified by the use of a event-based R factor derived from rainfall intensity at each time step. For each time step, this model computes the parameters for RUSLE3D – an event-based erosivity factor, the slope of the topography, the flow accumulation, and the 3-D topographic factor – and then solves the RUSLE3D equation for the rate of soil loss (i.e., the net soil erosion rate). The soil erosion rate is then used to simulate landscape evolution in a detachment-capacity-limited soil erosion regime.

2.3.1 Erosivity factor

The erosivity factor R in USLE and RUSLE is the combination of the total energy and peak intensity of a rainfall event, representing the interaction between the detachment of sediment particles and the transport capacity of the flow. It can be calculated as the product of the the kinetic energy of the rainfall event E and its maximum 30 min intensity I30 (Brown and Foster1987; Renard et al.1997; Panagos et al.2015, 2017). In this model, however, the erosivity factor is derived at each time step as a function of kinetic energy, rainfall depth, rainfall intensity, and time. First, rain energy is derived from rainfall intensity (Brown and Foster1987; Yin et al.2017):

(18)ere0=1.-bexpiri0,

where er is unit rain energy [MJ ha−1 mm−1], ir is rainfall intensity [mm h−1], b is empirical coefficient, i0 is reference rainfall intensity [mm h−1], and e0 is reference energy [MJ ha−1 mm1]. The parameters for this equation were derived from observed data published for different regions by Panagos et al. (2017). Then the event-based erosivity index Re is calculated as the product of unit rain energy, rainfall depth, rainfall intensity, and time:

(19)Re=ervrirΔt,

where Re is the event-based erosivity index [MJ mm ha−1 h−1], vr is the rainfall depth [mm] derived from vr=irΔt, and Δt is the change in time [s].

2.3.2 Flow accumulation

The upslope contributing area per unit width a is determined by flow accumulation (the number of grid cells draining into a given grid cell) multiplied by grid cell width (Fig. 3d). Flow accumulation is calculated using a multiple flow direction algorithm (Metz et al.2009) based on AT least-cost path searches (Ehlschlaeger1989). The multiple flow direction algorithm implemented in GRASS GIS, as the module r.watershed is computationally efficient, does not require sink filling and can navigate nested depressions and other obstacles.

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f03

Figure 3Water and sediment flows modeled with spatially variable land cover for Patterson Branch, Fort Bragg, NC: (a) water depth [m] simulated by SIMWE for a 10 min event with 50 mm h−1 in the subwatershed; (b) flow accumulation for RUSLE3D in the subwatershed; (c) erosion and deposition [kg m−2 s−1] simulated by SIMWE in drainage area 1; and (d) erosion [kg m−2 s−1] modeled by RUSLE3D in drainage area 1.

Download

2.3.3 3-D topographic factor

The 3-D topographic factor LS3-D is calculated as a function of the upslope contributing area and the slope (Fig. 3e):

(20)LS3-D=(m+1)aa0msinββ0n,

where LS3-D is the dimensionless topographic factor, a is upslope contributing area per unit width [m], a0 is the length of the standard USLE plot [22.1 m], β is the angle of the slope [], m is an empirical coefficient, n is an empirical coefficient, and β0 is the slope of the standard USLE plot [5.14].

The empirical coefficients m and n for the upslope contributing area and the slope can range from 0.2 to 0.6 and 1.0 to 1.3, respectively, with low values representing dominant sheet flow and high values representing dominant rill flow.

2.3.4 Detachment-limited erosion rate

The erosion rate is a function of the event-based erosivity factor, soil erodibility factor, 3-D topographic factor, land cover factor, and prevention measures factor (Fig. 3d):

(21)E=ReKLS3-DCP,

where E is soil erosion rate (soil loss) [kg m−2 min−1], Re is the event-based erosivity factor [MJ mm ha−1 h−1], K is the soil erodibility factor [t ha h ha−1 MJ−1 mm−1], LS3-D is the dimensionless topographic (length–slope) factor, C is the dimensionless land cover factor, and P is the dimensionless prevention measures factor.

The detachment-limited erosion represented by RUSLE3D leads to the simulated change in elevation:

(22)Δzs=Dcρs-1=Eρs-1,

which is combined with Eq. (3) for gravitational diffusion.

2.4 Unit Stream Power Erosion Deposition (USPED)

USPED estimates net erosion–deposition as the divergence of sediment flow in a transport-capacity-limited soil erosion regime. The amount of soil detached is close to the amount of sediment that water flow can carry. As a transport-capacity-limited model, USPED predicts erosion where transport capacity increases and deposition where transport capacity decreases. The influence of topography on sediment flow is represented by a topographic sediment transport factor, while the influence of soil and land cover is represented by factors adopted from USLE and RUSLE (Mitasova et al.1996). Sediment flow is estimated by computing the event-based erosivity factor (Re) using Eq. (19), the slope and aspect of the topography, the flow accumulation with a multiple flow direction algorithm, the topographic sediment transport factor, and sediment flow at transport capacity. Net erosion–deposition is then computed as the divergence of sediment flow.

2.4.1 Topographic sediment transport factor

Using the unit stream power concept presented by Moore and Burch (1986), the 3-D topographic factor (Eq. 20) for RUSLE3D is modified to represent the topographic sediment transport factor (LST) – the topographic component of overland flow at sediment transport capacity:

(23)LST=amsinβn,

where LST is the topographic sediment transport factor, a is the upslope contributing area per unit width [m], β is the angle of the slope [], m is an empirical coefficient, and n is an empirical coefficient.

2.4.2 Transport-limited sediment flow and net erosion–deposition

Sediment flow at transport capacity is a function of the event-based rainfall factor, soil erodibility factor, topographic sediment transport factor, land cover factor, and prevention measures factor:

(24)T=ReKCPLST,

where T is sediment flow at transport capacity [kg m−1 s−1], Re is the event-based rainfall factor [MJ mm ha−1 h−1], K is the soil erodibility factor [t ha h ha−1 MJ−1 mm−1], C is the dimensionless land cover factor, and P is the dimensionless prevention measures factor.

Net erosion–deposition is estimated as the divergence of sediment flow, assuming that sediment flow is equal to sediment transport capacity:

(25)ds=Tccosαx+Tcsinαy,

where ds is net erosion–deposition [kg m−2 s−1], α is the aspect of the topography (i.e., the direction of flow) []. With USPED, the simulated change in elevation Δzs=ds is derived from Eq. (2) for landscape evolution and then Eq. (3) for gravitational diffusion.

3 Case study
Back to toptop

Military activity is a high-impact land use that can cause significant physical alteration to the landscape. Erosion is a major concern for military installations, particularly at training bases, where the land surface is disturbed by off-road vehicles, foot traffic, and munitions. Off-road vehicles and foot traffic by soldiers cause the loss of vegetative cover, the disruption of soil structure, soil compaction, and increased runoff due to reduced soil capacity for water infiltration (Webb and Wilshire1983; McDonald2004). Gullies – ephemeral channels with steep headwalls that incise into unconsolidated soil to depths of meters – are a manifestation of erosion common to military training installations like Fort Bragg in North Carolina and the Piñon Canyon Maneuver Site in Colorado. While the local development of gullies can restrict the maneuverability of troops and vehicles during training exercises, pervasive gullying across a landscape can degrade an entire training area (Huang and Niemann2014). To test the effectiveness of the different models in r.sim.terrain, we compared the simulated evolution of a highly eroded subwatershed of Patterson Branch on Fort Bragg, North Carolina, against a time series of airborne lidar surveys. The models – SIMWE, RUSLE3D, and USPED – were tested in steady-state and dynamic modes for design storms with constant rainfall.

3.1 Patterson Branch

With 650 km2 of land, Fort Bragg is the largest military installation in the US and has extensive areas of bare, erodible soils on impact areas, firing ranges, landing zones, and drop zones. It is located in the Sandhills region of North Carolina with a longleaf pine and wiregrass ecosystem (Sorrie et al.2006). The study landscape – a subwatershed of Patterson Branch (Fig. 4) in the Coleman Impact Area – is pitted with impact craters from artillery and mortar shells and has an active, approximately 2 m deep gully. It is a pine-scrub oak Sandhills community composed primarily of longleaf pine (Pinus palustris) and wiregrass (Aristida stricta) on Blaney and Gilead loamy sands (Sorrie2004). Throughout the Coleman Impact Area, frequent fires ignited by live munitions drive the ecological disturbance regime of this fire-adapted ecosystem. In 2016, the 450 m2 study site was 43.24 % bare ground with predominately loamy sands, 39.54 % covered by the wiregrass community, and 17.22 % forested with the longleaf pine community (Fig. 5a). We hypothesize that the elimination of forest cover in the impact zone triggered extensive channelized overland flow, gully formation, and sediment transport into the creek.

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f04

Figure 4Subwatershed with 2014 orthoimagery (a) draped over the 2016 digital elevation model and (b) drainage areas with 2014 orthoimagery, Patterson Branch, Fort Bragg, NC, USA.

Download

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f05

Figure 5Morphological change in the subwatershed of Patterson Branch, Fort Bragg, NC, USA: (a) land cover in 2014, (b) landforms in 2012, (c) elevation difference between 2012 and 2016 [m], and (d) landforms in 2016.

Download

Time series of digital elevation models and land cover maps for the study landscape were generated from lidar point clouds and orthophotography. The digital elevation models for 2004, 2012, and 2016 were interpolated at 1 m resolution using the regularized spline with tension function (Mitasova and Mitas1993; Mitasova et al.2005) from airborne lidar surveys collected by the NC Floodplain Mapping Program and Fort Bragg. Unsupervised image classification was used to identify clusters of spectral reflectance in a time series of 1 m resolution orthoimagery collected by the National Agriculture Imagery Program. The land cover maps were derived from the classified lidar point clouds and the classified orthoimagery. Spatially variable soil erosion factors – the k factor, c factor, Manning's coefficient, and runoff rate – were then derived from the land cover and soil maps. The dataset for this study is hosted at https://github.com/baharmon/landscape_evolution_dataset (last access: 3 July 2019) under the ODC Open Database License (ODbL) (Harmon2019b). The data are derived from publicly available data from the US Army, USGS, USDA, Wake County GIS, NC Floodplain Mapping Program, and the NC State Climate Office. There are detailed instructions for preparing the input data in the tutorial (https://github.com/baharmon/landscape_evolution/blob/master/tutorial.md, last access: 3 July 2019, Harmon2019c) and a complete record of the commands used to process the sample data in the data log (https://github.com/baharmon/landscape_evolution_dataset/blob/master/nc_spm_evolution/DATA.md, last access: 3 July 2019, Harmon2019d).

We used the geomorphons method of automated landform classification based on the openness of terrain (Jasiewicz and Stepinski2013) and the difference between the digital elevation models to analyze the changing morphology of the study area (Figs. 5 and 6). The 2 m deep gully – its channels classified as valleys and its scour pits as depressions by geomorphons – has multiple mature branches and ends with a depositional fan. The gully has also developed depositional ridges beside the channels. Deep scour pits have developed where branches join the main channel and where the main channel has sharp bends. A new branch has begun to form in a knickzone classified as a mix of valleys and hollows on a grassy swale on the northeast side of the gully. Between 2012 and 2016 a depositional ridge developed at the foot of this nascent branch where it would meet the main channel. The 2016 minus 2012 DEM of difference (DoD) – i.e., the difference in elevation (Figs. 5c and 6c) – shows a deepening of the main channel by approximately 0.2 m and scours pits by approximately 1 m, while depositional ridges have formed and grown up to approximately 1 m high. The DoD also shows that 244.60 m3 of sediment were deposited on the depositional fan between 2012 and 2016.

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f06

Figure 6Detailed morphological change for drainage area 1 of Patterson Branch, Fort Bragg, NC, USA: (a) land cover in 2014, (b) landforms in 2012, (c) elevation difference between 2012 and 2016 [m], and (d) landforms in 2016.

Download

3.2 Simulations

We ran a sequence of r.sim.terrain simulations with design storms for the Patterson Branch subwatershed study area to demonstrate the capabilities of the RUSLE3D, USPED, and SIMWE models (Table 2). To analyze the results of the simulations, we compared net differences in elevation morphological features, and volumetric change. While r.sim.terrain can use rainfall records, we used design storms to demonstrate and test the basic capabilities of the model. Our design storms were based off the peak rainfall values in records from the State Climate Office of North Carolina. We used RUSLE3D to simulate landscape evolution in a dynamic, detachment-capacity-limited soil erosion regime for a 120 min design storm with 3 min intervals and a constant rainfall intensity of 50 mm h−1 (Fig. 7). We used USPED to simulate landscape evolution in a dynamic, transport-capacity-limited soil erosion regime for a 120 min design storm with 3 min intervals and a constant rainfall intensity of 50 mm h−1 (Fig. 8). We used SIMWE to simulate landscape evolution in a steady state, variable erosion–deposition soil erosion regime for a 120 min design storm with a constant rainfall intensity of 50 mm h−1 (Fig. 9). In all of the simulations, a sink-filling algorithm – an optional parameter in r.sim.terrain – was used to reduce the effects of positive feedback loops that cause the overdevelopment of scour pits.

Table 2Landscape evolution simulations.

Download Print Version | Download XLSX

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f07

Figure 7Dynamic simulation with RUSLE3D for a 120 min event with a rainfall intensity of 50 mm h−1 for Patterson Branch, Fort Bragg, NC: (a) flow accumulation and (b) erosion [kg m−2 s−1] for the subwatershed in the final 3 min time step; (c) net difference [m] and (d) landforms for drainage area 1.

Download

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f08

Figure 8Dynamic simulation with USPED for a 120 min event with a rainfall intensity of 50 mm h−1 for Patterson Branch, Fort Bragg, NC: (a) flow accumulation and (b) erosion–deposition [kg m−2 s−1] for the subwatershed in the final 3 min time step; (c) net difference [m] and (d) landforms for drainage area 1.

Download

https://www.geosci-model-dev.net/12/2837/2019/gmd-12-2837-2019-f09

Figure 9Steady-state SIMWE simulations for a 120 min event with a rainfall intensity of 50 mm h−1 for Patterson Branch, Fort Bragg, NC: (a) depth of unsteady flow [m] and (b) erosion–deposition [kg m−2 s−1] for the subwatershed; (c) net difference [m] and (d) landforms for drainage area 1.

Download

The simulations were automated and run in parallel using Python scripts that are available in the software repository (https://github.com/baharmon/landscape_evolution, last access: 3 July 2019, Harmon2019a). The simulations can be reproduced using these scripts and the study area dataset by following the instructions in the Open Science Framework repository at https://osf.io/tf6yb/ (last access: 3 July 2019). The simulations were run in GRASS GIS 7.4 on a desktop computer with 64 bit Ubuntu 16.04.4 LTS, 8×4.20 GHz Intel Core i7 7700K CPUs, and 32 GB RAM. Simulations using SIMWE are far more computationally intensive than RULSE3D or USPED but support multithreading when compiled with OpenMP. Dynamic simulations of RUSLE3D and USPED took 2 min 36 s and 3 min 14 s, respectively, to run on a single thread, while the steady-state simulation for SIMWE took 44 min 51 s to run on six threads (Table 2).

3.3 Results

We used the difference in DEMs to compute volumetric changes between the lidar surveys and the simulations (Table 3). We applied a threshold of ±0.18 m to the lidar surveys since they had a vertical accuracy at a 95 % confidence level of 18.15 cm based on a 9.25 cm root mean square error in z (RMSEz) for non-vegetated areas in accordance with the National Digital Elevation Program guidelines (North Carolina Risk Management Office2018). Given the presence of the mature gully with ridges along its banks, we hypothesize that the study landscape had previously been dominated by a detachment-limited soil erosion regime but – given the net change of 654.77 m3 – had switched to a transport-capacity-limited or variable erosion–deposition regime during our study period.

Table 3Volumetric change.

Download Print Version | Download XLSX

The dynamic RUSLE3D simulation carved a deep incision in the main gully channel where water accumulated (Fig. 7). As a detachment-capacity-limited model, RUSLE3D's results were dominated by erosion and thus negative elevation change. It eroded 1480.75 m3 of sediment with no deposition.

The dynamic USPED simulation eroded the banks of the gully and deposited in channels causing the gully grow wider and shallower (Fig. 8). As a transport-capacity-limited model, USPED generated a distributed pattern with both erosion and deposition. Erosion far exceeded deposition with 1235.08 m3 of sediment eroded and 727.46 m3 deposited for a net change of −507.62 m3. While USPED's pattern of elevation change was grainy and fragmented, it captured the process of channel filling and widening expected with a transport-capacity-limited soil erosion regime.

The steady-state SIMWE simulation for a variable erosion–deposition regime predicted the morphological processes and features expected of its regime including gradual aggradation, channel widening, the formation of depositional ridges along the thalweg of the channel, and the development of the depositional fan (Fig. 9). SIMWE was the closest to the observed baseline volumetric change. It balanced erosion and deposition with 785.56 m3 of sediment eroded and 608.91 m3 deposited for a net change of −149.66 m3. Only the SIMWE simulation deposited sediment on the depositional fan. While the difference of lidar surveys showed that 244.60 m3 of sediment were deposited on the fan, SIMWE predicted that 54.05 m3 would be deposited.

SIMWE was unique in simulating unsteady flows (Fig. 9a) and fine-scale geomorphological processes such as the development of depositional ridges and a depositional fan. While USPED generated a grainy pattern of erosion and deposition, it was much faster than SIMWE (Table 2) and still simulated the key morphological patterns and processes – channel incision, filling, and widening. Given their speed and approximate modeling of erosive processes, RUSLE3D and USPED are effective for simulating landscape evolution on large rasters. RUSLE3D, for example, has been used to model erosion for the entire 650 km2 Fort Bragg installation at 9 m resolution (Levine et al.2018).

4 Discussion
Back to toptop

Limitations of this landscape evolution model include shallow overland flow, units, computation time, and raster size. r.sim.terrain only models shallow overland flows, not fluvial processes or subsurface flows. It requires data – including elevation and rainfall intensity – in metric units. The implementation of SIMWE in GRASS GIS is computationally intensive and may require long computation times even with multithreading. Because SIMWE uses a Green function Monte Carlo solution of the sediment transport equation, the accuracy, detail, and smoothness of the results depend on the number of random walkers. While a large number of random walkers will reduce the numerical error in the path sampling solution, it will also greatly increase computation time. A customized compilation of GRASS GIS is needed to run SIMWE with more than 7 million random walkers. This limits the size of rasters that can be easily processed with SIMWE, while RUSLE3D and USPED are much faster, computationally efficient, and can easily be run on much larger rasters.

In the future, we plan to assess this model by comparing simulations against a monthly time series of submeter-resolution surveys by unmanned aerial systems and terrestrial lidar. We also plan to develop a case study demonstrating how the model can be used as a planning tool for landscape restoration. Planned enhancements to the model include modeling subsurface flows, accounting for bedrock, and a reverse landscape evolution mode for backward modeling.

5 Conclusions
Back to toptop

The short-term landscape evolution model r.sim.terrain can simulate the development of gullies, rills, and hillslopes by overland water erosion for a range of hydrologic and soil erosion regimes. The model is novel for simulating landscape evolution based on unsteady flows. The landscape evolution model was tested with a series of simulations for different hydrologic and soil erosion regimes for a highly eroded subwatershed on Fort Bragg with an active gully. For each regime, it generated the morphological processes and features expected. The physics-based SIMWE model simulated morphological processes for a variable erosion–deposition regime such as gradual aggradation, channel widening, scouring, the development of depositional ridges along the thalweg, and the growth of the depositional fan. The empirical RUSLE3D model simulated channel incision in a detachment-limited soil erosion regime, while the semi-empirical USPED model simulated channel widening and filling in a transport-limited regime. Since r.sim.terrain is a GIS-based model that simulates fine-scale morphological processes and features, it can easily and effectively be used in conjunction with other GIS-based tools for geomorphological research, land management and conservation, erosion control, and landscape restoration.

Code and data availability
Back to toptop
Code and data availability. 

As a work of open science, this study is reproducible, repeatable, and recomputable. Since the data, model, GIS, and dependencies are all free and open source, the study can easily be reproduced. The landscape evolution model has been implemented in Python as a module for GRASS GIS, a free and open-source GIS. The source code for the model is hosted on GitHub at https://github.com/baharmon/landscape_evolution (last access: 3 July 2019) under the GNU General Public License version 2 (Harmon2019a). The code repository also includes Python scripts for running and reproducing the simulations in this paper. The digital object identifier (DOI) for the version of the software documented in this paper is https://doi.org/10.5281/zenodo.3243699 (Harmon2019a). There are detailed instructions for running this model in the manual at https://grass.osgeo.org/grass76/manuals/addons/r.sim.terrain.html (last access: 3 July 2019) (Harmon2019f) and the tutorial at https://github.com/baharmon/landscape_evolution/blob/master/tutorial.md (last access: 3 July 2019) (Harmon2019c). The geospatial dataset for the study area is available on GitHub at https://github.com/baharmon/landscape_evolution_dataset (last access: 3 July 2019) (Harmon2019b) under the Open Database License (https://opendatacommons.org/licenses/odbl/, last access: 3 July 2019) with the DOI: https://doi.org/10.5281/zenodo.3243700 (Harmon2019b). The data log (https://github.com/baharmon/landscape_evolution_dataset/blob/master/nc_spm_evolution/DATA.md, last access: 3 July 2019) has a complete record of the commands used to process the sample data. The source code, scripts, data, and results are also hosted on the Open Science Framework at https://osf.io/tf6yb/ (last access: 3 July 2019) (Harmon2019e) with the DOI https://doi.org/10.17605/osf.io/tf6yb (Harmon2019e).

Author contributions
Back to toptop
Author contributions. 

BH developed the models, code, data, case studies, and manuscript. HM contributed to the development of the models and case studies and revised the manuscript. AP and VP contributed to the development of the code. All authors read and approved the final manuscript.

Competing interests
Back to toptop
Competing interests. 

The authors declare that they have no conflict of interest.

Acknowledgements
Back to toptop
Acknowledgements. 

We acknowledge the GRASS GIS Development Community for developing and maintaining GRASS GIS.

Review statement
Back to toptop
Review statement. 

This paper was edited by Bethanna Jackson and reviewed by three anonymous referees.

References
Back to toptop

Barton, C. M., Ullah, I., and Mitasova, H.: Computational Modeling and Neolithic Socioecological Dynamics: a Case Study from Southwest Asia, Am. Antiquity, 75, 364–386, available at: http://www.jstor.org/stable/25766199 (last access: 3 July 2019), 2010. a, b

Bechet, J., Duc, J., Loye, A., Jaboyedoff, M., Mathys, N., Malet, J.-P., Klotz, S., Le Bouteiller, C., Rudaz, B., and Travelletti, J.: Detection of seasonal cycles of erosion processes in a black marl gully from a time series of high-resolution digital elevation models (DEMs), Earth Surf. Dynam., 4, 781–798, https://doi.org/10.5194/esurf-4-781-2016, 2016. a

Braun, J. and Sambridge, M.: Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretization, Basin Res., 9, 27–52, https://doi.org/10.1046/j.1365-2117.1997.00030.x, 1997. a

Brown, L. C. and Foster, G. R.: Storm Erosivity Using Idealized Intensity Distributions, Transactions of the American Society of Agricultural Engineers, 30, 0379–0386, https://doi.org/10.13031/2013.31957, 1987. a, b

Coulthard, T. J., Macklin, M. G., and Kirkby, M. J.: A cellular model of Holocene upland river basin and alluvial fan evolution, Earth Surf. Proc. Land., 27, 269–288, https://doi.org/10.1002/esp.318, 2002. a

Coulthard, T. J., Hancock, G. R., and Lowry, J. B. C.: Modelling soil erosion with a downscaled landscape evolution model, Earth Surf. Proc. Land., 37, 1046–1055, https://doi.org/10.1002/esp.3226, 2012. a

Dabney, S., Vieira, D., Bingner, R., Yoder, D., and Altinakar, M.: Modeling Agricultural Sheet, Rill and Ephemeral Gully Erosion, in: ICHE 2014. Proceedings of the 11th International Conference on Hydroscience & Engineering, Karlsruhe, 1119–1126, 2014. a

Ehlschlaeger, C.: Using the AT Search Algorithm to Develop Hydrologic Models from Digital Elevation Data, in: Proceedings of International Geographic Information Systems (IGIS) Symposium '89, Baltimore, MD, 275–281, MArch 1989. a

Flanagan, D. C., Frankenberger, J. R., Cochrane, T. A., Renschler, C. S., and Elliot, W. J.: Geospatial Application of the Water Erosion Prediction Project (WEPP) Model, T. ASABE, 56, 591–601, https://doi.org/10.13031/2013.42681, 2013. a, b, c

Foster, G., Meyer, L., and Onstad, C.: An erosion equation derived from basic erosion principles, T. ASAE, 20, 678–682, 1977. a

Goodwin, N. R., Armston, J., Stiller, I., and Muir, J.: Assessing the repeatability of terrestrial laser scanning for monitoring gully topography: A case study from Aratula, Queensland, Australia, Geomorphology, 262, 24–36, https://doi.org/10.1016/j.geomorph.2016.03.007, 2016. a

GRASS Development Team: GRASS GIS, available at: https://grass.osgeo.org, last access: 3 July 2019. a

Guertin, D. P., Goodrich, D. C., Burns, I. S., Korgaonkar, Y., Barlow, J., Sheppard, B. S., Unkrich, C., and Kepner, W.: Automated Geospatial Watershed Assessment Tool (AGWA), 120–130, https://doi.org/10.1061/9780784479322.012, 2015. a

Harmon, B. A.: r.sim.terrain, available at: https://github.com/baharmon/landscape_evolution, last access: 3 July 2019a. a, b, c, d

Harmon, B. A.: Landscape Evolution Dataset, available at: https://github.com/baharmon/landscape_ evolution_dataset, last access: 3 July 2019b. a, b, c

Harmon, B. A.: Landscape Evolution Tutorial, available at: https://github.com/baharmon/landscape_evolution/blob/master/tutorial.md, last access: 3 July 2019c. a, b

Harmon, B. A.: Landscape Evolution Data Log, available at: https://github.com/baharmon/landscape_evolution_dataset/blob/master/nc_spm_evolution/DATA.md, last access: 3 July 2019d. a

Harmon, B. A.: Landscape Evolution Repository, available at: https://osf.io/tf6yb/, last access: 3 July 2019e. a, b

Harmon, B. A.: r.sim.terrain, available at: https://grass.osgeo.org/grass76/manuals/addons/r.sim.terrain.html, last access: 3 July 2019f. a, b

Hobley, D. E. J., Adams, J. M., Nudurupati, S. S., Hutton, E. W. H., Gasparini, N. M., Istanbulluoglu, E., and Tucker, G. E.: Creative computing with Landlab: an open-source toolkit for building, coupling, and exploring two-dimensional numerical models of Earth-surface dynamics, Earth Surf. Dynam., 5, 21–46, https://doi.org/10.5194/esurf-5-21-2017, 2017. a, b

Hofierka, J., Mitášová, H., and Neteler, M.: Geomorphometry in GRASS GIS, in: Developments in Soil Science, edited by: Hengl, T. and Reuter, H. I., Elsevier, 33, 387–410, https://doi.org/10.1016/S0166-2481(08)00017-2, 2009. a

Huang, X. and Niemann, J. D.: Simulating the impacts of small convective storms and channel transmission losses on gully evolution, in: Military Geosciences in the Twenty-First Century, edited by Harmon, R. S., Baker, S. E., and McDonald, E. V., Geological Society of America, https://doi.org/10.1007/978-1-4020-3105-2_18, 2014. a

Jasiewicz, J. and Stepinski, T. F.: Geomorphons – a pattern recognition approach to classification and mapping of landforms, Geomorphology, 182, 147–156, https://doi.org/10.1016/j.geomorph.2012.11.005, 2013. a

Jeziorska, J., Mitasova, H., Petrasova, A., Petras, V., Divakaran, D., and Zajkowski, T.: Jeziorska, J., Mitasova, H., Petrasova, A., Petras, V., Divakaran, D., and Zajkowski, T.: Overland Flow Analysis Using Time Series of sUAS-derived Elevation Models, ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci., III-8, 159–166, https://doi.org/10.5194/isprs-annals-III-8-159-2016, 2016. a

Kasprak, A., Bransky, N. D., Sankey, J. B., Caster, J., and Sankey, T. T.: The effects of topographic surveying technique and data resolution on the detection and interpretation of geomorphic change, Geomorphology, 333, 1–15, https://doi.org/10.1016/j.geomorph.2019.02.020, 2019. a

Koco, Š.: Simulation of gully erosion using the SIMWE model and GIS, Landsurface Analysis, 17, 81–86, 2011. a

Landlab: http://landlab.github.io/, last access: 3 July 2019. a

Levine, J., Wegmann, K., Mitasova, H., Eads, C., Lyons, N., Harmon, B., McCarther, C., Peart, S., Oberle, N., and Walter, M.: Freshwater Bivalve Survey for Endangered Species Branch Fort Bragg, NC,North Carolina State University, Raleigh, NC, Tech. rep., 1–53, https://doi.org/10.13140/RG.2.2.17512.11521, 2018. a

Malik, I.: Dating of small gully formation and establishing erosion rates in old gullies under forest by means of anatomical changes in exposed tree roots (Southern Poland), Geomorphology, 93, 421–436, https://doi.org/10.1016/j.geomorph.2007.03.007, 2008. a

McDonald, K. W.: Military Foot Traffic Impact on Soil Compaction Properties, in: Studies in Military Geography and Geology, edited by: Caldwell, D. R., Ehlen, J., and Harmon, R. S., Springer Netherlands, Dordrecht, 229–242, https://doi.org/10.1007/978-1-4020-3105-2_18, 2004. a

Metz, M., Mitasova, H., and Harmon, R. S.: Fast Stream Extraction from Large, Radar-Based Elevation Models with Variable Level of Detail, in: Proceedings of Geomorphometry 2009, 237–242, 2009. a

Mitas, L. and Mitasova, H.: Distributed soil erosion simulation for effective erosion prevention, Water Resour. Res., 34, 505–516, https://doi.org/10.1029/97wr03347, 1998. a, b, c, d, e, f, g, h, i

Mitasova, H. and Mitas, L.: Interpolation by regularized spline with tension: I. Theory and implementation, Math. Geol., 25, 641–655, https://doi.org/10.1007/BF00893171, 1993. a

Mitasova, H. and Mitas, L.: Multiscale soil erosion simulations for land use management, in: Landscape erosion and evolution modeling, edited by: Harmon, R. S. and Doe, W. W., Springer, Boston, MA, chap. 11, 321–347, https://doi.org/10.1007/978-1-4615-0575-4_11, 2001. a

Mitasova, H., Hofierka, J., Zlocha, M., and Iverson, L. R.: Modelling topographic potential for erosion and deposition using GIS, Int. J. Geogr. Inf. Sci., 10, 629–641, https://doi.org/10.1080/02693799608902101, 1996. a, b, c, d, e, f

Mitasova, H., Thaxton, C., Hofierka, J., McLaughlin, R., Moore, A., and Mitas, L.: Path sampling method for modeling overland water flow, sediment transport, and short term terrain evolution in Open Source GIS, Dev. Water Sci., 55, 1479–1490, https://doi.org/10.1016/S0167-5648(04)80159-X, 2004. a, b, c, d, e

Mitasova, H., Mitas, L., and Harmon, R. S.: Simultaneous spline approximation and topographic analysis for lidar elevation data in open-source GIS, IEEE Geosci. Remote S., 2, 375–379, https://doi.org/10.1109/LGRS.2005.848533, 2005. a

Mitasova, H., Barton, M., Ullah, I., Hofierka, J., and Harmon, R.: 3.9 GIS-Based Soil Erosion Modeling, in: Treatise on Geomorphology, edited by: Shroder, J. F., Elsevier, San Diego, California, USA, chap. 3.9, 228–258, https://doi.org/10.1016/B978-0-12-374739-6.00052-X, 2013. a, b, c, d, e, f

Moore, I. and Burch, G.: Modeling Erosion and Deposition: Topographic Effects, Transactions of the American Society of Agricultural Engineers, 29, 1624–1640, 1986. a, b

North Carolina Risk Management Office: QL2 / QL1 LiDAR Collection, Tech. rep., 1–2, available at: https://sdd.nc.gov/sdd/docs/LiDARSummary.pdf (last access: 3 July 2019), 2018. a

Panagos, P., Ballabio, C., Borrelli, P., Meusburger, K., Klik, A., Rousseva, S., Tadić, M. P., Michaelides, S., Hrabalíková, M., Olsen, P., Aalto, J., Lakatos, M., Rymszewicz, A., Dumitrescu, A., Beguería, S., and Alewell, C.: Rainfall erosivity in Europe, Sci. Total Environ., 511, 801–814, https://doi.org/10.1016/j.scitotenv.2015.01.008, 2015. a

Panagos, P., Borrelli, P., Meusburger, K., Yu, B., Klik, A., Jae Lim, K., Yang, J. E., Ni, J., Miao, C., Chattopadhyay, N., Sadeghi, S. H., Hazbavi, Z., Zabihi, M., Larionov, G. A., Krasnov, S. F., Gorobets, A. V., Levi, Y., Erpul, G., Birkel, C., Hoyos, N., Naipal, V., Oliveira, P. T. S., Bonilla, C. A., Meddi, M., Nel, W., Al Dashti, H., Boni, M., Diodato, N., Van Oost, K., Nearing, M., and Ballabio, C.: Global rainfall erosivity assessment based on high-temporal resolution rainfall records, Scientific Reports, 7, 4175, https://doi.org/10.1038/s41598-017-04282-8, 2017. a, b

Perroy, R. L., Bookhagen, B., Asner, G. P., and Chadwick, O. A.: Comparison of gully erosion estimates using airborne and ground-based LiDAR on Santa Cruz Island, California, Geomorphology, 118, 288–300, https://doi.org/https://doi.org/10.1016/j.geomorph.2010.01.009, 2010. a

Renard, K. G., Foster, G. R., Weesies, G. A., McCool, D. K., and Yoder, D. C.: Predicting soil erosion by water: a guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE), US Government Printing Office, Washington, DC, Tech. Rep. No. 703, 1–384, available at: https://www.ars.usda.gov/ARSUserFiles/64080530/rusle/ah_703.pdf, (last access: 3 July 2019), 1997. a

Ad de Roo, A. P. J., Wesseling, C. G., Jetten, V. G., and Ritsema, C.: LISEM: A physically-based hydrological and soil erosion model incorporated in a GIS, in: Application of geographic information systems in hydrology and water resources management, edited by: Kovar, K. and Nachtnebel, H. P., Wallingford, UK, IAHS Publ. no. 235, 395-403, 1996. a

Salles, T.: eSCAPE: Regional to Global Scale Landscape Evolution Model v2.0, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2019-126, in review, 2019. a

Schoorl, J., Veldkamp, A., and Bouma, J.: Modeling Water and Soil Redistribution in a Dynamic Landscape Context, Soil Sci. Soc. Am. J., 66, 1610–1619, 2002. a

Schoorl, J. M., Sonneveld, M. P. W., and Veldkamp, A.: Three-dimensional landscape process modelling: the effect of DEM resolution, Earth Surf. Proc. Land., 25, 1025–1034, https://doi.org/10.1002/1096-9837(200008)25:9<1025::AID-ESP116>3.0.CO;2-Z, 2000. a

Shobe, C. M., Tucker, G. E., and Barnhart, K. R.: The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution, Geosci. Model Dev., 10, 4577–4604, https://doi.org/10.5194/gmd-10-4577-2017, 2017. a

Sorrie, B. A.: An Inventory of the Significant Natural Areas of Hoke County, North Carolina, North Carolina Natural Heritage Program, Tech. rep., 1–125, available at: http://digital.ncdcr.gov/cdm/ref/collection/p249901coll22/id/190173 (last access: 3 July 2019), 2004. a

Sorrie, B. A., Gray, J. B., and Crutchfield, P. J.: The Vascular Flora of the Longleaf Pine Ecosystem of Fort Bragg and Weymouth Woods, North Carolina, Castanea, 71, 129–161, https://doi.org/10.2179/05-02.1, 2006. a

Starek, M. J., Mitasova, H., Hardin, E., Weaver, K., Overton, M., and Harmon, R. S.: Modeling and analysis of landscape evolution using airborne, terrestrial, and laboratory laser scanning, Geosphere, 7, 1340–1356, https://doi.org/10.1130/GES00699.1, 2011. a, b

Telling, J., Lyda, A., Hartzell, P., and Glennie, C.: Review of Earth science research using terrestrial laser scanning, Earth-Sci. Rev., 169, 35–68, https://doi.org/10.1016/j.earscirev.2017.04.007, 2017. a

Temme, A., Schoorl, J., Claessens, L., and Veldkamp, A.: 2.13 Quantitative Modeling of Landscape Evolution, in: Treatise on Geomorphology, edited by: Shroder, J. F., Academic Press, San Diego, 2, 180–200, https://doi.org/10.1016/B978-0-12-374739-6.00039-7, 2013.  a

Thaxton, C. S.: Investigations of grain size dependent sediment transport phenomena on multiple scales, PhD thesis, North Carolina State University, available at: http://www.lib.ncsu.edu/resolver/1840.16/3339 (last access: 3 july 2019), 2004. a, b, c

Thomas, J. T., Iverson, N. R., Burkart, M. R., and Kramer, L. A.: Long-term growth of a valley-bottom gully, Western Iowa, Earth Surf. Proc. Land., 29, 995–1009, https://doi.org/10.1002/esp.1084, 2004. a

Tucker, G., Lancaster, S., Gasparini, N., and Bras, R.: The channel-hillslope integrated landscape development model (CHILD), in: Landscape erosion and evolution modeling, Springer, Boston, MA, 349–388, https://doi.org/10.1007/978-1-4615-0575-4_12, 2001. a, b

Tucker, G. E. and Hancock, G. R.: Modelling landscape evolution, Earth Surf. Proc. Land., 35, 28–50, https://doi.org/10.1002/esp.1952, 2010. a

Tucker, G. E. and Slingerland, R. L.: Erosional dynamics, flexural isostasy, and long-lived escarpments: A numerical modeling study, J. Geophys. Res., 99, 12229–12243, https://doi.org/10.1029/94JB00320, 1994. a

Webb, R. and Wilshire, H.: Environmental Effects of Off-Road Vehicles: Impacts and Management in Arid Regions, Environmental Management Series, Springer New York, https://doi.org/10.1007/978-1-4612-5454-6, 1983. a

Willgoose, G.: Mathematical Modeling of Whole Landscape Evolution, Annu. Rev. Earth Pl. Sc., 33, 443–459, https://doi.org/10.1146/annurev.earth.33.092203.122610, 2005. a

Wischmeier, W. H., Smith, D. D., Science, U. S., Administration, E., and Station, P. U. A. E.: Predicting Rainfall Erosion Losses: A Guide to Conservation Planning, Washington, D.C., Tech. rep., 537, 1–58, available at: https://naldc.nal.usda.gov/download/CAT79706928/PDF (last access: 3 July 2019), 1978. a

Yang, S., Guan, Y., Zhao, C., Zhang, C., Bai, J., and Chen, K.: Determining the influence of catchment area on intensity of gully erosion using high-resolution aerial imagery: A 40-year case study from the Loess Plateau, northern China, Geoderma, 347, 90–102, https://doi.org/10.1016/j.geoderma.2019.03.042, 2019. a

Yin, S., A. Nearing, M., Borrelli, P., and Xue, X.: Rainfall Erosivity: An Overview of Methodologies and Applications, Vadose Zone J., 16, 12, https://doi.org/10.2136/vzj2017.06.0131, 2017. a

Zahra, T., Paudel, U., Hayakawa, Y., and Oguchi, T.: Knickzone Extraction Tool (KET) – A New ArcGIS toolset for automatic extraction of knickzones from a DEM based on multi-scale stream gradients, Open Geosci., 9, 73–88, https://doi.org/10.1515/geo-2017-0006, 2017. a

Publications Copernicus
Download
Short summary
The numerical model, r.sim.terrain, simulates how overland flows of water and sediment shape topography over short periods of time. We tested the model by comparing runs of the simulation against a time series of airborne lidar surveys for our study landscape. Through these tests, we demonstrated that the model can simulate gully evolution including processes such as channel incision, channel widening, and the development of scour pits, rills, and depositional ridges.
The numerical model, r.sim.terrain, simulates how overland flows of water and sediment shape...
Citation