The MARSSIM model is a landform evolution model primarily focuses on relatively long temporal scales (relative to the timescale for noticeable landform change) through fluvial and mass wasting processes.

## Model introduction

The program is designed be computationally efficient such that individual runs can be done on a modern microcomputer in no more than a few tens of hours. The more recent additions to the model have focused on processes relevant to planetary landscapes, including lava flows, groundwater seepage and sapping, impact cratering, surface-normal accretion and ablation, and volatile redistribution by radiation-induced sublimation and recondensation. Individual process formulations vary from completely heuristic to modestly mechanistic. Important limitations for some potential applications are the assumption of a single representative bed material grain size in the fluvial system and no tracking of internal stratigraphy of sedimentary deposits.

## Model parameters

Parameter Description Unit
Input directory - -
Variable information filename - -
Parameter Description Unit
Random seed an integer random seed -
Switch for new simulation If set to 1, start a new simulation. Otherwise, continue from a saved state -
Switch for fluvial and slope modeling If set to 1, fluvial and slope processes are modeled -
Switch for impact crater modeling If set to 1, impact craters are modeled -
Switch for lava inundation If set to 1, lava inundation is modeled -
Switch for eolian deposition If set to 1, eolian transport is modeled -
Switch for ocean processes If set to 1, ocean processes is modeled -
Switch for accretion and ablation If set to 1, model accretion and ablation -
Switch for lake evaporation If set to 1, model runoff and evaporation with lakes -
Number of columns for DEM - -
Number of rows for DEM - -
Number of vertical cells - -
Grid cell size Grid cell size arbitrary units. See conversion factor below -
Vertical scale units Elevation units (vertical). Must as same as grid cell size units (horizontal) -
Convert-to-meters factor Factor to convert grid cell size and vertical scale units to meters -
Switch for variable erosion rate If set to 1, erosion rate will vary in time in the specified manner. Otherwise, model parameters are constant -
Boundary lowering rate Rate at which the boundary is lowering over time (meter/year). To be used with a nonperiodic vertical boundary -
Channel timestep scale factor channel process timestep scale factor -
Default channel time step channel process timestep years
Maximum timestep Maximum allowed timestep years
Minimum timestep Minimum allowed timestep years
Sediment timestep scale factor - -
Mass wasting timestep scale factor - -
I-coordinate for center of debug window Used for debug mode -
J-coordinate for center of debug window Used for debug mode -
Size of debug window, I-direction Used for debug mode -
Size of debug window, J-direction Used for debug mode -
Switch to print debug info If set to 1, print info in specified debug window, Affects performance -
Switch to calculate morphometric statistics If set to 1, calculate morphometric statistics -
Switch to use horizontal lower boundary If set to 1, southern (lower ) boundary is level. Do not use with a periodic lower boundary -
Switch to use non-eroding lower boundary If set to 1, use a non-eroding lower/southern boundary can be used in conjunction with a horizontal lower boundary -
Switch for periodic Y axis If set to 1, use periodic boundary condition along the Y-axis -
Switch for periodic X axis If set to 1, use periodic boundary condition along the X-axis -
Switch to allow flow across all boundaries All boundaries are potential flow exits and are non-eroding -
Parameter Description Unit
Starting iteration Generally this is zero for new simulatins -
Present time Generally this is zero for new simulatins years
Maximum iteration The number of timesteps or iterations to be used in the simulation -
Maximum simulation time The maximum time that the simulation is allowed to run years
Iterations between prints The number of iterations (or timesteps) between writing out a lot of info -
Simulation progress interval How often to summarize simulation progress. Generally should be the same as elevinterval -
Recalculate gradient interval Number of interations before recalculating gradients and contributing areas (Currently hardwired) -
Write change interval Number of iterations before printing out information on the rate of change in gradients and flow directions (Currently hardwired) -
Morphometric channel threshold This has to do with the definition of where channels start when printing out morphometric statistics -
Matrices to output What matrices to output at Elevation_print_interval times. Each digit represents a flags are (1) elevation; (2) relative elevation; (3) D8 gradient; (4) erode_slope; (5) cfne; (6) cfnw; (7) cfw; (8) cfn; (9) cfne; (10) drainage_area; (11) relative_gradient; (12) regolith and relative resistance. Generally, only first flag is set (elevation) -
Switch for absolute or relative elevations If set to 1, write out absolute vs. relative elevations -
Interval to print out shaded relief images Shaded relief images are typically generated and saved every 50 to 100 iterations -
Divergence Interval - -
Parameter Description Unit
Switch for wet or dry depressions If set to 1, all depressions fill with water and overflow. Otherwise, all depressions are infinite sinks for water. Most Mars simulations have used 0, but some simulatins with non-periodic BCs use 1 -
Switch for rescale discharges IF set to 1, then effective discharges are rescaled as a power function of contributing area (e.g., Q = K*A**E, where E is typically in [0.5,1.0]) -
Discharge constant Used with option to rescale discharges. In that case, Discharge = discharge_constant*Area**discharge_exponent m3 / s
Discharge exponent Used with option to rescale discharges. In that case, Discharge = discharge_constant*Area**discharge_exponent m3 / s
Area factor Generally not used -
Rain depth - -
Rain standard deviation - -
Switch for divergence-dependent runoff If set to 1, runoff depth functionally depends on topographic divergence (VARYIELDUSE; 0, or 1), usually not set -
Mean convergence runoff Only used if 'Divergence dependent runoff' is toggled on -
High convergence runoff Only used if 'Divergence dependent runoff' is toggled on -
Low convergence runoff Only used if 'Divergence dependent runoff' is toggled on -
Iterations between recalculating evaporation rates Only used if 'Lake evaporation' is toggled on -
Iterations between stocahstic changes in evaporation rates Only used if 'Lake evaporation' is toggled on -
Mean yearly evaporation depth - meters
Std. deviation of mean evaporation depth - meters
Switch to model pelagic deposition If set to 1, deposition of suspended sediment in enclosed basins is modeled (IMODEL_PELAGIC_DEPOSITION, 0 or 1) -
Wash load fraction Fraction of sediment delivered to basins that is deposited as pelagic sediment -
Amount of post-deposition A diffusivity that determines how much post-deposition diffusion occurs in recently-deposited pelagic sediment -
Switch to model influent rivers IF set to 1, incoming rivers are modeled -
Parameter Description Unit
Switch to model fluvial detachment This should only set to 0 when slope processes are being modeled and fluvial erosion, transport and deposition are not -
Switch for explicit or implicit channel bed If set to 1, use EXPLICIT method to model regolith depth (including bare rock). Otherwise use IMPLICIT method -
Model to use for bedrock channel 0 for Howard model, 1 for Sklar saltation, 2 for WWhipple saltation -
Rock tensile strength Tensible strength of bedrock if saltation model is used -
Bedrock discharge exponent In the Howard detachment model, the rate of bedrock erosion is computed as: bedrock erodibility * (Discharge ** bedrock discharge exponent) * ( Gradient ** bedrock gradient exponent) -
Bedrock gradient exponent In the Howard detachment model, the rate of bedrock erosion is computed as: bedrock erodibility * (Discharge ** bedrock discharge exponent) * ( Gradient ** bedrock gradient exponent) -
Bedrock erodibility In the Howard detachment model, the rate of bedrock erosion is computed as: bedrock erodibility * (Discharge ** bedrock discharge exponent) * ( Gradient ** bedrock gradient exponent) -
Manning roughness parameter Manning's channel roughness parameter, N -
Gravitational constant 3.7 for Mars, 9.8 for Earth m/s2
Detachment critical shear A critical stress that must be exceeded in the Howard detachment model before transport occurs, usually set to 0 kg/(m*s2)
Regolith erodibility factor Determines how much more erodible weathered regolith is as compare to bedrock, minimum value is 1. -
Regolith critical shear factor Determine how much less the critical shear stress for detachment of regolith is as compared to bedrock. Generally set to the same value as the regolith erodibility factor -
Channel width constant Channel width is computed as: Width = channel width constant * ( Discharge ** channel width exponent). Generally, mean annual discharge is used -
Channel width exponent Channel width is computed as: Width = channel width constant * ( Discharge ** channel width exponent). Generally, mean annual discharge is used -
Switch for variable vegetation resistance Switch to determine whether variable vegetation modeling is used (ITAUAREA, 0 or 1) -
Vegetation upland resistance Only used if 'variable vegetation resistance' flag is set -
Vegetation channel resistance Only used if 'variable vegetation resistance' flag is set -
Vegetation area minimum Only used if 'variable vegetation resistance' flag is set -
Vegetation area maximum Only used if 'variable vegetation resistance' flag is set -
Parameter Description Unit
Switch to write sediment diagnostics If set to 1, detailed sediment routing information is printed (ISEDDEBUGT, 0 or 1) -
Switch for sediment transport If set to 1, alluvial channels, delta and fans are modeled (ISEDIMENT, 0 or 1) -
Switch for sediment routing If set to 1, sediment routing is modeled (ISEDIMENT, 0 or 1) -
Switch for sediment diffusion If set to 1, sediment diffusion is modeled (ISEDIMENT, 0 or 1) -
Switch for no-flux lower boundary If set to 1, a no-flux lower boundary is used (IREFLECT, 0 or 1) -
Switch for alluvial smoothing If set to 1, alluvial depositional surfaces (fans and deltas) are smoothed during deposition (SMOOTHSEDUSE, 0 or 1) -
Alluvial smoothing factor only used if 'alluvial rerouting' switch is used -
Switch to read alluvial thickness at start If set to 1, alluvium thickness is read at beginning of simulation (IREADALLUV, 0 or 1), generally 0 -
Switch to use an ocean If set to 1, an ocean is used, as when modeling delta (WATERUSE, 0 or 1), set to 1 if USEWET is set to 1 -
Ocean elevation If modeling a global ocean, set this to your assumed ocean level. Otherwise, set it to a very low value -
Sediment 1 exponent Used in sediment transport relationship -
Sediment 2 exponent Used in sediment transport relationship -
Effective discharge ratio Used in sediment transport relationship -
Sediment gradient exponent Used in sediment transport relationship -
Sediment transport exponent Used in sediment transport relationship -
Transport factor Generally, 40.0 for sand bed, 8.0 for gravel bed -
Flow fraction Set to 1 when modeling erosion by groundwater seepage -
Sediment porosity - -
Sediment specific gravity - -
Grain size Typical value are 0.0002 for sand and 0.02 for gravel meters
Transport critical dim shear Generally not changed -
Switch to use sticky sediment routing If set to 1, use 'sticky' VS. normal sediment. Used in an experiment to model birdfoot delta. Generally not used -
Sticky routing critical value Use in an experiment to model birdfoot delta. Generally, not used -
Bedload fraction Fraction of eroded sediment that is transported as bedload. Used for calculation of alluvial channel gradients -
Parameter Description Unit
Rock weathering rate Rate of bare rock weathering meters/year
Weather decay rate Determines how rapidly weathering rate decreases with increasing regolith thickness, based on terrestrial values -
Initial regolith thickness - meters
Switch for two-term weathering If set to 1, a humped weathering function is used (WEATHER2USE, 0 or 1) -
Weathering term 2 only used if 'two term weathering' flag is set -
Weathering decay 2 only used if 'two term weathering' flag is set -
Switch for seepage weathering If set to 1, seepage-dependent weathering is modeled (ISSEPAGEWEATHER, 0 or 1) -
Seepage weathering scale factor only used if 'Seepage weathering' flag is set -
Seepage weathering exponent only used if 'Seepage weathering' flag is set -
Critical bedrock gradient Maximum steepness of bedrock slope -
Weather multiplier This scales the rate of bedrock mass wasting -
Weather divergence If greater than 0, then rate of bedrock weathering will depend on topographic divergence -
Switch to read initial regolith thickness If set to 1, read initial regolith thickness from file (READREOLITH, 0 or 1) -
Parameter Description Unit
Switch for slope erosion modelling If set to 1, use slope erosion modelling (ISLOPEUSE, 0 or 1) -
Slope diffusity Determines relative rate of downslope creep (m2/yr) try values between 0.0002 and 0.02 -
Proportion of slope constant in alluvial creep Relative diffusivity of alluvium mass wasting to slope mass wasting. Generally set to 1 -
Switch to use critical slope gradient If use to 1, use a critical slope gradient (Critical_gradient_use, 0 or 1) -
Switch to use Roering mass wasting If set to 1, use the Roering formula for mass wasting. Otherwise, use the Howard formula. The Roering formula allows a threshold maximum slope gradient in regolith. The Howard formula uses linear creep -
Slope gradient exponent Used for Howard two-term relationship -
Parameter Description Unit
Switch to use bistable fluvial erosion If set to 1, use bistable fluvial erosion (HIGHRATEUSE, 0 or 1) used for modelling gully development -
Switch to use bistable bedrock If set to 1, use bistable bedrock (BEDROCKHIGH, 0 or 1) used for modeling gully development -
Low erosion threshold Only used with bistable erosion and gully development -
High erosion threshold Only used with bistable erosion and gully development -
Erosion rate change lag Only used with bistable erosion and gully development -
Bistable critical shear Only used with bistable erosion and gully development -
Bistable runoff factor Only used with bistable erosion and gully development -
Bistable bedrock erodibility Only used with bistable erosion and gully development -
Parameter Description Unit
Switch to use random critical shear stress If set to 1, use random critical shear stress (RANDTHRESHUSE, 0 or 1) -
Switch to use random discharge If set to 1, use random discharge (RANDDISCHUSE, 0 or 1) -
Critical shear variablity Scale factor for random critical shear stress -
Discharge coeff variation Standard error of discharge, assuming a lognormal distribution -
Omega weight Governs inherientance from past value, similar to Hurst parameter set to 1 for no inheritance -
Parameter Description Unit
Switch to use 3D variable resistance weathering If set to 1, use 3D variable resistance weathering (SLOPEVARUSE, 0 or 1) -
Diffusivity variability - -
Switch to use a resistant surface layer If set to 1, use a resistant surface layer (CRUSTUSE, 0 or 1). Experimental, for simulating the effect of development of a "duricrust" -
Surface layer thickness Crust thickness, using the same units as cell size -
Surface layer resistance Crust relative shear stress factor, greater than 1 for a resistant crust -
Vertical resistance scale factor Ratio of horizontal to vertical scale's (1.0 if equal). Only used for 3D varing rock resistance -
Parameter Description Unit
Switch to use rock deformation DEFORMUSE, 0 or 1, not implemented in current version -
Deformation scale factor Not implemented in current version -
Parameter Description Unit
Switch to model groundwater seepage If set to 1, groundwater seepage is modeled (SEEPUSE, 0 or 1) -
Seepage iteration interval Number of iterations between recalculation of groundwater flow -
Switch for permeability rescaling If set to 1, permeability in exponential vertical decay of permeability is scaled to the present land surface elevation. Otherwise it is scaled to the original land surface elevation (at the start of the simulation) -
Switch to show groundwater calculations If set to 1, show groundwater calculations -
Switch for exponential permeability decay If set to 1, permeability is constant through an aquifer of constant thickness. Otherwise, it decays with depth below the surface -
Yearly recharge rate - meters/year
Viscosity of water - -
Permeability, in darcies Permeability darcies
Groundwater depth scale For exponential permeability, this is the depth to the half value of permeability. For constant permeability, it is the acquiter thickness meters
Groundwater flow fraction The fraction of surface water flow that becomes groundwater flow, set to 1 for all groundwater flow -
Initial groundwater depth Initial depth beneath the land surface for calculation of steady state groundwater surface meters
Depth decay of permeability power Set this to 1 for normal exponential decay -
Max. groundwater iterations Max number of iterations, permitted in calculating steady state groundwater table -
Max. groundwater error Max residual error permitted in calculating steady state groundwater table -
Groundwater relaxation coeff Coefficient used for SOR (Successive Over Relaxation) -
Switch to use groundwater flux If set to 0, the groundwater flow term used for seepage calculation of surface flows and weathering rate is the groundwater flux divergence (Seepage rate to surface). Otherwise, it is the groundwater flow rate per unit aquifer width -
Switch to use seepage averaging If set to 1, the groundwater flow term used in further calculations is a 9 point average value of raw calculated value -
Parameter Description Unit
Eolian event probability Probability of an eolian deposition/erosion event per unit time (year or iteration) -
Eolian time increment Scales the overall Eolian Erosion/Deposition rate -
Switch to use total exposure If set to 1, all cells within the calculation window are used to compute the exposure index. Otherwise, only cells visible to the local cell are used -
Switch to use default eolian process If set to 1, eolian erosion and deposition occur normal to the topographic surface. Otherwise, erosion and deposition are modeled as vertical additions or subtractions -
Minimum eolian deposition rate This can be negative (If eolian erosion of exposed locations) -
Maximum eolian deposition rate Generally unity and scaled by Eolian_time_increment determines Eolian_constant_1 and Eolian_constant_2 -
Eolian parameter code Code that determines how next 2 parameters are used: (1) Exposure_10_percent, Exposure_90_percent; (2) Exposure_50_percent, Exposure_90_percent; (3) Zero_percent_exposure, Exposure_90_percent; (4) Rate0, exposure_50_percent -
Exposure parameter 1 Used to determine Eolian_constant_3 that defines the shape of eolian erosion/deposition curve as related to the exposure index -
Exposure parameter 2 Used to determine Eolian_constant_3 that defines the shape of eolian erosion/deposition curve as related to the exposure index -
Distance decay factor In calculation of exposure index, this determines how rapidly the weighting of surrounding cells decays with distance. Also used for accretion/ablation modeling -
Weighting calculation distance The maximum distance (in number of cells) that elevations are used to calculate the exposure index. Also used for accretion/ablation modeling -
Parameter Description Unit
Number of lava sources Number of sources (vents, volcanoes) that are present on the surface (Multiple sources not yet supported) -
Lava source I coordinate I coordinate column of lava source in DEM -
Lava source J coordinate J coordinate column of lava source in DEM -
Lava event probability Probability of a lava flow event during a single iteration -
Minimum lava flow slope Minimum gradient for laval flow at the flow source -
Lava flow thickness Assumed thickness of individual lava flow deposits meters
Minimum lava flow thickness Minimum thickness of a lava flow that can flow into adjoining cells meters
New segment interval Number of interations before starting a new lava flow source -
Source segment interval Number of iterations between changeover to a different lava flow source -
Eruption stop probability Probability per iteration, that the existing flow will solidify and stop being active. If this happens, a new flow starts at the source -
No flow probability The lower limit of probability for a cell to be a source for a new flow segment. If the probability drops below this value then the cell is no longer considered to be a possible flow source -
Lava gradient weight Determines how much the gradient between the edge of a flow and the neighboring point determines the probability of flow in that direction -
Lava duration weight Determine how rapidly a new cell diminishes in probability that it can be the source of a flow into a neighboring cell -
Parameter Description Unit
Impact probability probability of an impact event, per year -
Switch for ejecta wraparound If set to 1 and if the domain is periodic in both X and Y, then ejecta deposition can carry over onto the opposite side (IFOLD, 0 or 1) -
Switch for regolith crater If set to 1, then crater slopes and ejecta are soft (regolith). Otherwise they are considered to be initially bedrock -
Large crater depth scale Based upon scaling of fresh impact craters in Forsberg et al. (2004) and Howard (2007) meters
Large crater depth exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Large crater rim scale Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) meters
Large crater rim exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007)
Transition diameter Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) meters
Small crater depth scale Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) meters
Small crater depth exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Small crater rim scale Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) meters
Small crater rim exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Large crater shape scale Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Large crater shape exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Small crater shape scale Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Small crater shape exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Crater frequency exponent Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Frequency cutoff scaling Based on scaling of fresh impact craters in Forsberg et al (2004) and Howard (2007) -
Smallest possible crater Depends on the scale of the simulation, generally in meters, Generally keep at about (4*cellscale) -
Smallest modelled crater Depends on the scale of the simulation, generally in meters, Generally keep at about (4*cellscale) -
Largest modeled crater Depends on the scale of the simulation, generally in meters. Generally keep at about (cell_scale*ncols/2) -
Ejecta thickness variability - -
Noisesd Standard deviation of random noise used for ejecta thickness -
Inheritance parameter - -
Parameter Description Unit
Switch for variable ocean elevation If set to 1, ocean elevation varies in time and the times and corresponding sea levels are read in from the file "oceanlevels.dat" -
Parameter Description Unit
Accretion rate Rate of non-fluvial and non-eolian surface accretion and degradation -
Switch for exposure dependent creep If set to 1, the mass wasting creep rate depends on exposure -
Switch for solar radiation If set to 1, sublimation due to reflected radiation is modeled -
Switch for top exposure - -
Switch for inverse exposure - -
Switch for exposure smoothing - -
Radiation dust factor Determines the relative amount of re-emitted thermal radiation from dust-covered surface relative to bedrock surface -
Radiation threshold convexity Critical value of exposure index for redeposition of sublimated ice -
Radiation deposit rate Rate scaling for redeposition of sublimated ice on less-exposed surfaces -

## Uses ports

This will be something that the CSDMS facility will add

## Provides ports

This will be something that the CSDMS facility will add

## Main equations

• Local vertical rate of land surface elevation change
 $\displaystyle{ {\frac{\partial y}{\partial t}} - U = {\frac{\partial z}{\partial t}} = - \nabla \ast q }$ (1)
• Weathering, Rain Slash, and Mass Movement

1) Potential erosion or deposition due to rain splash and regolith mass movement

 $\displaystyle{ {\frac{\partial z}{\partial t}}|_{m} = -\nabla \ast q_{m} }$ (2)

2) The rate of movement

 $\displaystyle{ q_{m} = [K_{s} \wp \left (S\right ) + K_{f}{\frac{1}{\left (1 - K_{x}|S|^a\right )}}] s }$ (3)
 $\displaystyle{ \wp \left (S\right ) = sin \theta }$ (4)

A linear dependency

 $\displaystyle{ \wp \left (S\right ) = tan \theta = S }$ (5)
• Fluvial Processes

1) Potential channel deposition or erosion

 $\displaystyle{ {\frac{\partial z}{\partial t}}|_{c} = - \nabla \ast q_{s} = - {\frac{\partial q_{s}}{\partial x}} }$ (6)

2) Nonalluvial channels (A) Detachment capacity

 $\displaystyle{ -{\frac{\partial z}{\partial t}}|_{c} = {\frac{\partial q_{s}}{\partial x}} = C_{d} = K_{t} \left (\tau - \tau_{c}\right ) }$ (7)

(B) equations for steady, uniform flow

 $\displaystyle{ \tau = \gamma RS }$ (8)
 $\displaystyle{ V = K_{n}R^\left ({\frac{2}{3}}\right ) S^\left ({\frac{1}{2}}\right ) / N }$ (9)
 $\displaystyle{ Q_{w} = K_{p} RWV }$ (10)

(C) Simple power law equations

 $\displaystyle{ Q_{w} = K_{a} A^e }$ (11)
 $\displaystyle{ W = K_{w}Q_{w}^b = K_{w} K_{a}^b A^ \left (be\right ) }$ (12)

(D)

 $\displaystyle{ {\frac{\partial z}{\partial t}}|_{c} = - K_{t} \left (K_{z} A^g S^h - \tau_{c}\right ) }$ (13)
 $\displaystyle{ g = 0.6 e \left (1 - b \right ) }$ (14)
 $\displaystyle{ h = 0.7 }$ (15)
 $\displaystyle{ K_{z} = \gamma \{{\frac{N K_{a}^\left (1 - b\right )}{K_{p} K_{w} K_{n}}}^ \left ({\frac{3}{5}}\right ) }$ (16)

(E) Channel elevation change (a) The net elevation change of the channel when sufficiently large amount of regolith is delivered to the cell (the channel never erodes bedrock)

 $\displaystyle{ {\frac{\partial z}{\partial t}} = {\frac{\partial z}{\partial t}}|_{m} + F \xi {\frac{\partial z}{\partial t}}|_{c} }$ (17)
 $\displaystyle{ \xi = W / \delta }$ (18)

(b) The channel erosion when the channel is capable of eroding all regolith delivered to it

 $\displaystyle{ \epsilon = {\frac{\partial z}{\partial t}}|_{m} \delta^2 - {\frac{\partial z}{\partial t}}\delta \left ( \delta - W \right ) }$ (19)

The fraction of time the channel spends eroding regolith

 $\displaystyle{ \eta = {\frac{\epsilon}{F \xi \delta^2 {\frac{\partial z}{\partial t}}|_{c}}} }$ (20)

The net erosion into bedrock for the channel

 $\displaystyle{ {\frac{\partial z}{\partial t}} = \left ( 1 - \eta \right ) {\frac{\partial z}{\partial c}}|_{c} = \beta \{{\frac{\partial z}{\partial t}}|_{m} + F \xi {\frac{\partial z}{\partial t}}|_{c} \} }$ (21)
 $\displaystyle{ \beta = {\frac{1}{F \xi + 1 - \xi}} }$ (22)

(c) Channel erosion when the volume of regolith delivered to the cell per unit time is negative (the removal of regolith by channel erosion pertains only to the volume surrounding the stream and below the level of removal by slope erosion)

 $\displaystyle{ \xi = \{{\frac{\partial z}{\partial t}}|_{m} - {\frac{\partial z}{\partial t}}\} \delta \left (\delta - W \right ) }$ (23)

The net erosion into bedrock

 $\displaystyle{ \{{\frac{\partial z}{\partial t}}|_{m} - {\frac{\partial z}{\partial t}}\} = \left ( 1 - \eta \right ) {\frac{\partial z}{\partial t}}|_{c} }$ (24)

3) Alluvial channels (A) Potential rate of fluvial erosion

 $\displaystyle{ {\frac{\partial z}{\partial t}}|_{c} = - \nabla \ast q_{sb} }$ (25)
 $\displaystyle{ \phi = \Im \left ({\frac{1}{\Psi}}\right ) = {\frac{q_{sb}}{\omega d \left ( 1 - \mu \right )}} }$ (26)
 $\displaystyle{ {\frac{1}{\Psi}} = {\frac{\tau}{\left ( \gamma_{s} - \gamma \right ) d }} }$ (27)

 $\displaystyle{ \phi = K_{e} \{{\frac{1}{\Psi}} - {\frac{1}{\Psi_{c}}}\} ^p }$ (28)

(C) Total bed sediment discharge

 $\displaystyle{ Q_{sb} = K_{q} A^ \left (eb \right ) [K_{v} A^ \left (0.6 e \left (1 - b \right ) \right ) S^ \left (0.7\right ) - 1/ \Psi_{c}]^ p }$ (29)
 $\displaystyle{ K_{q} = {\frac{K_{e} \omega \left (1 - \mu \right ) K_{w} K_{a}^b}{\left ( \gamma_{s} - \gamma \right )^p d^ \left ( p - 1 \right )}} }$ (30)
 $\displaystyle{ K_{v} = \gamma \{{\frac{N K_{a}^ \left (1-b\right) }{K_{w} K_{n} K_{p}}}\} ^ \left (0.6\right ) }$ (31)

(D) The net erosion of channel

 $\displaystyle{ {\frac{\partial z}{\partial t}}|_{c} = {\frac{\left (Q_{sbi} - Q_{sbo}\right )}{LW}} }$ (32)

(4) Overall volume V of the alluvial channel between the nonalluvial upstream node at elevation Z0 and the downstream unmodified node with elevation En+1

 $\displaystyle{ V = \delta^2 [\left (n + {\frac{1}{2}} \right ) Z_{0} - n S_{0} l_{0} - \Sigma_{i = 1} ^ \left ( n-1\right ) S_{i} l_{i} + {\frac{1}{2}}Z_{n+1} ] }$ (33)

The net deposited volume of sediment for an individual nonalluvial upstream contribution during the time step

 $\displaystyle{ Q_{sj} \Delta t = v - V }$ (34)

Upstream gradient of the alluvial channel node

 $\displaystyle{ \zeta_{0} = S_{0} - [\Sigma_{i=1} ^ \left (n-1 \right ) \left ( n - i \right ) \left ( S_{i} - \zeta_{i} \right ) l_{i} + \left (Q_{sj} \Delta t \right ) / \delta^2] / [n l_{0}] }$ (35)

Gradient of the alluvial channel node n

 $\displaystyle{ \zeta_{n} = [\left (Z_{0} - Z_{n+1} \right ) - \zeta_{0}l_{0} - \Sigma_{i=1}^ \left (n-1\right ) \zeta_{i}l_{i}] }$ (36)

The new elevations of the alluvial channel point i when all of the bed sediment supplied form upstream is deposited

 $\displaystyle{ J_{i} = Z_{0} - \zeta_{0}l_{0} - \Sigma_{k=1}^ \left (i-1 \right ) \zeta_{k}l_{k} }$ (37)

The new elevations of the alluvial channel point i when location n is a fixed boundary point with elevation Zf

 $\displaystyle{ J_{i} = Z_{f} + \Sigma_{k = n-1}^ \left (i+1\right ) \zeta_{k}l_{k} }$ (38)

## Notes

General features and assumptions Within each square matrix cell of dimension δ by δ both slope and channel processes occur. The width W of the active channel is specified to be less than δ, even for alluvial pediments and alluvial fans, corresponding to the observation that flow on such surfaces is generally channelized and confined to a few active channels at any time. Both mass wasting and fluvial transport/erosion occur in each cell. In this model, it is assumed that runoff erosion occurs only in channelized flow.On natural headwater slopes there might be several ephemeral rills within the area represented by one simulation grid cell, but it is assumed that all runoff becomes concentrated into a single permanent or ephemeral channel running the length of the cell with a gradient equal to the overall slope gradient and that all fluvial erosion within the cell occurs in this channel. Erosion by runoff will be somewhat exaggerated in the present model because of the scale efficiency associated with a single channel rather than multiple small channels. Nonetheless, channel erosion becomes negligible compared with mass wasting on matrix cells located near divides.

Weathering, Rain Splash, and Mass Movement

Weathering is not explicitly modeled, but is assumed to keep pace with erosion by rain splash, mass wasting, or fluvial processes such that everywhere except in some channels there is a weathered, cohesive regolith. Creep or slope failure is assumed to occur at a sufficiently shallow depth with regard to the thickness of the regolith that its thickness is not a limiting factor. The rate of movement calculated in equation (3) is expressed by two additive terms, one for creep and/or rain splash diffusion and one for near-failure conditions.

Fluvial Processes

The model incorporates both detachment-limited erosion on slopes and steep headwater channels as well as sediment transport in alluvial channels. The channels can be classified as nonalluvial and alluvial channels. The nonalluvial channel is defined as one in which the bed load sediment flux is less than a capacity load, while for alluvial channels the potential rate of fluvial erosion equals the spatial divergence of the volumetric unit bed sediment transport rate.

In the nonalluvial channels, the detachment capacity is assumed to be proportional to the shear stress exerted on the bed and banks by a dominant discharge, shown as equation (7). Shear stress can be related to channel gradient and drainage area through the use of equations of steady, uniform flow, as equation (8,9,10). However, the model also assume simple power law equations of hydraulic geometry for dominant discharge and channel width (equation (11,12) ) for runoff on slopes as well as in channels, which implies that little regolith or depression storage of precipitation occurs during erosion events, since runoff production is assumed to be areally uniform. On badland slopes this is probably a reasonable approximation. For landscapes with appreciable infiltration capacity, particularly where vegetated, the assumption is made that most runoff erosion on slopes occurs during infrequent, very intense rainfall events due to saturateion overland flow or shallow interflow through large macropores.

Channel maybe eroding both regolith material which is delivered into the channel by mass wasting as well as uneroded bedrock. Regolith material is assumed to be more erodible than bedrock by a factor F, where F >= 1. In landscapes with a shallow weathering regolith, such as in badlands and steep mountain slopes, bedrock is commonly exposed in headwater rills and channels, and a large ratio of bedrock to regolith erodibility (F >> 1) would be appropriate. However, in landscapes in deeply saprolitized bedrock, till, or uncemented alluvium, the parent material and the surface soil involved in mass wasting may have nearly equivalent erodibility (F ~ 1).

In calculating net erosion in matrix cells containing nonalluvial channels, during each time step, the channel is assumed to erode first regolith delivered to the channel by slope erosion and then, if it is capable of eroding all regolith, bedrock. The fraction of time during each time step that the channel is eroding regolith is η (0 <= η <= 1). There are three cases, the simplest of which occurs when sufficiently large amount of regolith is delivered to the cell that the channel never erodes bedrock (η = 1). The second case is that when the channel is capable of eroding all regolith delivered to it, and then the value of η must be determined. The final case explains the condition that the removal of regolith by channel erosion pertains only to the volume surrounding the stream and below the level of removal by slope erosion.

In alluvial channels, the bed load or total load formula is expressed as equation (28). For sand bed channels with high effective discharges, the exponent p in the formula is generally about 3, and 1/ψc is a negligible correction.

As to the sediment flux through channels, several adjacent cells may contribute to the sediment flux to a cell, but the outflux occurs to the adjacent cell with the steepest download gradient. In the case of a depression there is no outflux. Influxes may originate from alluvial streams debauching into the cell from adjoining cells or from nonalluvial streams and rills. In the latter case the incoming total sediment flux equals the total volume of upstream nonalluvial channel erosion less the proportion of load carried as wash load, corresponding to the assumption that downstream routing of sediment through nonalluvial channels is very rapid compared to the iteration time step.

Although alluvial transport is assumed to occur in well-defined channels whose width is lass than the simulation cell width, over long periods of time the channel is assumed to erode or deposit over the entire area of the cell. So that the divisor LW in equation (32) becomes the area of the cell δ2.

In the model implementation, a primary concern is selection of a temporal increment that assures numerical stability. Each of the processes has a characteristic maximum time step, and the overall model must utilize the minimum of these. However, the maximum time step for sediment transport divergence erosion is less by a factor of 103 to 104 than that for slope or bedrock channel erosion. This occurs primarily for channels with large drainage area and low gradient. Because the exponent of gradient in equation (29, 30, 31) is about 2, small changes in gradient create large variations in transport rate so that time steps must be very small to assure that numerical errors are not amplified. An alternative approach for alluvial channels with large drainage areas are employed in the mode, which permits simulation of alluvial channel evolution using workstations rather than supercomputers. The timescale for alluvial channel regarding due to changes in supply of sediment or water from upstream or due to local uplift base level changes is assumed to be short compared to that for overall evolution of basin relief. The short timescales of alluvial channel regarding is manifest in many drainage basins by the often multifold deposition and dissection of terraces in response to changes in hydraulic regime over time periods during which overall basin relief has changed very little. The approach is to utilize the longer time step which is appropriate for slope and bedrock channel erosion, thus ignoring the type of short-term fluctuations responsible for terraces.

At the head of each section of alluvial channel a discharge of sediment Qsj is delivered by upstream channel erosion, making a volumetric contribution QsjΔt during each time step; the subscript j is used because each channel or slope draining into the alluvial channel contributes individually to the overall sediment discharge. This volumetric contribution is routed downstream, in general, being deposited within n segments (elements) of the alluvial channel system. If Zi is the mean elevation of a given alluvial segment, then the overall volume of the alluvial channel between the nonalluvial upstream node at elevation Z0 and the downstream unmodified node with elevation En+1 is calculated as the equation (33). Equation (34,35,36) are used to calculate the volume contribution and gradients iteratively, increasing n until one of the following occurs: (1) the gradient of the final segment is less than or equal to the steady state alluvial gradient for the local drainage area and total sediment discharge, or (2) location n is a fixed boundary point with elevation. For the former case, the new elevation of the alluvial channel points are calculated with equation (37), and for the latter case, the channel is graded upstream from the fixed boundary point, with the assumption that undeposited sediment load is carried from the system, using equation (38).

Special treatment is required for deltaic deposition into a water-filled depression or sea. Submerged alluvial bed sediment deposits are assumed to be foreset beds with an angle of repose gradient ζr. The routing scheme takes into consideration that a channel segment with a subareal upstream node and a submerged downstream node will have a portion of the segment having gradient ζi and part with gradient ζr. Another special circumstance occurs for the first case when one or more of the downstream locations is nonalluvial. If one of these locations is nonalluvial and the predicted new elevation is less than the existing elevation, then the alluvial channel is backgraded from the existing elevation using a formula similar to equation (38), and the bed sediment that is not deposited is routed further downstream through the bedrock section. Since many nonalluvial slopes and channel may each contribute sediment during each time step, the calculation is repeated for each of these sources, with each previously calculated new elevation Ji, now becoming the initial elevation Zi for the new source. Thus alluvial channels may undergo several subiterations of deposition during each time step. Finally, the upstream location at Z0 is converted to an alluvial channel if the calculated new gradient ζ0 is less than the steady state alluvial gradient calculated using equation (29,30,31).

Transitions

All cells in the simulation matrix have a contribution to erosion or deposition from either an alluvial or a nonalluvial channel, but the two types of channels are mutually exclusive. Thus the model must provide criteria for temporal and spatial transitions. If the channel in a cell is currently nonalluvial, the sediment arriving from upstream and that locally eroded by the channel is considered to be transported without deposition and to be routed instantaneously (relative to the simulation timescale) through the nonalluvial network. During each iteration the total potential bedload transport capacity of the channel in each cell is calculated, and if the actual transport rate exceeds the potential rate, then the cell is converted to an alluvial channel. When such a transition occurs the base of the alluvial deposit is set equal to the local elevation. So long as erosion due to sediment divergence in an alluvial channel during an iteration is less than the thickness of the alluvium, the channel remains alluvial. However, the alluvial channels in nature and in the model may also occur under conditions of slow downcutting (such that the base of the alluvial deposit is constantly reset during each iteration to the local elevation). The maximum amount of erosion that an alluvial channel on shallow alluvium over bedrock can accomplish is considered to equal the rate for nonalluvial channels into bedrock. If the erosion due to sediment divergence or the routing procedure exceeds this value within a channel on shallow alluvium, then the channel is converted to a nonalluvial channel.

## Examples

An example run with input parameters, BLD files, as well as a figure / movie of the output

Follow the next steps to include images / movies of simulations: