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

Algorithm 6: Extracting a depression-free water surface

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

Output: Elevation of the water surface $w$

1  $w \leftarrow \text{copy of } z$ 
2 for  $i \leftarrow 1$  to  $\log_2(|\mathcal{T}|)$  do
3   foreach  $\text{cell } c \in \mathcal{T}$  in parallel do
4      $w[c] \leftarrow \max(z[c], w[\text{recipient of } c] + 2^{i-1}\epsilon)$ 
5      $\text{recipient of } c \leftarrow \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 $\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.

Figure 6: A 3D visualization of an ecosystem simulation. It shows a landscape with a river network. The river channels and immediate banks are colored in shades of green, black, and purple, indicating different species or vegetation types. The surrounding terrain is a light tan color. The perspective is from an elevated angle, showing the flow of water through the landscape.

Figure 6: Ecosystem simulation with water provided by flow routing. Certain species (colored in shades of green, black, and purple) favor river banks due to higher levels of available water.

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

$\frac{\partial z}{\partial t} = u - kQ^m \frac{\partial z}{\partial x}, \quad (2)$

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 $\partial z/\partial x$ . The discharge provided by flow routing is scale-independent, so to compensate we scale the discharge $Q$ by accumulating $\Delta x^2 p$ , where $\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:

$z_i[c] = z_{i-\Delta t}[c] - K(z_i[c] - z_i[\text{recipient of } c]), \quad (3)$

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

$z_i[c] = \alpha z_i[\text{recipient of } c] + \beta, \quad (4)$

Figure 7: Four sequential 3D terrain renderings showing the stages of an interactive authoring session. The first stage shows a user-drawn uplift. The second stage shows sediment deposition in dark brown and lakes/rivers in blue. The third stage shows the overall simulation duration increased. The fourth stage shows the final result.

Figure 7: Different stages in an interactive authoring session. The user paints uplift to provide an initial mountain structure, then changes the deposition parameters (sediment in dark brown and lakes and rivers in blue), and, finally, increases the overall simulation duration.

where $\alpha = K/(1+K)$ and $\beta = z_t - \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.