Flow Routing

Flow routing consists of finding the paths along which some material (e.g., water, sediment) “flows” over the topographic surface, from upslope to downslope areas. It is an important component of most Landscape Evolution Models (LEMs) ; some of the main driver processes like river channel erosion depend on it.

Flow routing is also a complex problem that requires advanced algorithms and (graph-based) data structures in order to compute the flow paths more-or-less accurately, which often represents the main bottleneck in landscape evolution simulations.

Fastscapelib implements a collection of efficient, state-of-the-art algorithms for flow routing (e.g., [BW13], [BLM14], [CBB19]). Those algorithms are optimized for LEMs (i.e., require updating the flow paths repeatedly) and are accessible through a flexible and user-friendly API.

Flow Graph

Flow paths are contained in a “flow graph” data structure (more precisely, a Directed Acyclic Graph or DAG) that connects each grid node to its flow receiver(s) (neighbors).

The flow graph data structure is distinct although closely related to a grid (Figure 1):

• each graph node correspond to a grid node

• graph edges together form a subset of the grid (implicit or explicit) edges

Figure 1: Two kinds of flow graphs (left: single flow direction, right: multiple flow direction) on top of a raster grid with 8-node connectivity. Dots are grid (and graph) nodes, dots with a different color at the borders are the base level nodes, thin dashed lines represent the grid connectivity and thick arrows are the flow paths (graph edges).

A flow_graph (C++) or FlowGraph (Python) object can be created from any grid object (see Section Grids), e.g., in the example below from a raster grid:

C++

#include "fastscapelib/flow/flow_graph.hpp"
#include "fastscapelib/flow/flow_router.hpp"
#include "fastscapelib/grid/raster_grid.hpp"

namespace fs = fastscapelib;

fs::raster_boundary_status boundaries{ fs::node_status::fixed_value };
fs::raster_grid grid({ 100, 100 }, { 200.0, 200.0 }, boundaries);

fs::flow_graph<fs::raster_grid> graph(grid, { fs::single_flow_router() });

PYTHON

import fastscapelib as fs

grid = fs.RasterGrid([100, 100], [200.0, 200.0], fs.NodeStatus.FIXED_VALUE)

graph = fs.FlowGraph(grid, [fs.SingleFlowRouter()])

A sequence of flow operators must be passed as second argument to the flow graph constructor in order to define the flow routing strategy (see Sections Flow Operators and Flow Routing Strategies).

Base Level Nodes

The flow routing algorithms implemented in Fastscapelib require one or more nodes of the graph to be set as “base level” (or outlet) nodes (Figure 1).

These specific nodes may collect input flow but cannot release any output flow (it could also mean that the flow is leaving the modeled domain or sub-domain, e.g., surface land water entering the sea).

When creating a new flow graph, all grid nodes with status FIXED_VALUE are set as base level nodes by default.

It is possible to set alternative base level nodes via the set_base_levels() method (C++) or the base_levels property (Python). This can be done each time prior to updating the flow paths in the cases where the base levels are evolving over time (e.g., inner lake, sea or ocean dynamics, planet craters).

Mask

A binary mask can be set via the set_mask() method (C++) or the mask property (Python) to exclude grid nodes from the computation of the flow graph. By default, all grid nodes are included in the flow graph (no mask).

Setting a mask is useful if the modeled domain encompass multiple sub-domains such as land vs. ocean. Those sub-domains may each have their own flow graph with different flow routing strategies, base level nodes and masks.

Like for the base level nodes, a flow graph’s mask can be set or updated each time prior to updating the flow paths, e.g., according to sea level change over time.

Initializing / Updating the Flow Paths

Creating a new flow graph object doesn’t compute any flow path yet. To (re)compute the flow paths, the update_routes() method must be called with an input topographic surface (i.e., an elevation field defined on the grid):

C++

#include "xtensor/xarray.hpp"
#include "xtensor/xrandom.hpp"
#include "fastscapelib/flow/flow_graph.hpp"
#include "fastscapelib/flow/flow_router.hpp"
#include "fastscapelib/grid/raster_grid.hpp"

namespace fs = fastscapelib;

fs::raster_boundary_status boundaries{ fs::node_status::fixed_value };
fs::raster_grid grid({ 100, 100 }, { 200.0, 200.0 }, boundaries);

fs::flow_graph<fs::raster_grid> graph(grid, { fs::single_flow_router() });

xt::xarray<double> elevation = xt::random::rand<double>(grid.shape());

// update_routes returns a const reference!
const auto new_elevation = graph.update_routes(elevation);

PYTHON

import numpy as np
import fastscapelib as fs

grid = fs.RasterGrid([100, 100], [200.0, 200.0], fs.NodeStatus.FIXED_VALUE)

graph = fs.FlowGraph(grid, [fs.SingleFlowRouter()])

elevation = np.random.uniform(size=grid.shape)

new_elevation = graph.update_routes(elevation)

update_routes() returns another elevation field, which may differ from the input elevation field depending on the applied flow operators (e.g., some operators fill the closed depressions found in the input topography).

Note

update_routes() never updates in-place the values of the input elevation field.

Flow Operators

Flow can be routed over the topographic surface in many different ways ; choosing one approach over another highly depends on the case studied. Fastscapelib relies on the concept of “flow operators” that provide a flexible and convenient solution for implementing simple to advanced flow routing strategies.

A flow operator is a “routing unit” that:

• may read and/or modify in-place the flow graph instance to which it has been attached

• may read and/or update (a copy of) the elevation values passed to the update_routes() method

• may expose its own parameters

There are currently three categories of operators (see the C++ and Python API reference for a full list of available operators).

Flow Routers

Flow router operators generally compute new flow paths from scratch and fully update the flow graph. The updated flow graph has a defined FlowDirection: either SINGLE with one unique receiver per node or MULTI where the flow is partitioned among multiple receiver nodes with fixed or variable weights (Figure 1).

C++

#include "fastscapelib/flow/flow_router.hpp"

namespace fs = fastscapelib;

fs::single_flow_router::graph_updated  // true
fs::single_flow_router::out_flowdir    // fs::flow_direction::single

fs::multi_flow_router::graph_updated   // true
fs::multi_flow_router::out_flowdir     // fs::flow_direction::multi

PYTHON

import fastscapelib as fs

fs.SingleFlowRouter.graph_updated   # True
fs.SingleFlowRouter.out_flowdir     # fs.FlowDirection.SINGLE

fs.MultiFlowRouter.graph_updated    # True
fs.MultiFlowRouter.out_flowdir      # fs.FlowDirection.MULTI

Sink Resolvers

Sink resolver operators may either update the flow graph or the topographic surface (or both) so that no flow is trapped in closed depressions (i.e., every flow path is ensured to reach one of the base level nodes).

Some operators like PFloodSinkResolver only update the topographic surface and do not require any pre-existing flow paths, while other operators like MSTSinkResolver require single flow paths that will be updated in-place.

C++

#include "fastscapelib/flow/sink_resolver.hpp"

namespace fs = fastscapelib;

fs::pflood_sink_resolver::graph_updated      // false
fs::pflood_sink_resolver::elevation_updated  // true
fs::pflood_sink_resolver::in_flowdir         // fs::flow_direction::undefined
fs::pflood_sink_resolver::out_flowdir        // fs::flow_direction::undefined

fs::mst_sink_resolver::graph_updated         // true
fs::mst_sink_resolver::elevation_updated     // true
fs::mst_sink_resolver::in_flowdir            // fs::flow_direction::single
fs::mst_sink_resolver::out_flowdir.          // fs::flow_direction::single

PYTHON

import fastscapelib as fs

fs.PFloodSinkResolver.graph_updated       # False
fs.PFloodSinkResolver.elevation_updated   # True
fs.PFloodSinkResolver.in_flowdir          # fs.FlowDirection.UNDEFINED
fs.PFloodSinkResolver.out_flowdir         # fs.FlowDirection.UNDEFINED

fs.MSTSinkResolver.graph_updated       # True
fs.MSTSinkResolver.elevation_updated   # True
fs.MSTSinkResolver.in_flowdir          # fs.FlowDirection.SINGLE
fs.MSTSinkResolver.out_flowdir         # fs.FlowDirection.SINGLE

Flow Snapshots

FlowSnapshot is a special operator that can be used to read and save intermediate states of the flow graph and/or topographic elevation between two other operators.

C++

#include "fastscapelib/flow/flow_router.hpp"

namespace fs = fastscapelib;

fs::flow_snapshot::graph_updated      // false
fs::flow_snapshot::elevation_updated  // false

PYTHON

import fastscapelib as fs

fs.FlowSnapshot.graph_updated       # False
fs.FlowSnapshot.elevation_updated   # False

Flow Routing Strategies

A flow routing strategy is defined by an ordered sequence of flow operators passed to the FlowGraph constructor. Calling update_routes() will execute all the operators one after each other.

Here below are a few examples from simple to more advanced strategies (grid setup is skipped for more clarity, see the flow graph full example above).

Single Direction Flow Routing

The simplest example is to compute the flow paths along the steepest descent. It only requires the SingleFlowRouter operator.

C++

fs::flow_graph<fs::raster_grid> single_graph(grid, { fs::single_flow_router() });

PYTHON

single_graph = fs.FlowGraph(grid, [fs.SingleFlowRouter()])

Applied on a raster grid with 8-node connectivity, this is equivalent to the so-called “D8” algorithm [OCallaghanM84].

Flow Routing Across Closed Depressions

The priority flood sink resolver PFloodSinkResolver [BLM14] can be used to fill the closed depressions in the topography before computing the flow paths.

C++

#include "fastscapelib/flow/sink_resolver.hpp"

fs::flow_graph<fs::raster_grid> single_graph_nosink(
    grid, { fs::pflood_sink_resolver(), fs::single_flow_router() });

PYTHON

single_graph_nosink = fs.FlowGraph(
    grid, [fs.PFloodSinkResolver(), fs.SingleFlowRouter()]
)

Flow Routing Across Closed Depression 2

Alternatively, MSTSinkResolver [CBB19] can be used after computing the flow paths to re-route the flow trapped in closed depressions.

C++

#include "fastscapelib/flow/sink_resolver.hpp"

fs::flow_graph<fs::raster_grid> single_graph_nosink2(
    grid, { fs::single_flow_router(), fs::mst_sink_resolver() });

PYTHON

single_graph_nosink2 = fs.FlowGraph(
    grid, [fs.SingleFlowRouter(), fs.MSTSinkResolver()]
)

Compared to PFloodSinkResolver, MSTSinkResolver can be much faster to execute when the number of closed depressions found in the topography is small (often the case in LEM simulations already after a few time steps if no process is creating new closed depressions). However, MSTSinkResolver requires pre-computed single flow paths while PFloodSinkResolver has no specific requirement.

Multiple Direction Flow Routing (with Snapshots)

Here is a more advanced example where depression-free multiple direction flow paths are obtained after applying several operators. Two graph snapshots are also saved: the single direction flow paths before and after resolving closed depressions, respectively.

C++

#include <memory>
#include "fastscapelib/flow/flow_snapshot.hpp"

auto mrouter_ptr = std::make_shared<fs::multi_flow_router>(1.0);

fs::flow_graph<fs::raster_grid> multi_graph(
    grid,
    { fs::single_flow_router(),
      fs::flow_snapshot("single")
      fs::mst_sink_resolver(),
      fs::flow_snapshot("single_nosink"),
      mrouter_ptr });

PYTHON

mrouter = fs.MultiFlowRouter(1.0)

multi_graph = fs.FlowGraph(
    grid,
    [
        fs.SingleFlowRouter(),
        fs.FlowSnapshot("single"),
        fs.MSTSinkResolver(),
        fs.FlowSnapshot("single_nosink"),
        mrouter,
    ],
)

Note that flow operators may be instantiated outside of the flow graph constructor, like MultiFlowRouter in the example above. This is useful for changing parameter values between two successive calls to update_routes(), e.g.,

C++

multi_graph.update_routes(elevation);

mrouter->m_slope_exp = 1.5;

multi_graph.update_routes(elevation);

PYTHON

multi_graph.update_routes(elevation)

mrouter.slope_exp = 1.5

multi_graph.update_routes(elevation)

The flow partitioning method used in MultiFlowRouter is equivalent to Quinn et al. [QBCP91] method (slope_exp = 1) or Holmgren [Hol94] method.

Graph snapshots are separate instances of FlowGraph accessible via the graph_snapshot() method of the main graph instance. They are read-only.

C++

auto snapshot = multi_graph.graph_snapshot("single");

snapshot.update_routes(elevation)  // Error (read-only)!

PYTHON

snapshot = multi_graph.graph_snapshot("single")

snapshot.update_routes(elevation)  # Error (read-only)!

Flow Accumulation

After computing the flow paths, a FlowGraph object is ready to accumulate some locally produced quantity or flux along the network via the accumulate() method. This is handy for computing a range of internal variables, model outputs and/or diagnostics during a simulation. Flow accumulation requires a source term, which may be spatially uniform or variable.

Note

The input source term passed to accumulate() must expressed in units per-area. It will be integrated uniformly over the area surrounded by each grid node given by nodes_areas().

A few examples:

• drainage area (or upslope contributing area)

C++

auto drainage_area = graph.accumulate(1.0);

PYTHON

drainage_area = graph.accumulate(1.0)

Where the source term is equal to 1 (unit-less) so that drainage_area has dimensions \([L^2]\) and only results from the accumulation of the node (cell) areas along the flow paths.

• water discharge computed from local surface runoff rate

C++

xt::xarray<double> runoff_rate = xt::random::rand<double>(grid.shape());

auto discharge = graph.accumulate(runoff_rate);

PYTHON

runoff_rate = np.random.uniform(size=grid.shape)

discharge = graph.accumulate(runoff_rate)

Where runoff_rate has dimensions \([L/T]\) and discharge has dimensions \([L^3/T]\).

• sediment volume computed from vertical (local) erosion

C++

double dt = 1e3;
xt::xarray<double> erosion = channel_eroder.erode(elevation, drainage_area, dt);

auto sediment_vol = graph.accumulate(erosion);

PYTHON

dt = 1e3
erosion = channel_eroder.erode(elevation, drainage_area, dt)

sediment_vol = graph.accumulate(erosion)

Where erosion has dimensions \([L]\) and sediment_vol has dimensions \([L^3]\).

Advanced Usage

For advanced use cases it is possible to have read-only access to the graph internal data via the impl() method, which returns the FlowGraphImpl instance attached to the flow graph.

Warning

Relying on graph internal data for implementing custom logic in 3rd-party code should be done only when there is no alternative. Always prefer accumulate() when possible. Feedback and/or feature requests are also welcome.

FlowGraphImpl is an implementation detail and shouldn’t be taken as stable API. In C++ this class is under the fastscapelib::detail namespace.

Graph internal data should never be updated directly, or at your own risk! Only flow operators may do this internally. If direct update is possible this is probably a bug.