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/containers/xarray.hpp"
#include "xtensor/generators/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]\).

Flow Kernels

Flow accumulation is a simple way to compute values by traversing the flow graph from upstream to downstream nodes. This could be limited for more advanced use cases, though.

Fastscapelib also supports “flow kernels”, which consist of user-provided (Python) functions in which one can compute the value of one or more variables at one node of the graph - and/or at its receiver or donor nodes - depending on the values of variables evaluated at those nodes and also depending on the value of some arbitrary parameters. The so-called “flow kernel” functions are executed for each node of a flow graph traversed in a given direction (see FlowGraphTraversalDir), updating the output variable values between two executions of the kernel at adjacent nodes of the graph.

Flow kernels are efficient and extremely flexible. The kernel functions written in pure-Python are jit-compiled using numba before being called from within C++, therefore making their application very fast. In certain cases it is even possible to execute those functions in parallel (multi-thread), which may yield some nice speed-up for functions that are expensive to evaluate at a single node.

Warning

This feature is EXPERIMENTAL and for advanced usage!

Flow kernels are currently only available for the Python bindings and require numba.

Applying flow kernel functions in parallel (n_threads > 1) should be safe in most cases, although race conditions may still occur depending on how the kernel function is implemented.

Applying a kernel function in parallel may yield poorer performance than applying it sequentially, especially when the kernel function is cheap to evaluate at a single node. From our experience the benefit of parallelization may be difficult to evaluate.

numba has some complex internal optimizations to speed up the execution of the jit-compiled code. As a consequence, subtle differences in the kernel function’s code may significantly affect the performance.

As a simple example, basic flow accumulation (computation of drainage area) is implemented below as a flow kernel function. See also here for a more advanced example that re-implements the Stream-Power law as a flow kernel function (using the convenient FlowKernelEroder built on top of flow kernels).

Kernel function

Let’s first create the kernel function. It consists of an arbitrary function that accepts a single argument.

C++

// Flow kernels are not yet available in C++.

PYTHON

def accumulate_kernel_func(node):
    r_count = node.receivers.count
    if r_count == 1 and node.receivers.distance[0] == 0.0:
        # base level node
        return
    for i in range(r_count):
        weight = node.receivers.weight[i]
        node.receivers.drainage_area[i] += node.drainage_area * weight

The node argument passed to the kernel function above is a special object used for getting or setting the value(s) of any arbitrary variable at the current visited node with node.<variable_name>, at its receivers with node.receivers.<variable_name> or at its donors with node.donors.<variable_name>. In the example above there is only one user-defined variable drainage_area.

In addition to user-defined variables, the following attributes are defined:

• node.receivers.count: the number of receiver nodes.

• node.receivers.distance: the distance between the current node and each of its receivers (array of size=node.receivers.count).

• node.receivers.weight: the flow partitioning weight for each of the receivers (array of size=node.receivers.count).

• node.donors.count: the number of donor nodes.

Kernel and Kernel Data Objects

Next, create_flow_kernel() is used to compile the kernel function and create the flow kernel (data) objects.

C++

// Flow kernels are not yet available in C++.

PYTHON

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

kernel, kernel_data = create_flow_kernel(
    flow_graph,
    accumulate_kernel_func,
    spec=dict(drainage_area=nb.float64[::1]),
    outputs=["drainage_area"],
    n_threads=1,
    apply_dir=fs.FlowGraphTraversalDir.BREADTH_DOWNSTREAM,
)

A few arguments are passed to create_flow_kernel(): an existing flow graph object, the kernel function, the variable specifications (i.e., here above the “drainage_area” variable and its array-like numba type), the list of output variables (here “drainage_area” is computed by the kernel) as well as the direction in which to traverse the flow graph. Additionally, n_threads=1 means that the kernel will be applied sequentially along the graph nodes.

create_flow_kernel() returns two objects:

• kernel: a NumbaFlowKernel object that is a proxy to the jit-compiled kernel function.

• kernel_data: a NumbaFlowKernelData dict-like object that is a proxy to the kernel’s data (input/output variables and parameters).

Apply the Kernel

Finally let’s apply the kernel with apply_kernel(), which requires passing both the kernel and kernel data objects created at the previous step.

Prior to applying the kernel, it is required to create the input/output variables (below drainage_area is initialized with the grid’s node area values) and bind them to the kernel data via it’s bind() method.

C++

// Flow kernels are not yet available in C++.

PYTHON

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

drainage_area = grid.nodes_areas()
kernel_data.bind(drainage_area=drainage_area)
flow_graph.apply_kernel(kernel, kernel_data)

After applying the kernel, the drainage_area variable should have the correct values resulting from the accumulation / traversal.

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() and flow kernels 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.