FastFlow: GPU Acceleration of Flow and Depression Routing for Landscape Simulation

While the task of computing water flow over a terrain is a mainstay of hydrology and geomorphology, it is also directly applied in com-

puter graphics to the modeling of natural phenomena, in particular the generation of terrains and placement of vegetation.

Terrain Simulation. Methods for generating terrain for use in computer graphics can be broadly categorised [GGP⁺19] into procedural (using algorithmic rules to mimic emergent properties), example-based (learning structure from existing terrain data), and simulation (mathematically emulating natural processes). Within terrain simulation, fluvial erosion is recognized as a primary force in shaping the topography of mountains. Early erosion methods in computer graphics [MKM89] simulated water dynamics directly using shallow water equations [Ben07] or smoothed particle hydrodynamics [KBKv09]. Unfortunately, water dynamics are very short term and need to be applied many thousands of times to capture long-term erosion processes.

Instead, large-scale terrain erosion models [CBC⁺16] take their inspiration from geomorphology [BW13] and directly compute total water discharge by accumulating precipitation from high (mountain ridges) to low elevations (the sea). This flow accumulation is communicated across large stretches of the terrain and is thus less local and less amenable to parallelisation [VBHS11] than direct water dynamics. However, this is more than offset by the sheer number of iterations required for a dynamics solution to reach steady-state. Schott et al. [SPF⁺23] provide an approximate variant of discharge-based Stream Power erosion that propagates discharge by a few cells on each time step. While this strategy is trivially parallelizable, it requires a stable river network and thus precludes outside terrain forces such as time-dependent tectonics or sediment deposition. Furthermore, the underlying explicit time-stepping scheme requires many timesteps, while our method is amenable to an implicit scheme that overcomes this constraint.

Ecosystem Simulation. In ecosystem simulation, soil moisture is one of the viability criteria for plant growth. Here, flow maps based on monthly precipitation patterns [GLCC17] are used to supply moisture inputs for plant growth. This has a particularly strong effect in riparian areas. The most commonly used per-cell moisture proxy, the Topographic Wetness Index [RKK121], is based on a combination of slope and catchment area, with the latter derived directly from flow routing. Unfortunately, ecosystem simulations suffer from long run times on the order of hours and one of the roadblocks is flow calculation, particularly if the simulation is run at a weekly or daily time-step granularity.

Flow Routing. The problem of calculating water discharge over a terrain is well-studied and has given rise to many competing algorithms. Some methods are adapted to regular grid [GM97, OM84], and others to triangulated irregular network (TIN) [Ban07]. In the case of grids, the connections between cells can be either 4- or 8-way, depending on whether or not the diagonals are included. While our method is equally well-suited to TINs and 4- and 8-connectivity grids, for implementation purposes we currently focus on optimizing for 4-connectivity grids.

A further distinction lies in how these techniques direct flow to neighboring destination cells from a given source. The most common choice is Single Flow Direction (SFD), in which water always flows toward one of the lower neighbors of a given cell. The mechanism for selecting among lower cells matters in landscape simulation [TGS24]. A stochastic selection leads to more natural-looking results than, for instance, always selecting the neighbor with the lowest elevation. In contrast, Multiple Flow Directions (MFD) splits water across multiple destinations in function of the difference in elevation [QBCP91, Coa20]. Our method handles all SFD variants and is agnostic to the choice of the single destination. However, achieving $\backslash mathcal\{O\}(\backslash log(n))$ complexity is enabled by a tree data structure formed by SFD connections and this is not trivially extendable to the Directed Acyclic Graph (DAG) required by MDF.

Figure 2: (a) The discharge $q\_c$ corresponds to the integral of the precipitation $p$ upstream of a cell $c$ (Image courtesy of [CBC⁺ 16]). (b) Notation for the connectivity between terrain cells. (c) Two levels of graph abstraction. The stream trees cover the terrain grid and abstract the connections from each node to its recipient. At a higher level, the depression graph connects basins, which are regions of the terrain sharing a common local minimum. This depression graph is integral to the depression routing algorithm, which removes local minima: we update the recipients guided by the depression tree rooted in an outflow basin, defined as the minimum spanning tree of the depression graph.

There have been previous efforts to accelerate this problem by exploiting parallelism, either on CPU [Bar17] or GPU [Bar19, SPF²³]. These approaches parallelize the propagation of flow among independent flow paths but do not accelerate propagation within the paths themselves, leading to $\backslash mathcal\{O\}(l)$ iterations, where $l\; \backslash approx\; \backslash sqrt\{n\}$ is the length of the longest channel and $n$ the number of cells. In contrast, our approach requires $\backslash mathcal\{O\}(\backslash log(n))$ iterations.

Depression Routing. Crucially, previous GPU algorithms do not account for depressions or local minima in the topography, which are prevalent in terrains that have not been subject to depression removal [WQZ19]. In nature, rather than breaking the flow path, such depressions progressively fill and eventually overflow. Depression routing, therefore, aims to find paths out of all depressions, such that water can flow from any point on the terrain to a boundary.

Barnes et al. [BLM14] proposed the use of a priority queue in progressively flooding from the boundary to the interior, filling depressions in the process, leading to a $\backslash mathcal\{O\}(n\; \backslash log(n))$ algorithm. Subsequent work reduced dependence on the priority queue data structure [ZSF16, WZF18] and investigated CPU parallelization [DLM11, Bar16, ZLFS17, BCW20]. Ultimately, the priority queue was dispensed with altogether by exploiting the observation that depression routing is an instance of a Minimum Spanning Tree search, which can be found in linear time on a planar graph [CBB19]. We build on this idea and adapt it to the GPU.