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

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.

Figure 8: A 300,000 year-old mountain, generated on a $2048\; \backslash times\; 2048$ grid in under 10 seconds, incorporating fluvial erosion, sediment deposition, and river and lake formation.

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].

Figure 9: River and lake identification on an alpine terrain dominated by a flat glacial valley (left), and a complex estuary (right).

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.

Figure 10: Sediment deposition acts as a shield against erosion leading to more pronounced slopes.

Figure 11: The visual results of implicit and semi-implicit deposition are on par (with 100 iterations at $\backslash Delta\; t\; =\; 2000$, $\backslash Delta\; x\; =\; 16$, and a resolution of $2048^2$).

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.

Figure 12: Approximate flow routing [SPF*23] (left) versus ours (right), subject to noisy uplift with varying amplitude (top to bottom: $\backslash pm\; 0$, 1, and 2m per timestep). Compared to ours, the elevation with approximate routing depends on the magnitude of that noise.

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.

Figure 13: Our implicit solution with $\backslash Delta\; t\; =\; 20,000$ (right) provides a result similar to the explicit scheme with $\backslash Delta\; t\; =\; 1,000$ (left). The latter is at its upper stability limit.

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.

Figure 14: Erosion carving on an escarpment, without (left) and with (right) depression routing.

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.

Figure 15: The visual differences between our Single Flow Direction (SFD) algorithm and Multiple Flow Directions (MFD), which we are unable to handle.