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

Rain falling on a landscape flows downstream and progressively accumulates into streams and rivers. The volume of water flowing through a river cross-section per unit of time is called the discharge. Equivalently, the discharge is the upstream integral of the precipitation map. Flow routing involves the computation of this discharge for each cell of the terrain grid.

The direction of the water flow is discretized on the basis of a per-cell recipient, $r\_c$, which is chosen among the neighboring cells of the terrain and obeys the condition $z\_c\; >\; z\_\{r\_c\}$. This allows several mechanisms for the choice of neighbour. While the different choices have their place and are compatible with our algorithm, we prefer a random selection, with a probability based on the elevation difference between the source and potential recipient [TGSC24].

Let us consider the cells of a terrain as nodes in a graph, and interpret the connection between a cell and its recipient as a directed edge. The strict monotonically decreasing elevation of recipients guarantees that the graph does not contain any cycles. Since each cell has a single recipient (equivalently, a single edge leaves each node), this graph represents a set of trees (an algorithmic forest).

We use the term stream tree for each connected subgraph, since it represents a set of water channels. This allows flow routing to be recast as an accumulation (or scan). For a given stream tree, the scan proceeds from its leaves (ridges of the terrain) to a single root (an outflow cell). We propose employing rake and compress [SAF05] algorithms to perform this task in parallel (Section 4).

We handle the user-defined outflow mask by excluding recipient edges for its cells so that they become root nodes in the stream forest. However, other cells may also lack recipients if none of their neighbours meet the $z\_c\; >\; z\_r$ condition, i.e., local minima and flat areas. In reality, instead of interrupting flow, these areas would fill and overflow. This necessitates a solution for depression routing.

We are now in a position to define what a depression means in our context. As a preliminary, we use the term basin for the set of terrain cells that share the same stream tree, effectively all belonging to the same catchment area and channeling water to a single root cell. If this root cell is an outflow, then the basin is classified as an outflow basin, otherwise, it constitutes a depression (Figure 2 (c)).

Depression routing is the task of finding reasonable water paths out of all depression basins. Unfortunately, this problem is non-local and therefore not amenable to trivial parallelization. For instance, an obvious first thought would be to connect a depression to an adjacent basin at the point of the lowest neighboring elevation, but this may give rise to a chain of connections that form a cycle, and therefore is unable to route the water toward outflow cells.

Two families of solutions have been proposed to tackle this issue. One option is to progressively flood the terrain from the outflow cells towards the interior, keeping track of cells to be processed through a priority queue [BLM14]. While this enables a simple and efficient single-threaded solution, it is not suited to parallelization on the GPU, because it imposes a sequential order, with each processed cell potentially changing the subsequent state of the queue. The other option is to build a depression graph by representing depressions as nodes and their adjacency relationships as edges. There are two subtleties. First, all outflow basins are grouped into a single node marked as the start. Second, edges are weighted according to the lowest altitude along their shared boundary. The path out of all depressions is then given by the Minimum Spanning Tree (MST) of the depression graph, which connects all depressions to the collection of outflow basins.

In Section 5, we detail our GPU implementation of depression routing and propose two strategies for correcting the flow path: either creating a jump connection between a local minimum in the current depression and a cell in the next basin or preserving the continuity by inverting flow for a subset of cells in the depression.

Input : Terrain cells $\backslash mathcal\{T\}$ with initial discharge (precipitation) $q\; =\; p$ and their donors
Output: Updated accumulated discharge

1 for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{T\}|)$ do
2   foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
3     foreach donor $d$ of cell $c$ do
4       if $d$ is a leaf then
5         $q\_c\; \backslash leftarrow\; q\_c\; +\; q\_d$
6         remove edge $d\; \backslash rightarrow\; c$
7       end
8     else if $d$ has a single donor then
9       $q\_c\; \backslash leftarrow\; q\_c\; +\; q\_d$
10      donor of $c\; \backslash leftarrow$ donor of $d$
11    end
12  end
13  if all donors are removed then
14    tag $c$ as leaf
15  end
16 end
17 end

With the exception of some cells marked by the outflow mask, stream trees cover all of the terrain, which makes the design of an efficient parallel solution to Equation 1 non-trivial.

Fortunately, this problem can be cast as an instance of Parallel Tree Accumulation, which aims to compute the cumulative sum of quantities in tree nodes, from leaf to root. This can be achieved with straightforward pointer jumping, but this may result in a $\backslash mathcal\{O\}(n)$ span due to the required atomic operations. An alternative is to use an Euler tour [SAF05] along the tree edges combined with pointer jumping, which reduces the span to $\backslash mathcal\{O\}(\backslash log\; n)$, since the number of edges is less than the number of nodes. However, this method typically requires $2\; \backslash log\_2\; n$ iterations, as each edge is traversed twice.

A still more efficient approach is rake-compress, which combines pointer jumping with leaf pruning to achieve a span of $\backslash mathcal\{O\}(\backslash log\; n)$. This method typically needs only $\backslash log\_2\; n$ iterations, making it faster in practice than the Euler tour. Sevilgen et al. [SAF05] provide a starting point in the form of a general rake-compress solution to graph accumulation problems for distributed architectures. This relies on the observation that leaves constitute half the nodes in a balanced binary tree. So, if we iteratively transfer accumulated discharge from a leaf to its recipient and then prune the leaf and do this in parallel across all leaves, we only need to iterate $\backslash log\_2\; n$ times where $n$ is the number of nodes in the terrain.

In general, river trees might be neither binary nor balanced. In the worst case, where a river tree is composed entirely of single donors, the iterations degrade to $\backslash mathcal\{O\}(n)$. To compensate, we use pointer jumping by re-routing the edges of cells whose donor itself has a single donor (parent). Effectively, the parent (donor) is skipped and an edge is formed directly from the grandparent (donor of the donor). Thereby, worst-case performance is restored to $\backslash log\_2(n)$ parallel iterations.

As shown in Algorithm 1, with an accompanying example in Figure 3, on each iteration and for each terrain cell $c$, we parse each donor $d$ in parallel, with up to 4 donors in our implementation. If $d$ is a leaf, we add its discharge to $c$ and remove it from the donor matrix. Alternatively, if $d$ has only a single donor, we again add its discharge to $c$ but in this case reconfigure the donor matrix so that $c$ instead points at the donor of $d$. Distinguishing between these cases requires a single lookup in the donor count list, e.g., leaves have a donor count of zero. Straightforward extensions of Algorithm 1 can be used for other flow accumulation functions, including those involving maxima or weighted sums of upstream values.

Note that some care must be taken in the implementation of Algorithm 1, which exhibits potential read-after-write (RAW) hazards when updating the discharge and recipients. An immediate solution is to write the updated $q$ and $r$ to temporary arrays and copy them back at the end of each iteration. Unfortunately, this forces the unnecessary transfer on each iteration of a large number of cells that are unmodified. Instead, we propose a per-cell ping pong scheme. We first allocate sufficient space to hold the state of two full terrains, $T\_A$ and $T\_B$. We then introduce an additional source buffer with one 8-bit integer per cell, initialized to 0. The sign of the integer at cell $c$ in the source indicates whether data should be read from $T\_A$ and written to $T\_B$, or vice-versa. When cells are updated,

the sign is flipped. To avoid RAW hazards within the same iteration, the remaining bytes of each source entry store the iteration at which $c$ was last updated: a thread will check the source to negate the read sign if the stored iteration is the current one. At the end of each iteration, the latest state of each cell is either in $T\_A$ or $T\_B$, as indicated by source. Hence, a final merge step is required before returning the result as a single tensor. Incorporating this optimization leads to a $14\backslash times$ speed up on a $1024\; \backslash times\; 1024$ terrain.

Input : Terrain cells $\backslash mathcal\{T\}$, non-outflow local minima $\backslash mathcal\{L\}$
Output: Per-cell basin identifiers bid

1 for $i\; \backslash leftarrow\; 0$ to $|\backslash mathcal\{L\}|\; -\; 1$ in parallel do
2   $bid$ of $\backslash mathcal\{L\}[i]\; \backslash leftarrow\; i$
3 end
4 for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{T\}|)$ do
5   foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
6     $bid$ of $c\; \backslash leftarrow\; bid$ of recipient of $c$
7     recipient of $c\; \backslash leftarrow$ recipient of recipient of $c$
8   end
9 end

graph to progressively detect the path of water and update recipients. The algorithm proceeds through the stages of basin identification, saddle edge selection, and donor-recipient re-routing.

Basin identification. Our algorithm begins by segmenting out each stream tree and associating with its cells a unique basin identifier (hereafter referred to as a basin-id, or bid in the pseudo-code of Algorithm 2). To this end, we assign to each non-outflow local minima a unique identifier – leaving a common outflow identifier for all outflow basins – and use pointer jumping to copy this identifier in parallel to all upstream nodes of its stream tree (Algorithm 2).

Saddle edge selection. We use the per-cell basin annotation to isolate cells lying along the border with neighboring basins. To do so, we tag as a border cell, any cell in the current basin with at least one adjacent cell belonging to another basin. Ultimately, one of these border cells is chosen as a saddle edge in the depression graph. We do not need to construct these edges explicitly as Boruvka's algorithm only requires the edge with the lowest weight.

For each border cell, this weight is set as the lowest elevation that will permit an overflow from it into an adjacent basin. Specifically, this is the maximum of the border cell's altitude and the lowest altitude among its neighboring cells. Then we choose the border cell with minimum weight per basin.

Although we could compute the minimum altitude using a variant of rake-compress in $\backslash mathcal\{O\}(\backslash log\; n)$, we found it more efficient to instead compute it in a single step with atomic operations. We also store the saddle, which is the source cell in the current basin with minimum weight, and the outlet, which is its lowest neighbouring cell across the basin border (see Figure 4).

Close-coupled loops can occur when two basins have saddle edges pointing at each other. We resolve this by simply deleting the edge with the higher saddle basin-id. We also apply a similar strategy to choose among saddles of identical elevation, which could lead to larger cycles. This is prevented by selecting saddles in lexicographic order (elevation concatenated with the basin-id of the outlet). This algorithm is summarized in Algorithm 3.

Re-routing recipients. Ordinarily, the next step in Boruvka's algorithm would be to collapse the newly incorporated edges. Due to

Input : Terrain cells $\backslash mathcal\{T\}$ with elevation $z$, set of depression basins $\backslash mathcal\{D\}$
Output: Per-depression saddle and outlet

1 // Compute border cells and store the border altitude in $z\text{'}$
2 foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
3   if $bid\_c\; \backslash neq\; bid\_n$ of any 4-neighbour $n$ of $c$ ($n\; \backslash in\; \backslash mathcal\{N\}\_c$) then
4     tag $c$ as a border cell
5     $znb\; \backslash leftarrow\; \backslash min(z\_n\; \backslash text\{\; for\; \}\; n\; \backslash in\; \backslash mathcal\{N\}\_c\; \backslash text\{\; such\; that\; \}\; bid\_c\; \backslash neq\; bid\_n)$
6     $z\text{'}\_c\; \backslash leftarrow\; \backslash max(z\_c,\; znb)$
7   end
8 end

9 // Saddles and outlets (with atomic lexicographic argmin)
10 foreach basin $d\; \backslash in\; \backslash mathcal\{D\}$ in parallel do
11   saddle of $d\; \backslash leftarrow\; \backslash text\{argmin\}((z\text{'}\_c,\; bid\_n)$ for each border cell $c$
12   and $n\; \backslash in\; \backslash mathcal\{N\}\_c$ such that $bid\_c\; =\; d$ and $bid\_n\; \backslash neq\; d$)
13   outlet of $d\; \backslash leftarrow\; \backslash text\{argmin\}(z\_n$ for each neighbor $n$ of the
14   saddle such that $bid\_n\; \backslash neq\; d$)
15 end

16 // Remove cycles
17 foreach basin $d\; \backslash in\; \backslash mathcal\{D\}$ in parallel do
18   $d\text{'}\; \backslash leftarrow$ basin of outlet of $d$
19   if $bid$ of outlet of $d\text{'}\; =\; bid$ of saddle of $d$ then
20     if $bid$ of outlet of  of saddle of $d$ then
21       Delete the saddle and outlet of $d$
22     end
23   end
24 end

our two-layer structure, with an implicit depression graph and explicit flow routing, this needs to be interpreted differently. Instead, we merge stream trees by adding and re-routing recipient links, using one of two strategies (Figure 5 (b)).

The first strategy, depression jumping, involves adding a recipient link directly from the depression's local minima to the newly established outlet cell. This is straightforward and efficient but leads to a jump discontinuity in the flow path, which, depending on the application, may be undesirable.

The second strategy, depression carving [Rie98], mimics the course of water channeling down from the saddle. We first compute the flow path that connects a saddle and outlet to the associated local minimum and then reverse the direction of all recipient links. In this way, we allow water to flow upstream from the local minimum to the outlet. The pseudo-code for these two strategies is presented in Algorithm 4.

A single iteration of basin identification, saddle edge selection, and re-routing of recipients ends up being insufficient. It will not necessarily connect all depressions to the set of outflow basins. However, on any given iteration this process removes at least half of the remaining depressions, meaning that a total of $\backslash log(L)$ iterations are required, where $L$ is the number of local minima (Algorithm 5).

Input : Terrain cells $\backslash mathcal\{T\}$, set of depression basins $\backslash mathcal\{D\}$
Output: Updated recipients

1 if Depression jumping variant then
2   foreach depression $d\; \backslash in\; \backslash mathcal\{D\}$ in parallel do
3     $l\; \backslash leftarrow$ local minimum of $d$
4     recipient of $l\; \backslash leftarrow$ outlet of $d$
5   end
6 else
7   // Depression Carving Variant
8   tag $\backslash leftarrow$ False for all cells
9   foreach depression $d\; \backslash in\; \backslash mathcal\{D\}$ in parallel do
10    tag[saddle of $d$] $\backslash leftarrow$ True
11    end
12    // recipient is local to next two loops
13    for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{T\}|)$ do
14      foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
15        if tag[$c$] then
16          tag[recipient of $c$] $\backslash leftarrow$ True
17          recipient of $c\; \backslash leftarrow$ recipient of recipient of $c$
18        end
19      end
20      foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
21        if tag[recipient of $c$] and $c$ is not a local minimum
22        then
23          recipient of recipient of $c\; \backslash leftarrow\; c$
24        end
25      end
26    end
27    foreach depression $d\; \backslash in\; \backslash mathcal\{D\}$ in parallel do
28      recipient of saddle of $d\; \backslash leftarrow$ outlet of $d$
29    end

Input : Terrain nodes $\backslash mathcal\{T\}$ with depressions; non-outflow local minima $\backslash mathcal\{L\}$
Output: Updated terrain graph with rerouted depressions

1 for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{L\}|)$ do
2   Construct basins $\backslash mathcal\{D\}$ and annotate the cells (Algorithm 2)
3   Find min edges, saddles and outlets (Algorithm 3)
4   Reroute recipients (Algorithm 4)
5 end

Each iteration, in turn, requires $\backslash log\_2(n)$ sub-iterations, where $n$ is the number of cells, yielding an $\backslash mathcal\{O\}(\backslash log(L)\; \backslash log(n))$ algorithm.

Implementation. As with flow routing, we prevent read-after-write hazards while minimizing the extent of copying. For basin identification, we store the basin-id only for the local minima and not for all cells, and we use the updated recipients at the end of Algorithm 2 as a per-cell pointer to the downstream local minima, which we interrogate to obtain the basin-id. This change relaxes the need for any copy during basin identification, as read-after-write no longer prevents the correct convergence of the algorithm. Furthermore, we observe that this operation is required only for the first iteration. In subsequent iterations, we need only update the pointers to reflect the new connections out of the local minima, yielding an overall complexity of $\backslash mathcal\{O\}(\backslash log(n)\; +\; \backslash log^2(L))$.

Another RAW hazard can occur in the depression carving variant for re-routing flow paths. To prevent this, we use a localized ping-pong strategy (similar to the one presented in Section 4), which is simplified by using the tag variable as a condition for copying.

Input : Terrain cells $\backslash mathcal\{T\}$ with elevation $z$, small slope $\backslash epsilon$

Output: Elevation of the water surface $w$

1 $w\; \backslash leftarrow\; \backslash text\{copy\; of\; \}\; z$
2 for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{T\}|)$ do
3   foreach $\backslash text\{cell\; \}\; c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
4     $w[c]\; \backslash leftarrow\; \backslash max(z[c],\; w[\backslash text\{recipient\; of\; \}\; c]\; +\; 2^\{i-1\}\backslash epsilon)$
5     $\backslash text\{recipient\; of\; \}\; c\; \backslash leftarrow\; \backslash text\{recipient\; of\; recipient\; of\; \}\; c$
6   end
7 end

River and lake modeling. To identify the surface of rivers and lakes we begin by filling depressions in the terrain. This is sometimes referred to as pit removal and has value as an algorithm in its own right. Our input is the terrain elevation grid, and a set of recipients identified during depression routing to provide the flow path with minimum energy cost. We call the water surface the elevations obtained after filling the identified depressions. This surface should obey two conditions: it should be as close as possible to the terrain, and at the same time be monotonically non-decreasing along all stream trees, which we intuitively expect since water flows downhill. To achieve this, in Algorithm 6 we introduce a variant of the basin identification algorithm. We follow the water path of each stream tree upwards from destination to sources, retaining as the water surface the maximum elevation previously reached. For completeness, we optionally allow a small slope $\backslash epsilon$ on the water surface as this is sometimes required in hydrology applications. Note that this variant is only possible with the depression carving strategy, as it requires a continuous path. Next, cells with water discharge above a specified threshold can be tagged as belonging to bodies of water, as illustrated in Figure 9 where the water surface and its discharge are highlighted.

Ecosystem simulation. The water routes derived through river and lake modeling have direct application to ecosystem simulation. Water, along with sunlight and warmth, is a key requirement for plant growth, with different species having varied adaptation to water availability. Modeling water discharge is therefore a crucial element in ecosystem simulations. To demonstrate applicability, Figure 6 shows the outcome of an ecosystem simulation for a biome in the Pyrenees mountains [PGG*24] at the 50 year mark, with water accumulation provided by our flow and depression routing. For clarity, pioneer species with strong drought tolerance have been removed and what remains are species, such as Sessile Oak and European Beech, that require the greater moisture found in riparian areas. Most ecosystem simulations use a month as the timestep granularity, due in part to the overhead of calculating water flow. The acceleration provided here creates an opportunity for shorter, weekly, or even daily, timesteps with non-uniform rainfall across the landscape [PMG*22] and hence greater simulation accuracy.

Terrain erosion. Flow routing is a critical component of terrain erosion, as it provides the mechanism by which material is worn away and transported, ultimately over geological time shaping the mountains themselves. Furthermore, even the ordering of recipient cells plays a role in the implicit integration of the stream power law erosion equation. This is expressed as [CBC*16]:

where the first term $u$ is the tectonically induced uplift (or growth rate of the mountain), and the second term is the erosion, dependent on the discharge $Q$, erosion coefficients $k$ and $m$, and the slope $\backslash partial\; z/\backslash partial\; x$. The discharge provided by flow routing is scale-independent, so to compensate we scale the discharge $Q$ by accumulating $\backslash Delta\; x^2\; p$, where $\backslash Delta\; x$ is the cell size and $p$ is the precipitation. The uplift is usually considered time-independent and therefore applied as a pre-process. This leaves an implicit solution to the second erosive part of Equation 2 as:

where $\backslash Delta\; t$ is the time step and $K\; =\; kQ^m\; \backslash Delta\; t/\backslash Delta\; x$. Eq. 3 can be further simplified to:

where $\backslash alpha\; =\; K/(1+K)$ and $\backslash beta\; =\; z\_t\; -\; \backslash Delta\; t[c](1+K)$.

As with water surface extraction, we can compute $z\_t$ using a simple modification of basin identification (see Algorithm 7). The outcome of applying our accelerated method of fluvial erosion is illustrated in Figure 8. Note that Equation 2 is a special case of the Stream Power Law with the slope exponent set to one. Further research is required to adapt our solution to the general case, for instance with a global Newton-Raphson's algorithm.

Input : Terrain cells $\backslash mathcal\{T\}$ with elevation $z$, per cell coefficients $\backslash alpha$ and $\backslash beta$

Output: Eroded terrain $z\text{'}$

1 $z\text{'}\; \backslash leftarrow$ copy of $z$
2 for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{T\}|)$ do
3   foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
4     $z\text{'}[c]\; \backslash leftarrow\; z\text{'}[c]\; +\; \backslash alpha[c]z\text{'}[\backslash text\{recipient\; of\; \}\; c]$
5     $\backslash alpha[c]\; \backslash leftarrow\; \backslash alpha[c]\backslash alpha[\backslash text\{recipient\; of\; \}\; c]$
6     recipient of $c\; \backslash leftarrow$ recipient of recipient of $c$
7   end
8 end

Fluvial erosion is not the only erosive process responsible for shaping terrain. Hillslope processes model the gradual accumulation of solid material at the base of mountains and hills [BS97] and is usually expressed as a diffusion equation. We follow Tzathas et al. [TGSC24] and approximate it by including additional terms in the Stream Power Equation, changing $kQ^m$ to $kQ^m\; +\; k\_t\; +\; k\_h\; A^\{-h\}$, where $k\_t$ and $h\; =\; 0.6$ are hillslope erosion parameters, $A$ is the drainage area (obtained via flow routing with precipitation set uniformly to $p\; =\; \backslash Delta\; x^2$). The parameter $k\_t$ accounts for slope-dependent effects (landslides, debris-flow), regrouped in computer graphics under the catch-all term thermal erosion [MKM89].

Sediment deposition. Finally, flow routing can more generally propagate any type of transportable material. A case in point is sediment deposition, which acts as a complement to hydraulic erosion. Incorporating sediment deposition requires an additional term in the Stream Power Law [YBG⁺19]:

where $k\_d$ is the sediment deposition coefficient and $Q\_s$ is sediment flux, obtained by accumulating the negative elevation balance $-\backslash Delta\; x^2\; \backslash frac\{\backslash partial\; z\}\{\backslash partial\; t\}$ downstream with flow routing. While Yuan et

al. [YBG⁺19] suggest a fully implicit solution that requires iterating between Equation 5 and the accumulation of $Q\_s$, we propose a semi-implicit variant in which $Q\_s$ is computed from the elevation balance of the previous time-step. While our solution introduces a time dependency into mass conservation, it does not require iterating, which leads to a faster algorithm with no visual difference.

We prototyped our algorithms in Python with PyTorch and optimized our implementation with custom CUDA kernels. In particular, we seek to strike a balance by executing high-level functions in PyTorch or TensorFlow and handling compute-intensive kernels with CUDA. The code for this paper is available at: https://gitlab.inria.fr/landscapes/fastflow. For ease of use, this includes our standalone unoptimized PyTorch implementation and a TensorFlow port, as well as our CUDA optimizations. The repository also includes examples of how to use our code to edit terrain with Houdini [Sid23], as demonstrated in the accompanying video.

In terms of optimization, we reduce memory allocations and copying and prevent CPU-GPU communication in the control flow. We found that allowing the GPU to run for the maximum number of iterations (e.g., $\backslash log\_2\; n$ for flow routing) was more efficient than checking for early convergence. The only exception is in depression routing, where we halt the algorithm when no depressions remain.

In this section, we showcase the application of our methods to interactive landscape simulation. In doing so, we compare against a recent approximate GPU solution for flow routing [SPF⁺23] and evaluate the impact of our new approach to sediment deposition.

As a starting point, Figure 1 shows a real digitized terrain, subsequently altered by fluvial erosion and sediment deposition, and overlaid with flow-dependent vegetation, lakes, and rivers, in under 1.5 seconds.

Interactive landscape authoring. To demonstrate the general applicability of our algorithms we provide extracts from an interactive landscape authoring session in Figure 7 and the accompanying video. Here, user-guided simulation is applied to a $512\; \backslash times\; 512$ terrain, which, thanks to GPU acceleration combined with implicit time-stepping, requires only 10 iterations and 0.1s to capture 700,000 years of geomorphological evolution. For comparison, a CPU implementation [CBC*16] requires 2.6s. The user first defines the primary mountains by progressively painting on an uplift map, and can then freely change simulation parameters (here the deposition constant). Finally, the user advances the age of the mountain by increasing the number of iterations to 100, producing the result on the far right in 0.7 seconds. Note that lakes disappear over time as a consequence of filling by deposition and uplift, combined with erosion, which gradually removes the obstructions between them.

Our method scales to higher-resolution landscapes, as illustrated in Figure 8 for a $2048\; \backslash times\; 2048$ terrain with 300,000 simulated years (200 iterations) generated in under 10 seconds.

Rivers and lakes In Figure 9, we demonstrate how our algorithms can be used on complex large-scale terrains to demarcate lakes and rivers that are registered correctly with existing erosion lines and depressions. On the left in Figure 9, a mountainous region in the Alps transitions into an urbanized flat glacial valley. On the right, a tidal estuary around Mont Saint-Michel, France, leads into the sea. They contain 383,918 and 8,445,644 basins, respectively. In both cases, the elevation data was obtained from the IGN RGE ALTI Digital Elevation 1m dataset [IGN22].

Sediment deposition. Figure 10 highlights the impact of sediment deposition. This increases the diversity of landscape patterns. In particular, sediment deposition, where present, dampens erosion, and leads to steeper slopes. Finally, we show in Figure 11 the visual impact of our semi-implicit deposition strategy. Our results are similar to a fully implicit variant [YBG*19]. We note, however, that our gain in performance comes at the cost of some material loss.

Comparative evaluation. An alternative strategy for landscape simulation on the GPU [SPF*23, JBC24] is to use explicit time-stepping for erosion combined with approximate flow routing, where the discharge is propagated by a single cell on each simulation step. This method requires many more iterations due to the stability conditions of the explicit scheme but does eventually converge to a solution qualitatively similar to ours. Crucially, there are several situations where this strategy cannot be applied. One case is the erosion of an existing terrain containing depressions. Another failure case (see Figure 12) occurs when noise is added to uplift to increase the diversity. The approximate solution relies on the assumption that the flow path weakly varies over time, which allows the approximation of the discharge to progressively improve. Adding noise invalidates this assumption and leads to a significant underestimation of drainage and, consequently, an underestimation of erosion. This is visible in Figure 12, where the maximal elevation depends on the amount of noise. In contrast, our combination of implicit-time-stepping and exact-flow does not suffer from this shortcoming.

Another consideration is the acceleration afforded by the larger timesteps of our implicit solution. We show in Figure 13 that a $512\; \backslash times\; 512$ terrain with $\backslash Delta\; x\; =\; 32\backslash text\{m\}$ obtained with an explicit scheme ($\backslash Delta\; t\; =\; 1,000$ years) is visually similar to the result of our implicit solution ($\backslash Delta\; t\; =\; 20,000$ years). Our implicit scheme allows an increase of the timestep by a factor of 20, which reduces, by the same factor, the iterations needed to achieve the same total geological timespan. With such timesteps, generating a 10 million year-old landscape requires 7.2s and .5s, for the explicit and our implicit scheme, respectively. Note that larger implicit timesteps are unconditionally stable, and only diverge slightly from the explicit solution. Note that explicit and implicit schemes in general do not yield identical solutions as they accumulate discretization error differently.

We further motivate the need for depression routing by recreating the escarpment scenario from Cordonnier et al. [CBB19]. This begins with a flat terrain in which the boundary edges are lowered by 100m. Figure 14 shows the results without uplift after 50 iterations of 5,000 years, both with and without depression routing. Depression routing prevents water loss, and therefore leads to a much deeper carving of the valleys.

Limitations. For efficiency, our flow routing relies on parallel operations on a tree data structure, which in turn requires a single recipient per cell (Single Flow Direction, or SFD). As a consequence, we cannot handle multiple recipients (Multiple Flow Directions, or MFD), unlike some other schemes [Bar19,SPF*23]. Figure 15 illustrates the differences in appearance between rivers formed by SFD and MFD, using a CPU variant of Barnes approach [Bar19]. From a visual standpoint, MFD is more blurred and SFD more sharply defined. From an application standpoint, SFD is favoured for terrain erosion and river modeling, while MFD is more used in ecosystem simulation. Furthermore, Algorithm 7 for implicit fluvial erosion assumes a linear correlation with slope, meaning it cannot be applied study non-linear behavior in geomorphology.

In this section, we compare our approach against other parallel GPU and CPU implementations of flow and depression routing.

Our performance tests were executed on a computer equipped with an Nvidia RTX A6000 GPU with 48GB of memory and used 20 cores of an Intel Xeon Gold CPU clocked at 2.10GHz with 128GB RAM.

For our test dataset, we incorporated a mix of real and synthetic terrains. We began by extracting a set of 10 real terrains from the IGN RGE ALTI Digital Elevation 1m dataset [IGN22] chosen on the basis of topographic variety. These 10 terrains were sampled at cell resolutions of $512\; \backslash times\; 512$, $1024\; \backslash times\; 1024$, $2048\; \backslash times\; 2048$ and $4196\; \backslash times\; 4196$ to support scaling experiments.

We also included 5 synthetic terrains at the same resolution levels, and each with a controlled proportion of depression coverage (at synth-1%, synth-5% and synth-15% levels). These were generated by erosion simulation (Section 6), followed by layering different amplitudes of uniformly distributed noise. Note that it is common practice to add this type of noise during erosion simulation to mimic natural stochastic processes.

Finally, we included a uniformly flat terrain as an extreme test case, since in this instance all cells are local minima and hence coded as depressions.