Quoted from: Coulthard, T. J. , M. G. Macklin , and M. J. Kirkby . "A cellular model of Holocene upland river basin and alluvial fan evolution." Earth Surface Processes & Landforms 27(2010). https://onlinelibrary.wiley.com/doi/epdf/10.1002/esp.318 

Channel flow 

Within a catchment, water flows from the highest to the lowest point, so within these models any precipitation or flow must be routed from the highest grid cell or node, down to the lowest. Computationally, the simplest method is to sort all the grid cells or nodes within the catchment by elevation, and then work down through the list. For example, if we have a simple grid of 100 by 100 cells, 10 000 points have to be sorted,and if we allow the elevations to change through erosion and deposition, these cells have to be reordered for every iteration. Computationally this sorting can be time-consuming, restricting the number of grid cells, and therefore the area and number of model iterations that can be studied. To resolve this, authors have used a technique such as the ‘bucket passing algorithm’ (Braun and Sambridge, 1997) where each node is given a ‘bucket’ or set amount of water and asked to pass it to its lowest neighbour. After this, all nodes that have not received any water are local ‘maxima’ and are placed at the top of a list. This process continues until all nodes have been accounted for, resulting in an ordered list. Whilst being significantly quicker, these solutions are only partial, as they require all water to flow in only one direction – that of the steepest descent. At first glance, this may be considered a reasonable approximation, but rivers and drainage networks are full of features characterized by divergent channel flow, for example, mid-channel bars, alluvial fans and deltas. AsHoward (1994) acknowledges when simulating fan and escarpment development: ‘owing to the limitation offlow in one of the eight directions at each cell, the fan assumes a prismatic shape during its growth’. In order to fully capture the dynamics of fluvial behaviour a multiple flow algorithm that permits flow in all possible directions is required.

There are two notable exceptions to these single-flow-direction models. First, Moglen and Bras (1995), use a multiple-flow algorithm, but only on a small (40 by 40) grid. Secondly, Murray and Paola (1994, 1997) developed a novel solution for simulating braided rivers. They used a square grid to represent a reach or section of braided river and assumed that there was flow in one main direction (the direction of the valley floor). Water was then added to cells at the top row of the reach and routed to the three cells immediately downstream in the second row, according to the local bed slope. The model then moved onto this second row and routed flow to the three downstream cells in row 3 and so on, in effect pushing the flow down through the reach. This simulated the flow in a braided reach with both divergent and convergent flow. Unfortunately, there are problems associated with their approach, as it failed to calculate a flow depth and only routes water in one main direction, making it unsuitable for drainage basin simulations or complex channel sections (e.g.meanders).

The model described here uses a ‘scanning’ algorithm that incorporates the simplicity of the Murray and Paola method but allows calculation of flow in all directions. It operates on a square grid and is illustrated schematically in Figure 1. Here there is a small V-shaped valley where the dark dots represent precipitation added to each grid cell. For every iteration, the procedure uses four scans. The first scan (box 1) operates from left to right, pushing water from grid cells to the three lower cells on the right. When the base of the valley is reached the scan proceeds up the right valley wall, but as the cells are all upslope, no flow is routed to them (as per box 2). On the second scan (box 3), the same process is repeated but here working from right to left, leaving all the water from the catchment in the base of the valley. This is routed by a third scan from top to bottom, pushing all the water out of the catchment (box 4). A fourth scan then operates from bottom to top, but for clarity is not shown here. For each scan, the maximum flow through each point is recorded and taken as the discharge for that point. If the three ‘downstream’ cells are all higher than the contributing cell, but the combined water depth and elevation of the contributing cell is higher, water is retained in the contributing cell up to the depth of the obstruction whilst the rest is routed on. This process has the effect of filling hollows that the model may create with water ready for the next iteration. These hollows or sinks frequently have to be removed from DEM data (Goodchild and Mark, 1987; Hutchinson, 1989). This water trapped within the topography also allows this scanning method to simulate flow around a series of complex bends or channel geometries. For instance, in a meander sequence, water may be routed around a first corner, but be trapped by a second. However, in the next iteration, this water is still there, to be released in the next scan, and replaced by more water from upstream, allowing the continuum of flow. This method gives very similar results to those from multiple flow algorithms but takes a fraction of the time. This speed has allowed CAESAR to be applied to large grids with over a million cells but still, use a multiple flow algorithm.

Process representation

Within a drainage basin there are a large number of processes operating over a wide range of spatialand temporal scales, and integrating these within one model can prove difficult. For instance, the complex topography of a braided river channel will require a far higher level of spatial detail than a bare hill slope. Furthermore, parameterization of these processes is especially difficult as their importance can vary over both different time and space scales. For instance, soil creep may be inconsequential when simulating a series of floods, but very important when modelling temperate catchment evolution over thousands of years.

Conversely, the location of channel bars may prove important for the single flood but less so over longer periods. Previous workers have used a variety of techniques to address these problems. Braun and Sambridge(1997) and Tucker et al. (2001) have used an irregular grid with nodes and links to represent the topography, so areas where more detail is required (e.g. river channels) have more nodes and thus more detail than relatively uniform areas (e.g. a hill slope). Other authors (Howard, 1994; Tucker and Slingerland, 1994) have used sub-grid-scale representation of processes, where for example a small river channel (10 m wide) is represented within a larger (50 m) grid cell, and erosion and deposition averaged across the whole cells area. To cope with long-term studies, most models use averages of erosion rates and long time steps that range from years to decades. The danger of such averaging is that the importance of large individual events may be blurred.

Here, a different approach is taken. The previously described routing algorithm allows a large number of smaller grid cells to be used, which reduces the need for sub-grid-scale representation. This lets a river channel occupy a space one or more cells wide, permitting a more detailed simulation of the channel bed topography, including alluvial features such as channel bars, terraces and alluvial fans. CAESAR is optimized by concentrating on cells where fluvial processes dominate, yet periodically checking hill slopes. This procedure involves routing an initial set discharge down the catchment, noting which cells are ‘wet’ and setting a five-cell buffer around this area. Then, for every iteration, fluvial erosion, deposition and local slope processes are calculated within this ‘buffer zone’, and every 100 iterations the entire catchment is checked. If during a flood the discharge exceeds 50 per cent of the set discharge used to calculate the buffer zone, this discharge is doubled and the buffer area recalculated. This allows the drainage network to expand or contract during a flood whilst optimizing the number of cells examined. In addition, a variable time step is used which is related to the volume of erosion and/or deposition. This can reduce the time step to fractions of a second during a flood yet increase it to a maximum of 30 minutes during low-flow situations. This means that the model simulates entire floods with rising and falling stages, instead of averaging over a single storm or several events. CAESAR also counts both iterations and time elapsed during the simulation, which allows certain processes (e.g. soil creep) to be calculated every simulated month, while others that may be more dependent on the erosion activity (e.g. landslides) are determined every set number of iterations.