# Surface processes model¶

The whole landscape evolution can basically be represented in a single equation:

All landscape must obey this fundamental statement about sediment transport. In this equation \(\frac{\partial z}{\partial t}\) is the change in surface elevation, \(U\) is the uplift rate, \(E\) is the erosion rate and \(\nabla \cdot q_s\) is the sediment flux divergence.

The erosion rate \(E\) corresponds to sediment production from weathering and bedrock erosion by glacier, wind, water. The sediment fluxes is transported by hillslope and fluvial transport processes.

Note

Many geomorphological applications have demonstrated the usefulness of these models, whose predictions help researchers test simple to complex hypotheses on the nature of landscape evolution. Surface processes models (SPMs) also provide connection between small scale, measurable processes and their long-term geomorphic implications.

## Preamble¶

The roots of landscape evolution theory can be found in the pioneering work of Gilbert (1877), who proposed a set of hypotheses to relate various landforms to the mechanisms of weathering, erosion and sediment transport. The first quantitative models appeared later in the 1960s (e.g., Culling, 1960; Scheidegger, 1961; Ahnert, 1970; Kirkby, 1971). These models formalise the concepts of Gilbert (1877) to the development of hillslope profiles. A few years later, these models were extended to two dimensions, although still focused on hillslope morphology.

During the last two decades, as computers continued to get faster, a number of sophisticated numerical SPMs have been developed, mainly focusing on watershed and mountain belt evolution. Both hillslope and fluvial processes are involved in these models, which differ from each other by the parameterisation of these processes and their numerical resolution.

## Continuity of mass¶

In the case where there is no distinction between a regolith layer and the bedrock underneath, the mass continuity equations for a column of soil or rock is expressed as:

where the elevation \(z\) (m) is measured vertically, \(q_s\) is the total downhill soil flux, \(\nabla \cdot\) is the spatial divergence operator and \(U\) (m/yr) is a source term that can either represent the rate of incision of channel streams at the hillslope boundaries or uniform uplift.

## Sediment transport¶

To describe the rates of sediment erosion/transport/deposition, several approaches have been proposed. In its simplest formalism a **detachment-limited equation** is often used.

### Detachment-limited model¶

The soil transport rate per unit width by flowing water, \(q_r\), is modelled as a power function of topographic gradient \(\nabla z\) and surface water discharge per unit width \(q_w\) (m2/yr):

This detachment-limited incision rate, which is calculated as a power law function of fluvial discharge only applies where channel slope is positive. This brings the following relatioship:

This expression corresponds to a simplified form of the usual expression of sediment transport by water flow, in which the transport rate is assumed to be equal to the local carrying capacity, which is itself a function of boundary shear stress or stream power per unit width. We consider additionally no threshold for particle entrainment. Generally, the exponents m and n have values between 1 and 2.

### Hillslope processes¶

In its most simple form, the parameterisation of hillslope transport is based on a **linear dependence to the topographic gradient**. This linear law has in fact been used to represent a variety of transport processes such as **creep**, **biogenic activity** or **rain splash**.

Downslope simple creep is commonly regarded as operating in a shallow superficial layer and is defined as:

Note that because of the multi-process parameterisation of soil transport, the coefficient \(\kappa_d\) is also scale-dependent, like the \(\kappa\)-scale parameters of the other stream power law defined above.

## Incision laws overview¶

Important

Several formulations of river incision have been proposed to account for long term evolution of fluvial system. These formulations describe different erosional behaviours ranging from **detachment-limited incision**, governed by bed resistance to erosion, to **transport-limited incision**, governed by flow capacity to transport sediment available on the bed.

As we already discussed, mathematical representation of erosion processes in these formulations is often assumed to follow a stream power law. These relatively simple approaches have two main advantages. First, they have been shown to approximate the first order kinematics of landscape evolution across geologically relevant timescales (> \(10^4\) years). Second, neither the details of long term catchment hydrology nor the complexity of sediment mobilisation dynamics are required. However, other formulations are sometimes necessary when addressing specific aspects of landscape evolution.

### Detachment-limited¶

The simplest law to simulate fluvial incision is based on the detachment-limited stream power law (option 1, in the above figure), in which erosion rate depends on drainage area \(A\), net precipitation \(P\) and local slope \(S\) and takes the form:

\(\kappa_d\) is a dimensional coefficient describing the erodibility of the channel bed as a function of rock strength, bed roughness and climate, \(l\), \(m\) and \(n\) are dimensionless positive constants.

Default formulation assumes \(l = 0\), \(m = 0.5\) and \(n = 1\). The precipitation exponent \(l\) allows for representation of climate-dependent chemical weathering of river bed across non-uniform rainfall. In this model sediment deposition occurs solely in topographically closed depression and offshore.

### Transport-limited¶

Here, volumetric sediment transport capacity (\(Q_t\)) is defined using a power law function of unit stream power:

where \(\kappa_t\) is a dimensional coefficient describing the transportability of channel sediment and \(m_t\) and \(n_t\) are dimensionless positive constants. In this equation, the threshold of motion (the critical shear stress) is assumed to be negligible.

An additional term is now introduced in the stream power model:

with \(f(Q_s)\) representing a variety of plausible models for the dependence of river incision rate on sediment flux \(Q_s\). In the standard detachment-limited, \(f(Q_s)\) is equal to unity, which corresponds to cases where \(Qs << Qt\). All sediment is dispersed downstream and the incision limiting factor is bedrock erodibility.

Addition of the transport-limited function results in the fact that, where sediment flux equals or exceeds transport capacity (\(Q_s/Q_t \ge 1\)) the system becomes transport-limited and depositional if \(Qs/Qt > 1\). In this model the time-evolving distribution of erosion and sedimentation, is affected by the distribution of detachment-limited and transport-limited reaches, which is controlled by the respective values of \(\kappa_d\) and \(\kappa_t\).

The transition from one behaviour to the other can be treated either abruptly, progressively, through the use of one of the following formulations:

**Linear decline**: This model belongs to the undercapacity family of models: it assumes that stream incision represents the expenditure of the energy in excess of that needed to transport the bypassing sediment load. Stream incision potential decreases linearly from a maximum where sediment flux is negligible, to zero where sediment flux equals transport capacity (option 2).**Almost parabolic**: Both qualitative and experimental observations have shown that sediment flux has a dual role in the river bed incision. First, when sediment flux is low relative to carrying capacity, erosion potential increases with sediment flux (tool effect: bedrock abrasion and plucking). Then, with increased sediment flux, erosion is inhibited (cover effect: sediments protect the bed from impacts by saltating particles) (option 3).**Dynamic cover**: Typically gravel-river beds have an armoured layer of coarse grains on the surface, which acts to protect finer particles underneath from erosion. To account for sediment and spatial heterogeneity in the armouring of the river bed, Turowski et al. proposed a modified form of the ‘almost parabolic’ model that better estimates the original Sklar and Dietrich experiments (option 4).

## Step-by-step approach to landscape evolution model¶

### Step 1: Flow directions¶

Important

Landscape evolution applications generally require computing the **drainage network** of a terrain, consisting of the **flow direction** and **flow accumulation**. Intuitively, they are the path that water flows through the terrain and the amount of water that flows into each terrain cell supposing that each cell receives a rain drop

The drainage network of a terrain delineates the path that water flows through the terrain (the flow direction) and the amount of water that flows into each terrain cell (the flow accumulation). The flow direction problem is to assign flow directions to all cells in the terrain such that the following three conditions are fulfilled:

- Every cell has at least one flow direction;
- No cyclic flow paths exist; and
- Every cell in the terrain has a flow path to the edge of the terrain.

The flow direction can be modelled considering **single flow direction** (**SFD**) or **multiple flow directions** (**MFD**). In SFD, each terrain cell is assigned a direction towards the **steepest downslope neighbour**, while in MFD, each cell has directions to **all downslope neighbours**. The use of SFD or MFD is essentially a modeling choice since the computational complexity of the flow routing problem is the same in both models.

### Step 2: Pit filling¶

Important

The major challenge in the process is the flow routing in **local minimum** and **flat areas**. A local minimum is a cell with no downslope neighbour and a flat area is a set of adjacent cells with a same elevation.

A neighbour cell of c is called a **downslope neighbour** if it has a strictly lower elevation than c. A cell in a flat area that has a downslope neighbour is called a **spill-point**. Also, a flat area can be classified as a plateau or a sink where the plateau has a spill point and a sink doesn’t. Intuitively, water will accumulate in a sink until it fills up and water flows out of it, while in the plateau the water should flow towards spill points.

Usually, most drainage network computation methods use a **preprocessing step to remove the sinks and the flat areas**. Initially, the elevation of the cells belonging to a sink are increased to transform it into a plateau. Next, the flow directions on the plateau are assigned to ensure that there is a path from each cell to the nearest spill point.

Pit-filling exercise

Starting from the **source** cell and considering a single flow direction (SFD) approach answer these questions:

- There is a pit in this DEM, identify it, which value is required to fill it and allow the flow to keep moving downslope?
- Where is the water from the source entering the lake?

### Step 3: Flow accumulation and erosion¶

After obtaining the flow direction, the next step consists in computing the **flow accumulation** in each terrain cell, that is, the amount of water flowing to each cell supposing that all cells receive a drop of water and this water follows the direction obtained in the previous step.

Once the flow accumulation has been computed for a particular topography, the erosion is then estimated using one of the incision laws defined in the previous section and requires at least the estimation of the slope based on the flow direction. The erosion values are finally used to change the topography elevations and the model moves forward in time.

At the next iteration, steps 1 to 3 are applied on the new elevation grid allowing to simulate landscape evolution over time.