Sub-catchment Distributed Landscape Evolution Model


Abstract here, ca. 100 words.ff ffff


Introduce the scientific background and the motivation for developing this code.

Ordering and sub-catchment partitioning

Here, the topography is represented with a triangulated irregular network (TIN) data structure. Using a point selection criterion and constrained Delaunay triangulation, dense or high-resolution DEMs obtained from remote sensing are automatically sampled to construct initial triangulated topography. Taking advantage of the Triangle library, it is possible to define higher resolution domains within the entire simulation area.

Various factors motivate the use of irregular triangular elements to represent landscape evolution. The primary advantage is the variable resolution offered by the irregular domain, as opposed to the single resolution inherent in raster grids. Regions of high terrain variability can be modeled more precisely than areas of lower variability. Multiple resolutions could translate into computational savings as the number of nodes could be reduced in areas of low terrain variability (􏰴Goodrich et al., 1991)􏰪. While the TIN data structure is more complex to implement (􏰴Tucker et al., 2001b)􏰪, the reduction achieved in the number of model nodes results in a significant savings that can allow TIN-based landscape evolution models to operate over large regions􏰪. A second advantage is that TINs permit linear features to be preserved within the model mesh. For example, this allows the terrain to mimic natural terrain stream networks without introducing any raster artifacts inherent in grid methods. For a landscape model, TINs allow the stream network to be precisely depicted through the geomorphologic evolution.

Stream network algorithm

Independently of the mathematical complexity one wishes to use to simulate river evolution and associated sediment transport, all models efficiency will be mainly limited to their ability to compute river discharge and sediment yield at each point on the landscape. Thus, key to any implementation is to find the order in which rivers are progressively flowing from one node to the other.

The computation of river flow is typically an \(\mathcal{O}(n_p^2)\) problem where \(n_p\) is the number of points of the TIN on which it is calculated. Such a dependency greatly limits the applicability of landscape evolution model for large values of \(n_p\). Most of the flow routing algorithms, which have been routinely used, are \(\mathcal{O}(n_p\log{}n_p)\) (Wilson et al., 2008). Freeman (1991) developed a recursive algorithm that is \(\mathcal{O}(n_p)\). More recently, Braun & Willett (2013) proposed a new algorithm to compute drainage area that is also \(\mathcal{O}(n_p)\) and adapted to parallel computing architecture. According to the authors, their approach can be used to solve basic landscape equation (stream power law) on meshes that are several orders of magnitude finer than those commonly used with previous ordering algorithms. Here, the flow routing takes advantage of this new algorithm.

This algorithm falls in the SFD’s family (single flow direction approximation) where it is assumed that water goes down to the path of the steepest slope (O’Gallaghan and Mark, 1984; Gallant and Wilson, 2000). Assuming that the landscape is discretized in \(n_p\) nodes each node \(i\) has a single receiver \(r(i)\) except base level or local minima nodes which are their own receiver. Then, the approach consists in inverting the list of receiver nodes to find the list of donor nodes. Nodes with no donors are thus the ones to start the computation of the flow routing. From the receiver information and the number of neighbors to a given node \(i\), the procedure is based on inverting \(r(i)\) into donor array \(D(i,j)\) and number of donors per node array, \(d_i\). This information is then stored in a single dimension array \(D_i\) and an index table/array \(\delta_i\). Next a stack (\(s(i)\)) that contains the list of nodes per catchment is built recursively. Note only are the nodes within the stack or vector grouped by catchment but they are also ordered from downstream to upstream nodes. By inverting the order of the nodes on the stack, one go through the \(n_p\) nodes in an order that is appropriate to compute river discharge and associated sediment transport. The succession of operations used in the algorithm is fully explained in Braun and Willett (2013).

Once the stack is built, the water received from net precipitation \(\nu(i)\) for each node \(i\) is given by: \[\phi(i) = a(i) \nu(i) , \text{ for } i=1,\dots,n_p\] with \(a(i)\) the Voronoi area associated with node \(i\). The computation of the flow discharge is then expressed in the following manner: \[\phi(r(s(i))) = \phi(r(s(i)))+\phi(s(i)), \text{ for } i=n_p,\dots,1,-1\] The algorithm can take into account any spatial variation in rainfall/precipitation such as that resulting from an orographic model. The algorithm could also take into account more complex situations, involving water storage and/or release between adjacent cells.

Sub-catchment delineation

For numerical experiments with simple initial conditions (boundary, topography and tectonic history), the computation of the stack could be used to partition the nodes in a parallel architecture directly by distributing the base level nodes across the processors. However for more complex simulations, using this approach will often result in unbalanced computational load and therefore in inefficiently parallelized models.

To overcome this issue, one can increase the level of granularity in the representation of the stack by identifying not only the different catchments but also by defining a series of sub-catchments for each of them. Based on the constructed stack (\(s(i)\)) from Braun & Willett (2013) ordering process, one can rapidly determine all of the nodes that drain through a given point (outlet) \(i\). This in turn can be used to define a stream order by assigning a number to each section of the flow network. Several types of stream classification have been developed over the years (Horton (1945), Strahler (1952), Shreve (1966)). Here, the Strahler stream classification system is used. In this classification, waterways are given an order according to the number of additional tributaries associated with each waterway (Strahler, 1952).

Using similar notation as the one proposed in Braun & Willett (2013), one can define \(o(i)\) the Strahler order of node \(i\) and first initialized the list as: \[o(i)=1, \text{ for } i=1,\dots,n_p\] Then using the inverted stack order (e.g. proceeding downstream), the operation consists in verifying for every nodes that there are at least two equal tributaries (donors) with same order value. If not, a order value corresponding to the highest order \(o_{h}\) is attributed to the considered node; if this is the case, the node’s order is increased by 1 and continues downstream with new order. Which is built for \(i=n_p,\dots,1,-1\) according to: \[o(s(i))= \begin{cases} 1,& \text{if } d_{s(i)}=0\\ o_{h}, & \text{if } \max_{0 \leq j \leq d_{s(i)}} o(D(s(i),j)) = o_{h} , \text{ for 1 donor } \\ o_{h} + 1, & \text{ otherwise } \end{cases}\] The computation of the stream order is performed \(n_p \times\) a set of simple operations and is thus \(\mathcal{O}(n_p)\).

Along the computation of Strahler numbers, sub-catchments are delineated based on increasing stream order values. Each sub-catchment is therefore defined by its outlet node. Sub-catchment outlets list \(L_{sc}\) is incrementally filled during the stream order computation process. It is worth noting that the sub-catchment outlets list is ordered in the same way as the inverted stack from Braun & Willett (2013) and therefore the outlets list follows an upstream to downstream organization. The total number of outlets and consequently the number of sub-catchment \(Nsc\) is deduced from the \(L_{sc}\) dimension.

The next operation consists in finding the set of nodes for each sub-catchment. From the SFD approximation, a node belongs to an unique sub-catchment number. From the stack \(s\) (e.g., containing the list of nodes from downstream to upstream), one can find for each node \(i\) its sub-catchment number \(n_{sc}(i)\) by looping through \(Nsc\). Nodes, for a given sub-catchment \(k\), are within the range \(s(L_{sc}(k)),\dots,n_{ds}-1\) where \(n_{ds}\) is the first node within the stack \(s\) which has an elevation lower or equal to the considered outlet elevation. In situation where a sub-catchment number has already been attributed to a given node, its value is kept unchanged. The following procedure summarizes the approach for \(i=s(L_{sc}(k)),\dots,n_{ds}-1\): \[n_{sc}(i)=k \rightarrow \begin{cases} \text{if } n_{sc}(i)\gt k, \\ \text{if } h(i) \gt h(s(L_{sc}(k))), \\ \text{if } n_{ds} \leq s(n_p) \end{cases}\] where \(h\) is the elevation for a considered node.

The channel network represents an acyclic, directed graph, such that reaches and junctions are ordered from upstream to downstream. With the additional effort in defining not only base-level induced catchment but also sub-catchments, we will show in the next section how it is now possible to perform a domain partitioning which balances the number of TIN nodes across processors.