r.sim.terrain: a dynamic landscape evolution model

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, GIS-based topographic evolution model uses empirical models for soil erosion at watershed to regional scales and a physics-based model for shallow overland water flow and soil erosion at subwatershed scales to compute short-term topographic change. This either steady state 5 or dynamic model simulates 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 Patterson Branch subwatershed on the Fort Bragg military installation in North Carolina, r.sim.terrain can realistically simulate the development of fine-scale morphological features including ephemeral gullies, rills, and hillslopes. Applications include land management, erosion control, landscape planning, and landscape restoration. 10 Copyright statement. ...


Introduction
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 (Temme et al., 2013). Numerical landscape evolution models such as 15 the Channel-Hillslope Integrated Landscape Development (CHILD) model (Tucker et al., 2001) and SIBERIA (Willgoose, 2005) simulate steady state flows over long temporal scales. Landlab, 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 20 soil erosion models or landscape evolution models (see Tables 1-2). Furthermore 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 a. DEM derived from 2016 airborne lidar survey b. DEM generated by a transport limited simulation with SIMWE Figure 1. The digital elevation model (DEM) before (a.) and after (b.) simulated landscape evolution with r.sim.terrain. This simulation used the SIMWE model for a 120 min rainfall event with 25 mm hr −1 in a transport limited soil erosion regime at steady state. In the evolved DEM (b.) the gully channel has widened with depositional ridges forming along its thalweg. peak flows at steady state (see Table 2), short-term erosional processes like gully formation can be dynamic with significant morphological changes happening within minutes before flows reach steady state. A dynamic landscape evolution model is needed to study fine-scale spatial and short-term temporal erosional processes such as gully formation and the development of microtopography.
At the beginning of a rainfall event the overland water flow regime is dynamic -its depth changes at a variable rate over time 5 and space. If the intensity of rainfall continues to change throughout the event then the flow regime will remain dynamic. If, however, the overland flow reaches a peak rate then the hydrologic regime is considered to be at steady state. At steady state: ∂h(x, y, t) ∂t = 0 where: (x, y) is the position [m] t is the time [s] 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 waterflow 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 5 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 until water flow reaches steady state. 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. After the initial formation of the gully the soil 10 erosion regime may switch to erosion-deposition if the intensity of rainfall decreases. 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, the gully may fill with sediment.
Gully erosion rates and evolution can be monitored in the field or modeled on the computer. Field methods include den-15 drogeomorphology (Malik, 2008) 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 (UAS), airborne lidar, and terrestrial lidar (Starek et al., 2011;Bechet et al., 2016). With terrestrial lidar, airborne lidar and UAS photogrammetry there is now sufficient resolution topographic data to morphometrically analyze and numerically model fine-scale landscape evolution in GIS including processes 20 such as gully formation and the development of microtopography. Gully erosion has been simulated with the Revised Universal Soil Loss Equation Version 2 (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 (Koco, 2011;Mitasova et al., 2013). Now numerical landscape evolution models that can simulate steady state and dynamic flow regimes and can dynamically switch between soil erosion regimes are needed to study 25 fine-scale spatial and short-term temporal 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 dynamic flow regimes, dynamically switch between soil erosion regimes, and simulate the evolution of fine-scale 30 morphological features such as ephemeral gullies (Figure 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 dynamic flow regimes, and studying the role of dynamic flow regimes in fine-scale morphological change. r.sim.terrain was tested with a subwatershed scale (450 m 2 ) case study and the simulations were compared against a time-series of airborne lidar surveys.
2 r.sim.terrain 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 (RUSLE 3D) model, or the Unit Stream Power Erosion Deposition (USPED) model. SIMWE is a physics-based simulation that uses a Monte Carlo path sampling method to 5 solve the water and sediment flow equations for detachment limited, transport limited, and variable erosion-deposition soil erosion regimes (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 sediment flows in detachment capacity limited soil erosion regimes (Mitasova et al., 1996). With RUSLE3D r.sim.terrain uses an event-based erosivity factor, the slope, the flow accumulation, and a 3D topographic factor to model sediment flow.

10
USPED is an empirical equation for net erosion and deposition in transport capacity limited soil erosion regimes. With USPED r.sim.terrain uses an event-based erosivity factor, the slope and aspect, the flow accumulation, and a 3D topographic factor to model erosion-deposition as the the divergence of sediment flows. For each of the models topographic change is derived at each time step from the sediment flow or net erosion-deposition rate and gravitational diffusion. The r.sim.terrain model can simulate either steady state or dynamic flow regimes. During simulations with SIMWE r.sim.terrain can switch between 15 detachment limited, transport limited, and variable erosion-deposition soil erosion regimes.
The r.sim.terrain model can simulate the evolution of gullies including processes such as knickpoint migration, channel incision, channel widening, aggradation, and scour pool and depositional ridge formation along the thalweg 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 water-20 sheds less than ten square kilometers with SIMWE to regional watersheds of hundreds of square kilometers with USPED or that uses a path sampling method to solve the continuity and momentum equations with a 2D diffusive wave approximation 30 (Mitas and Mitasova, 1998;Mitasova et al., 2004). SIMWE has been implemented in GRASS GIS as the modules r.sim.water and r.sim.sediment.   (Mitasova et al., 1996) USPED watershed continuous raster map algebra (Mitasova et al., 1996) SIMWE watershed event -raster GRASS modules (Mitas and Mitasova, 1998) continuous GeoWEPP watershed continuous raster ArcGIS module (Dennis C. Flanagan et al., 2013) AGWA watershed event -vector ArcGIS module (Guertin et al., 2015) continuous In SIMWE mode for each time step r.sim.terrain determines the soil erosion regime, simulates water and sediment flows, and then evolves the topography. In an variable erosion-deposition regime the model computes the partial derivatives of the topography, simulates shallow water flow and erosion-deposition, and then evolves the topography based on the erosiondeposition rate and gravitational diffusion. The same process is used in a transport capacity limited regime, except that the topography is evolved based on the transport limited erosion-deposition rate and gravitational diffusion. In a detachment 5 capacity limited regime the model instead computes the partial derivatives of the topography, simulates shallow water flow and sediment flow, and then evolves the topography based on the sediment flow rate and gravitational diffusion. The model simulates dynamic landscape evolution when the time step is less than the travel time for a drop of water or a particle of sediment to cross the landscape. With longer time steps the model simulates steady state dynamics. depends on soil and landcover properties. The detachment capacity is the maximum potential detachment rate by overland flow, while the sediment transport capacity is the maximum potential sediment flow rate. When rainfall intensity is very high , then the regime is detachment capacity limited. When rainfall intensity is not very high (i r < 60mm hr −1 ) and σ is high (σ ≥ 100m −1 ), then the regime is transport capacity limited. When rainfall intensity is not very high (i r < 60mm hr −1 ) and σ is neither high nor low (0.01m −1 < σ < 100m −1 ), then there is an variable Tc is the sediment transport capacity [kg m −1 s −1 ]

Shallow water flow
The SIMWE model simulates shallow overland water flow controlled by spatially variable topographic, soil, landcover, and 15 rainfall parameters by solving the continuity and momentum equations for steady state water flow with a path sampling method ( Fig. 2a). Shallow water flow q can be approximated by the bivariate form of the St. Venant equation: where: x Diffusive wave effects can be approximated so that water can flow through depressions by integrating a diffusion term ∝ ∇ 2 [h 5/3 ] into the solution of the continuity and momentum equations for steady state water flow. This equation is solved using a Green's function Monte Carlo path sampling method. where: ε is a spatially variable diffusion coefficient.

Sediment flow
In SIMWE the sediment flow rate q s is estimated as a function of water flow and sediment concentration (Mitas and Mitasova,5 1998) (Fig. 2b): where: qs is the sediment flow rate per unit width [kg m −1 s −1 ] ρs is sediment mass density [kg m −3 ].

Erosion-deposition
In SIMWE the net erosion-deposition rate is estimated using the bivariate form of sediment continuity equation to model sediment storage and flow based on effective sources and sinks (Fig. 2c). Net erosion-deposition d s -the difference between sources and sinks -is approximated by the steady state sediment flow equation with diffusion: where: ds is net erosion-deposition [kg m −2 s −1 ].

Landscape evolution 20
The simulated change in elevation ∆z due to water erosion and deposition is a function of the change in time, the net erosiondeposition rate, and the sediment mass density (Mitasova et al., 2013): In a detachment limited erosion regime the simulated change in elevation ∆z is a function of the change in time, the sediment flow rate, and the mass of water carried sediment per unit area (Mitasova et al., 2013): Gravitational diffusion is then applied to the evolved topography to simulate the settling of sediment particles. The simulated change in elevation ∆z due to gravitational diffusion is a function of the change in time, the sediment mass density, the gravitational diffusion coefficient, and topographic divergence -i.e. the sum of the second order derivatives of elevation (Thaxton, 2004): where: ∇ is the topographic divergence [m −1 ].

10
The Revised Universal Soil Loss Equation for Complex Terrain (RUSLE3D) is an empirical equation 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 Universal Soil Loss Equation (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. As 15 an empirical equation the predicted soil loss is spatially and temporally averaged. 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, 3D topography with converging and diverging flows. In RUSLE3D the topographic potential for erosion at any given point is 20 represented by a 3D topographic factor LS 3D , 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 the rainfall intensity at each time step. For each time step this model computes the parameters for RUSLE3D -an eventbased erosivity factor, the slope of the topography, the flow accumulation, and the 3D topographic factor -and then solves the RUSLE3D equation for sediment flow. The sediment flow is used to simulate landscape evolution in a detachment capacity 25 limited soil erosion regime.

Event-based 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 I 30 (Brown and 30 Foster, 1987;Renard et al., 1997). In this model, however, the erosivity factor is derived at each time step as a function of kinetic energy, rainfall volume, rainfall intensity, and time. First rain energy is derived from rainfall intensity (Brown and where: Then the event-based erosivity index R e is calculated as the product of unit rain energy, rainfall volume, rainfall intensity, and time: where: Re is the event-based erosivity index [MJ mm ha −1 hr −1 ] vr is the rainfall volume [mm] derived from vr = ir tr tr is the time interval [s].

Flow accumulation
The upslope contributing area per unit width is determined by flow accumulation times grid cell width (Fig. 2d). Flow ac-15 cumulation is calculated using a multiple flow direction algorithm (Metz et al., 2009) based on A T least cost path searches (Ehlschlaeger, 1989). The multiple flow direction algorithm implemented in GRASS GIS as the module r.watershed is computationally efficient and can navigate nested depressions and other obstacles.

3D topographic factor
The 3D topographic factor LS 3D (x, y) is calculated as a function of the flow accumulation, representing the upslope contribut-20 ing area, and the slope (Fig. 2e). 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.

Sediment flow
Sediment flow is a function of the event-based erosivity factor, the soil erodibility factor, the 3D topographic factor, cover 5 factor, and the prevention measures factor (Fig. 2f): where: Re is the event-based erosivity factor [MJ mm ha −1 hr −1 ]

10
K is the soil erodibility factor [ton ha hr ha −1 MJ −1 mm −1 ] LS3D is the dimensionless topographic (length-slope) factor C is the dimensionless land cover factor P is the dimensionless prevention measures factor.

15
For RUSLE3D the simulated change in elevation ∆z is derived from equation 8 for landscape evolution in an detachment limited soil erosion regime and then equation 9 for the settling of sediment particles due to gravitational diffusion.

Unit streampower erosion deposition model
The Unit Stream Power Erosion Deposition (USPED) model estimates net erosion-deposition as the divergence of sediment flow in transport capacity limited soil erosion regimes. At transport capacity shallow flows of water are carrying as much 20 sediment possible -more sediment is being detached than can be transported. As a transport capacity limited model USPED predicts erosion where transport capacity increases and deposition where transport capacity decreases. The influence of topography on erosion and deposition in USPED is represented by a topographic sediment transport factor, while the influence of soil and landcover are represented by factors adopted from USLE and RUSLE (Mitasova et al., 1996). Net erosion-deposition is estimated by computing the event-based erosivity factor (R e ) using Eq. 11, the slope and aspect of the topography, the flow 25 accumulation with a multiple flow direction algorithm, the topographic sediment transport factor, the sediment flow at transport capacity, and the divergence of the sediment flow.
The 3D topographic factor (Eq. 12) for RUSLE3D is adapted to represent the topographic sediment transport factor (LST ) -the topographic component of overland flow at sediment transport capacity: 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 5 n is an empirical coefficient.
The sediment flow at transport capacity is a function of the event-based rainfall factor, the soil erodibility factor, the topographic component of overland flow, the landcover factor, and the prevention measures factor: where: T is sediment flow at transport capacity [kg m −1 s −1 ] Re is the event-based rainfall factor [MJ mm ha −1 hr −1 ] K is the soil erodibility factor [ton ha hr ha −1 MJ −1 mm −1 ] C is the dimensionless land cover factor 15 P is the dimensionless prevention measures factor.
Net erosion-deposition at transport capacity is estimated as the divergence of sediment flow: α is the aspect of the topography [ • ].
With USPED the simulated change in elevation ∆z is derived from equation 7 for landscape evolution and then equation 9 for the settling of sediment particles due to gravitational diffusion. 25

Case study
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 Wilshire, 1983;30 McDonald, 2004). Gullies -ephemeral channels with steep headwalls that incise into unconsolidated soil to depths of meters -  To test the effectiveness of the different models in r.sim.terrain we compared the simulated evolution of a highly eroded , and USPED -were tested in steady state and dynamic modes for constant rainfall, design storms, and recorded rainfall.

Patterson Branch Creek
With 650 km 2 of land Fort Bragg is the largest military installation in the US and has extensive areas of bare, erodible soils on (Aristida stricta) on Blaney and Gilead loamy sands (Sorrie, 2004). 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 m 2 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 (Figure 4c). 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.

5
Timeseries of digital elevations models and landcover maps for the study landscape were generated from lidar pointclouds and orthophotography (Figure 4a-c). The digital elevations models for 2004, 2012, and 2016 were interpolated at 0.3 m resolution using the regularized spline with tension function (Mitasova and Mitas, 1993;Mitasova et al., 2005)  Spatially variable soil erosion factors -k-factor, c-factor, mannings, and runoff rates -were derived from the landcover and soil maps. The dataset for this study is hosted at https://github.com/baharmon/landscape_evolution_dataset under the ODC Open Database License (ODbL). The data is derived from publicly available data from the US Army, USGS, USDA, Wake County GIS, NC Floodplain Mapping Program, and the NC State Climate Office. 15 We used the geomorphons method of automated landform classification based on the openness of terrain (Jasiewicz and Stepinski, 2013)

Simulations
We ran a sequence of r.sim.terrain simulations for the Patterson Branch Creek subwatershed study area to test dynamic and steady state flow regimes in the SIMWE, RUSLE3D, and USPED models (Table 3). RUSLE3D was used to simulate 120 min events with rainfall intensities of 50 mm hr −1 for detachment capacity limited soil erosion regimes for both dynamic and steady state flow regimes using RUSLE3D (Figure 5a-c). USPED was used to simulate 120 min events with rainfall intensities  (Table 3).

10
The dynamic RUSLE3D simulation deepened the main channel of the gully, while the dynamic USPED simulation eroded the banks of the gully and deposited in channels causing the gully grow wider and shallower ( Figure 5). As a detachment capacity limited model RUSLE3D's results were dominated by erosion and thus negative elevation change. RUSLE3D carved a deep incision in the main gully channel where water and sediment flow accumulated (Figure 5c). As a transport capacity limited model USPED generated a distributed pattern with both erosion and deposition and thus negative and positive elevation change.

15
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 (Figure 5f).
The steady state SIMWE simulations predicted more realistic patterns of landscape evolution ( Figure 6). Given the presence of an active gully with ridges along its banks, this landscape is dominated by a detachment limited soil erosion regime. The detachment limited SIMWE simulation generated the morphological features -the deeply incised gully 25 channels, scour pits, and ridges along the channels -characteristic of its erosion regime, realistically simulating landscape evolution at the scale of a subwatershed. The erosion-deposition and transport limited SIMWE simulations also generated the morphological processes and features that would be expected in these regimes -gradual aggradation and the formation of a depositional ridge along the thalweg of the channel.
While RUSLE3D and USPED produced less realistic patterns of landscape evolution than SIMWE, these models were much 30 faster and still generated 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 at regional scales, i.e. for landscapes greater than 10 km 2 . RUSLE3D for example has been used to model erosion for the entire 650 km 2 Fort Bragg installation at 9 m resolution (Levine et al., 2018).

Conclusions
The short-term landscape evolution model r.sim.terrain can realistically simulate the development of gullies, rills, and hillslopes by overland water erosion for a range of hydrologic and soil erosion regimes. The landscape evolution model was tested with a series of simulations for different hydrologic and soil erosion regimes for a highly eroded sub-watershed on Fort Bragg with an active gully. For each regime it generated the morphological processes and features expected. The physics-based 5 SIMWE model realistically simulated short-term topographic change for steady state hydrologic regimes at sub-watershed to watershed scales. For detachment limited soil erosion regimes it simulated morphological processes including channel incision, channel widening, and the development of knickzones, rills, and scour pits. For transport limited and variable erosiondeposition regimes, it simulated processes such as channel aggradation, scouring, and the development of depositional ridges along the thalweg. The empirical RUSLE3D and USPED models approximated short-term topographic change at watershed to 10 regional scales. For detachment limited soil erosion regimes RUSLE3D simulated channel incision, while for transport limited regimes USPED simulated channel widening and filling. Since it is a GIS-based model that realistically simulates fine-scale morphological processes and features, r.sim.terrain 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.