Model help:MARSSIM
MARSSIM
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
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
<math> {\frac{\partial y}{\partial t}} - U = {\frac{\partial z}{\partial t}} = - \nabla \ast q </math> (1)
- Weathering, Rain Slash, and Mass Movement
1) Potential erosion or deposition due to rain splash and regolith mass movement
<math> {\frac{\partial z}{\partial t}}|_{m} = -\nabla \ast q_{m} </math> (2)
2) The rate of movement
<math> q_{m} = [K_{s} \wp \left (S\right ) + K_{f}{\frac{1}{\left (1 - K_{x}|S|^a\right )}}] s </math> (3)
<math> \wp \left (S\right ) = sin \theta </math> (4)
A linear dependency
<math> \wp \left (S\right ) = tan \theta = S </math> (5)
- Fluvial Processes
1) Potential channel deposition or erosion
<math> {\frac{\partial z}{\partial t}}|_{c} = - \nabla \ast q_{s} = - {\frac{\partial q_{s}}{\partial x}} </math> (6)
2) Nonalluvial channels (A) Detachment capacity
<math> -{\frac{\partial z}{\partial t}}|_{c} = {\frac{\partial q_{s}}{\partial x}} = C_{d} = K_{t} \left (\tau - \tau_{c}\right ) </math> (7)
(B) equations for steady, uniform flow
<math> \tau = \gamma RS </math> (8)
<math> V = K_{n}R^\left ({\frac{2}{3}}\right ) S^\left ({\frac{1}{2}}\right ) / N </math> (9)
<math> Q_{w} = K_{p} RWV </math> (10)
(C) Simple power law equations
<math> Q_{w} = K_{a} A^e </math> (11)
<math> W = K_{w}Q_{w}^b = K_{w} K_{a}^b A^ \left (be\right ) </math> (12)
(D)
<math> {\frac{\partial z}{\partial t}}|_{c} = - K_{t} \left (K_{z} A^g S^h - \tau_{c}\right ) </math> (13)
<math> g = 0.6 e \left (1 - b \right ) </math> (14)
<math> h = 0.7 </math> (15)
<math> K_{z} = \gamma \{{\frac{N K_{a}^\left (1 - b\right )}{K_{p} K_{w} K_{n}}}^ \left ({\frac{3}{5}}\right ) </math> (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)
<math> {\frac{\partial z}{\partial t}} = {\frac{\partial z}{\partial t}}|_{m} + F \xi {\frac{\partial z}{\partial t}}|_{c} </math> (17)
<math> \xi = W / \delta </math> (18)
(b) The channel erosion when the channel is capable of eroding all regolith delivered to it
<math> \epsilon = {\frac{\partial z}{\partial t}}|_{m} \delta^2 - {\frac{\partial z}{\partial t}}\delta \left ( \delta - W \right ) </math> (19)
The fraction of time the channel spends eroding regolith
<math> \eta = {\frac{\epsilon}{F \xi \delta^2 {\frac{\partial z}{\partial t}}|_{c}}} </math> (20)
The net erosion into bedrock for the channel
<math> {\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} \} </math> (21)
<math> \beta = {\frac{1}{F \xi + 1 - \xi}} </math> (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)
<math> \xi = \{{\frac{\partial z}{\partial t}}|_{m} - {\frac{\partial z}{\partial t}}\} \delta \left (\delta - W \right ) </math> (23)
The net erosion into bedrock
<math> \{{\frac{\partial z}{\partial t}}|_{m} - {\frac{\partial z}{\partial t}}\} = \left ( 1 - \eta \right ) {\frac{\partial z}{\partial t}}|_{c} </math> (24)
3) Alluvial channels (A) Potential rate of fluvial erosion
<math> {\frac{\partial z}{\partial t}}|_{c} = - \nabla \ast q_{sb} </math> (25)
<math> \phi = \Im \left ({\frac{1}{\Psi}}\right ) = {\frac{q_{sb}}{\omega d \left ( 1 - \mu \right )}} </math> (26)
<math> {\frac{1}{\Psi}} = {\frac{\tau}{\left ( \gamma_{s} - \gamma \right ) d }} </math> (27)
(B) Bed load or total bedload formulas
<math> \phi = K_{e} \{{\frac{1}{\Psi}} - {\frac{1}{\Psi_{c}}}\} ^p </math> (28)
(C) Total bed sediment discharge
<math> 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 </math> (29)
<math> 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 )}} </math> (30)
<math> K_{v} = \gamma \{{\frac{N K_{a}^ \left (1-b\right) }{K_{w} K_{n} K_{p}}}\} ^ \left (0.6\right ) </math> (31)
(D) The net erosion of channel
<math> {\frac{\partial z}{\partial t}}|_{c} = {\frac{\left (Q_{sbi} - Q_{sbo}\right )}{LW}} </math> (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
<math> 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} ] </math> (33)
The net deposited volume of sediment for an individual nonalluvial upstream contribution during the time step
<math> Q_{sj} \Delta t = v - V </math> (34)
Upstream gradient of the alluvial channel node
<math> \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}] </math> (35)
Gradient of the alluvial channel node n
<math> \zeta_{n} = [\left (Z_{0} - Z_{n+1} \right ) - \zeta_{0}l_{0} - \Sigma_{i=1}^ \left (n-1\right ) \zeta_{i}l_{i}] </math> (36)
The new elevations of the alluvial channel point i when all of the bed sediment supplied form upstream is deposited
<math> J_{i} = Z_{0} - \zeta_{0}l_{0} - \Sigma_{k=1}^ \left (i-1 \right ) \zeta_{k}l_{k} </math> (37)
The new elevations of the alluvial channel point i when location n is a fixed boundary point with elevation Zf
<math> J_{i} = Z_{f} + \Sigma_{k = n-1}^ \left (i+1\right ) \zeta_{k}l_{k} </math> (38)
| Symbol | Description | Unit |
|---|---|---|
| U | local tectonic uplift rate | - |
| q | spatial divergence of the vector of eroded material flux | - |
| y | local land surface elevation | - |
| z | land surface elevation caused by fluvial erosion relative to a bedrock-fixed reference frame | - |
| m | Potential erosion or deposition due to rain splash and regolith mass movement | - |
| qm | spatial divergence of the vector rate of movement | - |
| \wp (S) | an increasing function of slope gradient | - |
| s | the unit vector in the direction of S | - |
| S | local slope gradient | - |
| Ks, Kx, Kf | constants assumed to be spatially and temporally invariant | - |
| a | exponent assumed to be spatially and temporally invariant | - |
| θ | slope angle | - |
| x | downstream direction of coordinate | - |
| c | potential channel deposition or erosion | - |
| qs | volumetric sediment transport rate per unit channel width | - |
| qsw | wash load (a part of sediment discharge, which is assumed to never be redeposited except in depressions) | - |
| qsb | bed load (including bed load and suspended load), which is carried in capacity amounts if an alluvial bed is present | - |
| τ | shear stress | - |
| τc | critical shear stress | - |
| C | actual detachment rate | - |
| Cd | intrinsic detachment capacity | - |
| γ | unit weight of water | - |
| R | the hydraulic radius | - |
| N | Manning's resistance coefficient | - |
| Kp | a form factor close to unity | - |
| Qw | hydraulic geometry for dominant discharge (Assumed to be equal for both bed erosion and sediment transport) | - |
| W | channel width | - |
| Kt | this constant includes both effects of substrate erodibility as well as magnitude of the dominant discharge | - |
| F | radio of bedrock to regolith erodibility | - |
| η | the fraction of time during each time step that the channel is eroding regolith,0=<η<=1 | - |
| δ | the length of channel in the cell | - |
| ε | channel erosion | |
| Φ, 1/Ψ | dimensionless parameters for bed load and total load sediment transport in alluvial channels | - |
| ω | fall velocity of the sediment grains | - |
| d | sediment grain size | - |
| μ | alluvium porosity | - |
| γs | unit weight of sediment grains | - |
| 1/Ψc | the threshold for transport | - |
| p | exponent in total sediment load formula, about 3 for sand bed channels | - |
| Qsbi | total sediment influx into a given cell | |
| Qsbo | outflow from a given cell | |
| L | length of the channel through the cell | - |
| Zi | mean elevation of a given alluvial segment | L |
| V | overall volume of the alluvial channel between the nonalluvial upstream node at elevation Z0 and the downstream unmodified node with elevation En+1 | - |
| v | sediment volume at the end of the timestep | - |
| Ji | final elevations | L |
| ζi | final gradients | - |
| n | number of the alluvial channel nodes | - |
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:
- Upload file: https://csdms.colorado.edu/wiki/Special:Upload
- Create link to the file on your page: [[Image:<file name>]].
See also: Help:Images or Help:Movies
Developer(s)
References
Howard, A.D., 1994. A detachment-limited model of drainage basin evolution. Water Resources Research, 30,7: 2261~2285
