Physically-based analytical erosion for fast terrain generation

Mountain formation is generally associated with tectonic plate convergence, mantle dynamics, and/or volcanic activity. In the context of tectonic convergence, the eventual collision of the plates leads to a thickening of the crust accommodated by visco-elasto-plastic deformation, brittle rupture (faults), and folding of the rock layers [Avo03, GCG⁺09]. The thickening leads to a vertical upward motion of rocks towards the surface, referred to as rock uplift,

which competes with erosion and can lead to a progressive increase in the surface altitude, called surface uplift [EM90]. Rates of rock uplift and erosion vary in space and time and achieve values up to a few millimeters per year in some mountain ranges.

The uplift is counteracted by erosion, which impacts the slopes of the mountain and therefore its maximal elevation. Erosion comes from many factors: water, glaciers, landslides, wind, and even anthropic or biological impact. Many models in geomorphology consider only erosion by water, also called fluvial erosion. Indeed, the fluvial network is considered the backbone of landscapes, and fluvial incision dictates the rate of landscape erosion [Whi04]. While simple to model, fluvial erosion explains the main topographical characteristics of most mountain ranges and has been the dominant erosion factor over many geological periods - with the notable exception of the last million years, where the Quaternary saw an important increase in glacial erosion that leaves specific marks in high altitude [PMD01, ENPL09, SHV⁺12].

A first common modeling approximation is to consider surface evolution as a detachment-limited process (as opposed to transport-limited) where the evolution of surface elevation $z$ is directly related to the competition between rock uplift and erosion rate, and not to the capacity of the rivers to transport or deposit sediments [How94]. In this setting, erosion is generally described using the stream power incision model [HK83, How94, WT99, Lag14], where erosion rate is a power law of the surface slope $\backslash |\backslash nabla\; z\backslash |$, and drainage area $A$ that acts as a proxy for the river discharge. The drainage area $A(\backslash mathbf\{x\})$ is defined at a position $\backslash mathbf\{x\}$ as the area of the drainage basin - or catchment - upstream of $\backslash mathbf\{x\}$. Coupled with the uplift $u$, the Stream Power Law expresses the rate of change of surface elevation:

where $k$, $m$ and $n$ are erosion coefficients. Throughout the paper, we will use some of the common values: $m\; =\; 0.4$ and $n\; =\; 1$. The choice of $n\; =\; 1$, also commonly used in geomorphology, makes the equation linear and therefore simplifies the derivation of the analytical solutions. While this choice barely impacts the result as the valley profiles are mostly directed by the ration $m/n$, we acknowledge that the actual values of $m$ and $n$ remain an open question in geomorphology [Lag14].

Eqn. 1 is a linear hyperbolic Partial Differential Equation (PDE), which can be interpreted as an advection of the initial terrain $z\_0\; =\; z(t\; =\; 0)$ along the rivers, at a velocity $kA^m$. This equation admits analytical solutions in 1D that can be computed thanks to the method of characteristics [Luk74, RTP13], and used to model 2D heightmaps [Ste21].

Throughout our implementation, we use a regular grid as our data structure – a standard for terrains that eases the multigrid-inspired technique that we introduce below. This grid initially stores the user-specified uplift $u$, initial elevation $z\_0$, and boundary mask $b$: each cell where $b$ is True will force its altitude to remain constant: $z\; =\; z\_0$. At least one cell should be specified as a boundary to ensure that the erosion problem is well-posed. The grid is progressively augmented with all the intermediate values needed by our implementation, including the drainage area $A$. Our output is a regular grid containing the final elevation $z$.

Analytical solutions were exposed in the geoscience literature [Luk74, RTP13, Ste21] but with some restrictions: either these solutions are in 1D [Luk74, RTP13] and the proposed algorithms are not straightforward to extend to 2D, or the solutions are in 2D [Ste21] but are restricted to a simplified case where the initial elevation $z\_0$ as well as the bounds $z(b)$ are assumed to be zero, and the uplift is constant in space. In Section 4, we detail our first contribution which is the derivation and implementation of the analytical solutions in the general case, compatible with a 2D setting.

We follow [Ste21] and order the computation along the river network. This network consists of a set of trees that covers the terrain and represents the progressive merging of high-altitude small streams down to the larger rivers. We obtain the 2D analytical solution by solving the 1D problem on each upstream path embedded in the tree structures. Similarly to the previous implicit time stepping schemes for the Stream Power Law [BW13, CBC*16], we separate the computation into two parts: first we accumulate the drainage area by following the river directions from the high elevations (ridges) of the terrain to the bounds, then we evaluate the analytical solution upstream from the bounds to the ridges.

This strategy leads to a second issue that hinders the optimal complexity of our algorithm on a stream tree: the elevations obtained from the implicit solution depend on the ordering of the nodes in the river network, which itself depends on the elevations. In practice, this leads to the formation of discontinuities at the boundaries of drainage basins. A solution is to use a fixed point algorithm: start with the initial elevation $z\_0$, compute the corresponding drainage $A\_0$, and then iterate over a progressively improving sequence of elevations $z\_i$ and $A\_i$, that eventually converge to the desired solution $z$. A straightforward implementation of this algorithm converges [Ste21] but requires at worst a number of iterations proportional to the length of the longest river in the terrain. In Section 5, we propose two solutions to reduce the number of iterations (Figure 2):

therefore introduce a method inspired by multigrid to move the bounds across different scales.

We show how we can combine these two approaches to allow the user to specify the desired accuracy of the final terrain in a continuum between the geologically accurate solutions provided by multigrid and a result that preserves the hydrological features of the input terrain, enforced by our optimization.

Our last challenge lies in the fact that the solutions to the Stream Power Law lead to singularities close to the ridges: the amount of water decreases when we approach the higher parts of the mountain, which are therefore less impacted by erosion and eventually degenerate into infinitely steep slopes. In practice, fluvial erosion is negligible or even not applicable in these locations, in which other erosion processes dominate. In Section 6, we propose some geomorphological-inspired modeling strategies to mitigate this issue, by adding terms that mimic other erosion factors, such as landslides (or thermal erosion) [MKM89], or hillslope erosion [BS97].

Eqn. 1, with $n\; =\; 1$, can be rewritten as:

which falls in the category of hyperbolic partial differential equations and can be solved by the method of characteristics [RTP13,

Ste21]. This method consists of finding a family of space-time curves, called characteristic curves, along which the equation becomes an Ordinary Differential Equation [ZT86].

Our first observation is that the spatial component of the characteristic curve follows the direction of the topographic gradient $\backslash frac\{\backslash nabla\; z\}\{\backslash |\backslash nabla\; z\backslash |\}$. Intuitively, this means that the elevation at any point $x$ in the terrain will only depend on elevations downstream of $x$, by following the path of steepest gradients. We call this path a river path as this is the trajectory naturally followed by the water, and parameterize it by $s$, the distance between the bound ($s\; =\; 0$) and any point on the path (Figure 3, left).

Along this path, Eqn. 2 simplifies as:

where $a(s)\; =\; kA(s)^m$. Note that we assume that $A$ does not depend on time – we observed that, after some time, the drainage $A$ stabilizes in the main river channels. A solution of Eqn. 3 was proposed in Earth sciences [Ste21] with the assumption that $u$ is constant in space and varies in time, which is important for geomorphologists who study the erosional response to tectonic perturbations. We prefer an orthogonal approach where the uplift varies in space but not in time, as we expect the uplift to be easier to control for the user as a function of space alone [SPF⁺23].

We use the following assumptions at the bound:

The constraint that the uplift vanishes at the boundary ensures the well-posedness of the problem, as otherwise the mountains would keep growing indefinitely.

The method of characteristic relates space $s$ and time $\backslash tau$ on a set of characteristic curves with the equation:

As we seek the elevation at time $t$ and position $x$, we only consider the curve that includes the point $x,\; t$. For this curve, the space-time relationship is expressed by $\backslash tau\_\{x,t\}(s)$ defined such as $\backslash tau\_\{x,t\}(x)\; =\; t$, which gives after integration:

Note that $\backslash tau\_\{x,t\}(s)$ parameterizes time, and therefore should remain positive. The function $\backslash tau\_\{x,t\}$ being strictly increasing, we define a point $D\_\{x,t\}$ by $\backslash tau\_\{x,t\}(D\_\{x,t\})\; =\; 0$, such that the characteristic curve is only defined when $s\; \backslash geq\; D\_\{x,t\}$. This is illustrated in the right of Figure 3 by the curve associated with $(x,\; t\_1)$. The other curve – going through $(x,\; t\_2)$ – shows a case where $\backslash tau\_\{x,t\}(D\_\{x,t\})\; =\; 0$ does not have a solution, and for which we set $D\_\{x,t\}\; =\; 0$.

Then, along the characteristic curve, Eqn. 3 is rewritten as:

which we integrate on the domain of definition of the curve $s\; \backslash geq\; D\_\{x,t\}$, to obtain the formula for the elevation:

Note that for large $t$ or small $x$, $D\_\{x,t\}$ can be negative, in which case the boundary assumptions (Eqn. 4) gives the solution:

If this happens for all cells in the terrain, this solution corresponds to the steady state of the stream power law.

Compared to previous simulation-based methods that required iterating over time, this solution directly expresses the elevation of the terrain from the uplift and drainage area. This is why we follow geology literature [Ste21] to call this an analytical solution (with respect to time), even if we need to resort to a numerical evaluation of the integral over space.

Solutions based on dimensionless variables in 1D [RTP13] evaluates the solution at dimensionless positions, with a non-trivial mapping to real positions. Instead, we first propose a numerical evaluation of the analytical solutions (Eqn. 8) on the 1D case, at the real positions along the river path. We will explain in Section 4.3 how we use this solution to model 2D heightmaps. For now, we assume that all the values ($z\_0$, $u$, $a$) are known and stored in a 1D array. We denote by $\backslash delta\; x$ the spacing between cells of the array.

Our goal is to compute all the elevations in the array, essentially fixing $t$ and varying $x$. A naive implementation of Eqn. 8 is inefficient: for each position $x$, we would need to parse the array from $D\_\{x,t\}$ to $x$ to compute the integral, resulting in a quadratic complexity. Instead, we observe that significant portions of the computation is similar among neighboring cells of the array, which means that the computation of the elevation at a position $x$ can use data already computed at the previous position $x\; -\; \backslash delta\; x$. This observation suggests an iterative algorithm, which parses the array and progressively computes the elevations in a single iteration.

We introduce the function $T(x,\; y)$:

which gives the time required for the elevation to be advected from position $x$ to $y$. The first step of our algorithm is to compute $D\_\{x,t\}$, defined implicitly by $\backslash tau\_\{x,t\}(D\_\{x,t\})\; =\; 0$, or

Eqn. 11 is not well defined for small $x$ or large $t$, which are cases where the information would need to be advected from beyond the boundaries of the domain. To detect this situation, we set $T(0,\; 0)\; =\; 0$ and parse the array from the bound, progressively evaluating $T(0,\; x)\; =\; T(0,\; x\; -\; \backslash delta\; x)\; +\; \backslash delta\; x/a(x)$ until we find the first $x\_0$ for which $T(0,\; x\_0)\; >\; t$. We set $D\_\{x,t\}\; =\; 0$ for , and, assuming that $1/a$ is locally constant near $x\; =\; 0$:

Then, we propose a recursive formulation to compute $D\_\{x,t\}$ elsewhere. For a given cell position $x\; >\; x\_0$, we assume that we know $D\_\{x-\backslash delta\; x,t\}$. Writing Eqn. 11 for both $x\; -\; \backslash delta\; x$ and $x$ results in:

We approximate the left-hand side of Eqn. 13 as $\backslash delta\; x/a(x)$. Starting from $D\_\{x-\backslash delta\; x,t\}$, we progressively parse the array until we find the highest $x\_i$ such that , which allows us to split the integral and rewrite Eqn. 13 as:

and, again using a locally constant approximation of $1/a$, we obtain:

Note that for each $x$, we parse a portion of the array to find $x\_i$ and compute $D\_\{x,t\}$. However, the sum of the size of all portions seen during the whole execution of the algorithm does not exceed $2N$, with $N$ being the size of the array, which leads to linear time complexity.

With $D\_\{x,t\}$ computed, the first term of Eqn. 8, $z\_0(D\_\{x,t\})$ comes naturally from the linear interpolation of the value of $z\_0$ at the two cell positions enclosing the position $D(x,\; t)$.

Next, we compute the second term of Eqn. 8. We rewrite Eqn. 8 as:

Where $S(x,\; y)$ represents the difference of elevation between positions $x$ and $y$ and results from the balance uplift and erosion:

We use the notation $S(x)$ as a shorthand for $S(D\_\{x,t\},\; x)$, and propose another recursive formulation to compute $S(x)$ from $S(x\; -\; \backslash delta\; x)$. The

boundary condition gives the initialization: $S(0)\; =\; 0$, and, by decomposition of the integral,

We have $S(x\; -\; \backslash delta\; x,\; x)\; \backslash approx\; \backslash delta\; x\; u(x)/a(x)$, and compute $S(D\_\{x,t\},\; x)$, $D(x\; +\; \backslash delta\; x,\; t)$ parsing the array between $D(x\; -\; \backslash delta\; x,\; t)$ and $D(x,\; t)$, progressively accumulating $\backslash delta\; x\; u/a$. Note that $D(x\; -\; \backslash delta\; x,\; t)$ and $D(x,\; t)$ do not necessarily coincide with cell boundaries, and therefore the evaluation of the integral close to these end-points needs to be scaled by the distance between them and the closest cell boundary.

Here again, the sum of the sizes of the sub-arrays that we parse at each subsequent step of the recursion is smaller than $2N$, resulting overall in a linear time algorithm.

In summary, our algorithm consists of the following steps:

We now show that our algorithm for the 1D solution of the analytical equation readily extends to the 2D solution by changing our computation domain from a line to a tree covering the terrain surface. We explained on Section 4.1 that the analytical solutions are computed along a 1D river path, that follows the steepest path on the terrain. In practice, rivers merges into an algorithmic tree structure, that we call a river tree, and we call the set of all river trees the river network or hydrology network. On a river tree, each node has a single downstream neighbor, on which we can apply directly the recursive formulations from Section 4.1. In practice, this requires ordering the computation so that it follows the nodes of the tree, which is suggested by the literature on implicit solutions to the stream power law [BW13, CBC*16] or previous attempts toward analytical solutions in the geology literature [Ste21].

The ordering is (weakly) defined by assigning to each cell $c$ of the terrain a receiver $c\_r$, which is one of the four direct neighbors to the cell, and whose elevation is strictly lower than the cell. In our case, we choose the receiver randomly, with a probability proportional to the difference in elevation between the cell and its lower neighbors. We prefer a random sampling over the deterministic choice of the lowest neighbor [BW13, CBC*16] that was designed for irregular grids and yields the formation of uniform axis-aligned valleys on regular grids. (Figure 8). Furthermore, the noise induced by the randomization results in more naturally diverse patterns. To ensure the deterministic behavior of our algorithm, we precompute and store the random variable used by each cell to sample the receiver.

This ordering results in a set of river trees whose roots correspond to the bounds of the terrain and whose leaves correspond to the ridges. The inset figure below shows the river network for a $32\; \backslash times\; 32$ terrain, different colors corresponding to different trees.

We parse this tree in two directions: first, we accumulate the area of the cells $\backslash delta\; x^2$ from the leaves to the root to compute the drainage area $A$. Here, a topological sort on the tree guarantees the correctness of the computation of $A$ from the invariant that any cell is evaluated after all upstream cells [BW13, CBC⁺16, Ste21].

Then, we parse the tree in the other direction, from the root (bound) to the leaves, to evaluate the analytical solution. The 1D recursive algorithm (Section 4.2) translates directly to that case: each cell at position $x$ has a single receiver that plays the role of the cell at $x\; -\; \backslash delta\; x$ in the recursive formula (Eqns. 15, 18). The main issue, which distinguishes us from previous work, is that our algorithm relies on operations located at the advection location $D(x,\; t)$, downstream of the evaluated cell. We tackle this issue by using another ordering for this step: with a depth-first search, we can record the sequence of cells between the bound and the current cell, and use this sequence to evaluate the integrals involving $D(x,\; t)$ in Eqns. 15 and 18.

Strategies based on ordering share a common pitfall when local minima are present on the terrain – and analytical solutions in particular are not well defined for that case. Following previous work [BW13, CBC⁺16], we rely on a depression carving algorithm [SD21] to reroute some of the cells as if breaches were carved in the terrain to allow for the water to flow out of all local minimum.

The choice of a single receiver results in an incorrect evaluation of the slope of the terrain (as $\backslash |\backslash nabla\; z\backslash |$ requires the spatial derivative from both the $x$ and $y$ directions). We compensate for this approximation by introducing a correction term: for a cell at position $x$ with a receiver at position $x\_r$, we change $a(x)$ to $a(x)\; =\; kA^m\; \backslash frac\{\backslash delta\; x\; \backslash |\backslash nabla\; z\backslash |\}\{z(x)\; -\; z(x\_r)\}$. Note that we still evaluate the spatial derivatives of $\backslash |\backslash nabla\; z\backslash |$ downstream, as the difference between the cell elevation and the lowest elevation in the $x$ and $y$ directions, respectively. The inset figure shows the difference between not using (top) and using (bottom) the slope correction, with highlights on the isolines of a terrain constructed by enforcing a constant slope upward a boundary circle. The correction removes the directional artifacts and yields the expected concentric isolines.

Note that we need an existing terrain to compute both the ordering and the aforementioned correction term, which yields an interdependency between the ordering and the analytical solution: the analytical solution requires an ordering to provide the elevation at time $t$, but the ordering requires the same elevations. We will discuss strategies to combine ordering and analytical solutions in the next section.

To solve this interdependency, [Ste21] suggested using a fixed point algorithm, that starts with $z\_0$, computes the ordering (and correction), evaluates the analytical solutions, and iterates until the elevation converges. However, this corrects the ordering slowly because only the leaves of the trees can change their connectivity at each iteration, leading to an upper bound for the number of iterations corresponding to the number of nodes of the longest river.

We observed that this algorithm converges for the steady state case (large $t$), but could sometimes oscillate, especially for small $t$, small uplift, and large discontinuities in the initial topography. We address this issue with an exponential moving average (at each iteration, we average the topography predicted by the steady state with the elevations resulting from the previous iteration.)

The large number of iterations required limits the applicability of analytical solutions to interactive applications and hinders their benefits compared to a full simulation. Therefore, we propose two strategies: one inspired by the theory of multigrid – which still performs this iterative process but with fewer iterations at different scales, and the second using an optimization to reduce the discontinuities provoked by a mismatch between the ordering and the elevations.

Multigrid is commonly used to solve elliptic PDEs with a coarse-to-fine approach, first on a coarse approximation of the problem, which is then progressively upsampled toward the final resolution. This strategy accelerates the propagation of low-frequency information, which is similar to our needs and motivates the design of a multi-scale algorithm to accelerate the convergence.

We use a mipmapping type of downsampling operator where the elevation is averaged between $2\; \backslash times\; 2$ neighboring nodes, which we apply to the initial terrain $t\_0$, the uplift $u$, and the boundary condition $b$, down to a user-specified minimal resolution. We then compute a few iterations where we subsequently order the nodes and compute analytical solutions, similarly to Section 5.1. We upsample the resulting elevation with bilinear interpolation to reiterate this process on the next scale. A simple upsampling tends to produce a regular spacing of the rivers that follows the blocks of the lower resolutions (inset, top), with sharp turns. We prevent these artifacts by adding some jittering (inset, bottom) during the upsampling process, where a uniform random offset in $[-0.25,\; 0.25]^2$ is added to the locations of the nodes before interpolation.

After experimenting with several cycles from the multigrid literature, we observed that an up/down-sample strategy inspired by the simplest V-cycle gives the best trade-off between speed and accuracy. In practice, we mostly use around 4 – 5 levels and 5 – 6 iterations, which are enough to remove the most noticeable artifacts. If artifacts remain, they are generally solved by adding more iterations, but we prefer using the optimization algorithm presented next as a postprocess as it reduces the need for tuning hyperparameters.

Our multigrid algorithm efficiently updates the river network to remove the discontinuities but is not applicable in all situations. For instance, a user might want to preserve the original drainage or use a procedurally generated river network [GGG⁺13, GBG⁺19]. For these cases, we assume that the river network is given and fixed, and we find the minimal correction to the analytical solution that removes the discontinuities. Figure 4 shows the discontinuities seen after applying the analytical solutions directly to the initial river network (left), and how the discontinuities were removed by our optimization method (right).

Let us define a non-connected neighbor $c\_n$ of a cell $c$ as one of its 4 direct neighbors such as $c$ is not a receiver of $c\_n$ and $c\_n$ is not a receiver of $c$. We denote by $\backslash text\{NCN\}(c)$ the set of non-connected neighbors of a cell $c$. We detect a discontinuity when the elevation difference between a cell $c$ and a non-connected neighbor at position $c\_n\; \backslash in\; \backslash text\{NCN\}(c)$ is larger than the elevation difference between the cell and its receiver ($c\_r$):

Reducing the discontinuity requires raising $c\_n$, lowering $c$, or some combination of both, and then propagating these changes downstream and upstream. Along a single river, multiple such situations may occur and they might also lead to conflicting changes. Therefore, we developed an optimization algorithm that is more prone to balance these different constraints, rather than a deterministic algorithm based on a parsing of the river network.

A key insight of our algorithm is that we work on the space of elevation differences $d\_c\; =\; z(c)\; -\; z(c\_r)$, which has two advantages. First, it simplifies the treatment of the non-negativity of the slope: to preserve the river network, we force $d\_c\; \backslash geq\; 0$. Second, the elevation is reconstructed by accumulating the elevation differences, this propagates the influence of each constraint along the entirety of the tree at each step of gradient descent, while an approach based on the elevations would only propagate the information by one cell at each step.

In this setting, the set of elevation differences $d\_c$ (or $\backslash mathbf\{d\}$ in vector form) is the optimization variable. During the optimization, we recompute the elevation $z\_c$ from $\backslash mathbf\{d\}$ by accumulating $d\_i$ for $i$ downstream of $c$, down to the elevation $z\_b$ of a boundary node $b$. We use the analytical solutions to compute a target elevation difference $\backslash hat\{d\}\_c$.

We build our objective function around two terms. The river term retains the elevation differences given by the analytical solution:

The discontinuity term reduces the elevation difference between a cell $c$ and a lower non-connected neighbor $c\_n$ so that $l\_c$ stays above the receiver of $c$:

We obtain $\backslash mathbf\{d\}$ by solving the following optimization problem

where $\backslash lambda\_r\; =\; 1/3$ and $\backslash lambda\_d\; =\; 2/3$ weigh the two components of the objective.

We remap $z$ to $[0,\; 1]$ and use gradient descent to solve Eqn. 22. The sequential nature of the computation of $z\_c$ from $\backslash mathbf\{d\}$ is not readily compatible with existing automatic differentiation solver, therefore we show in Appendix A how we accumulate the gradients with respect to $\backslash mathbf\{d\}$ throughout the river network. To prevent exploding gradients, we weigh each gradient of the discontinuity term $\backslash partial\; L\_d\; /\; \backslash partial\; d\_c$ by the number of nodes $N\_c$ upstream of $c$ - which can be computed from the drainage area as $N\_c\; =\; A\_c\; /\; 8\_c^2$. We use the step size for the gradient descent $l\_r\; =\; 0.01$, which we modulate locally to enforce the inequality $d\_c\; \backslash geq\; 0$.

Hillslope processes encompass the weathering of the mountain slopes and the diffusion that results from the creep flow of the eroded material. It is often simplified as a linear diffusion of the terrain elevation [BS97, SPF⁺23, CJP⁺23]:

Where $k\_h$ is the hillslope erosion coefficient.

An analytical solution including the hillslope erosion would require joining the stream power law (Eqn. 1) with Eqn. 23 which results in an advection-diffusion equation. Solving this equation needs a global 2D treatment, preventing our decomposition to a set of 1D solutions on the stream tree, and the complexity of the derivations and implementation of the solutions even in the 1D case challenges their usability in a terrain modeling framework.

Instead, we propose another modeling choice for the hillslope erosion, which is simpler to incorporate into our analytical solutions. We observe that, while a solution in the transient case seems out of reach, it is possible to propose a model that has the same behavior as the diffusion equation at steady-state.

Let us consider the 1D steady-state equilibrium between the hillslope and the uplift:

We integrate from $s$ to the ridge $r$, where the slope vanishes; and use Hack's law [Hac57] to associate the distance to the ridge with the drainage $A$: $r\; -\; s\; =\; CA(s)^h$, where $C$ is a constant in the range [1.4 – 2], and the exponent is usually set to $h\; =\; 0.6$:

Finally, we observe that Eqn. 3 converges to Eqn. 25 at steady state and close to the ridge with a simple modification to the term $a$:

Thermal erosion is commonly used in computer graphics to model the granular nature of mountains, which are particularly unstable above a critical talus slope:

where $s\_c\; =\; \backslash text\{atan\}(30^\backslash circ)$ is the critical slope. We modify Eqn. 3 to add the contribution of the thermal erosion:

with a new expression for $a(s)$ (also including hillslope):

We implement the condition $\backslash frac\{\backslash partial\; z\}\{\backslash partial\; s\}\; >\; s\_c$ within our algorithm. When we compute Eqn. 18 to obtain the elevation at a position $x$, we first assume that there is no thermal erosion. If the resulting slope is

above the critical angle, we change $u$ and $a$ to incorporate the thermal erosion with Eqns. 28 and 29, and use the new values to re-estimate the elevation $z(x)$. Note that changing $a$ also change $D\_\{x,t\}$, which slightly change the algorithm summarized at the end of Section 4.2: instead of computing first $D\_\{x,t\}$ for all $x$, then $S(x)$ for all and finally $z(x)$, we interleave them to compute together $D\_\{x,t\}$, $S(x)$, and $z(x)$ from their values at $x\; -\; \backslash delta\; x$.

Analytical solutions and simulation. The purpose of our algorithm based on analytical solutions is to quickly generate terrains that are similar to the results of a simulation of the stream power law. In Figure 5, we compare between our method (with multigrid, top left), a simulation [CBC⁺16] (top right), and previous work in geology [Ste21] (Fixed point iteration, bottom left). The comparison is performed at steady-state ($t\; =\; 4.6\; \backslash text\{\; My\}$) as the last method does not handle initial topography and, also, in order to limit the integration error of the simulation. In this example, we use a constant uplift modulated by a subtle noise on a $512\; \backslash times\; 512$ terrain with $\backslash delta\; x\; =\; 50\; \backslash text\{\; m\}$. In the simulation, we use 460 iterations with $dt\; =\; 10000$, which we found to be the maximal time-step that did not produce visible artifacts. The fixed-point algorithm required 43 iterations to converge.

The strength of analytical solutions is that they do not need to explicitly model the intermediate time states of the terrain, but this

means that the river network at convergence is one of the many possible networks that agree with the stream power law. Therefore, it is impossible to measure a one-to-one difference between all three approaches - note that this is a common issue on models involving the stream power law, which strongly depend on initial conditions and implementation choices. Instead, we rely on the hypsometry [Str52] (distribution of elevations) to highlight the similarity between the three results.

Table 2 shows the time required to reach these results for different resolutions. We keep the total extent of the terrain constant, therefore decreasing $\backslash delta\; x$ from 50 to 25 and 12m. We observe that we need to decrease the timestep proportionally to the cell size to prevent artifacts in the simulation, therefore the number of iterations increases in both the fixed point (43, 82, 156 iterations) and the simulation (230, 460, 980 iterations). In contrast, the multigrid methods only require the addition of one level of down/up-sampling, which leads to a complexity almost linear to the number of cells. We additionally show the performance of the optimization algorithm, in the worst case, which is when we adjust the elevations to the initial drainage (disabling the iterative approach). In practice, we observed that the optimization cleans all visible discontinuities after 50 iterations for all resolutions. Overall, we did not observe significant changes in performance with other erosion parameters.

GPU simulations of the stream power law [SPF*23] might in some cases be faster than our method even though they require many iterations. There are, however, caveats inherent to the GPU architecture that prevent them from being used in all cases. First, GPU simulations use an explicit time-stepping scheme, which bounds the admissible time step and can yield a prohibitive number of iterations for small $\backslash delta\; x$. Second, depressions in the topography lead to local minima that interrupt the river network. Similarly to other CPU algorithms, we use depression breaching [CBC*16, SD21] to enforce the continuity of the river across the depressions. The absence of such an algorithm on GPU implementations is particularly visible in cases where we erode without uplift - all the water is trapped within the depressions and the erosion only occurs in the vicinity of the topographic gradients.

We illustrate this problem in the inset figure, where we use a simulation without depression filling to erode an escarpment for 500ky. While the result should be similar to Figure 12 (right), the erosion remained local to the initial cliff and did not expand toward wide canyons.

Hydrology-based and physically-based. Many procedural methods targeted the instant generation of terrains, often through the use of noise functions [EMP*02]. Closer to us, Gnevaux et al. [GGG*13] proposed a procedural approach based on hydrology. They generate a river network whose formation and elevation is guided by a slope map. The elevation of areas between rivers is computed using a second slope map. Since our method, similarly, computes the elevations progressively in the order defined by the river network, we alter the $a(s)$ term of Eqn. 3 so as to enforce a small slope for main rivers ($\backslash |\backslash nabla\; z\backslash |\; =\; 0.06$ if the drainage is above $2500m^3y^\{-1\}$), and a stronger slope for the sides of the mountain ($\backslash |\backslash nabla\; z\backslash |\; =\; 0.6$), corresponding to the two slope maps of [GGG*13]. In Figure 6, we compare our method (right) with this hydrology-based approach. While the latter convincingly arranges rivers and surrounding mountains, our physically-based solutions yield more diverse patterns at all scales, self-emerging from the combination of uplift, stream power law, and hillslope [CMA*16].

Here, we explore our implementation choices and show the impact of the erosion models we introduced in Section 6.

Choice of the receiver. In Section 4.3, we suggested choosing the receiver randomly among the neighbors below each node. In Figure 8, we compare this strategy with the choice of the lowest neighbor (minimum receiver). Choosing the lowest neighbor leads to the formation of uniform axis-aligned slopes. In practice, this choice constrains the rivers to follow a straight line, reducing the possible convergence of small streams. This is the reason for both the axis-aligned appearance of the slopes and their steepness, as non-merging streams reduce the drainage area locally and therefore the erosion.

Hillslope erosion stabilizes the slopes of the mountains and acts as the main erosion process at low drainage where the stream power law becomes negligible [LD03]. We show the impact of hillslope in Figure 9, where we show the analytical solutions of the stream power law only at $\backslash delta\; x\; =\; 50\; m$ (left), compared with our modified formulation that includes the hillslope erosion (Section 6.1, right of Figure 9). Without hillslope, we observe the emergence of unrealistically sharp ridges and peaks.

Thermal erosion is commonly used in computer graphics [MKM89] to dampen the main cliffs that arise in procedurally generated terrains. Strong cliffs also occur naturally, for instance, if the outlet of a river is below an average plateau altitude. We synthesize this behavior on a mostly flat initial terrain - with some Perlin noise to initiate the river profiles - where we significantly reduce the altitude of one of the boundary nodes. We see cliffs appearing

in Figure 10 (left), but the stream power law does not let us control the steepness of the cliff, only the speed at which they propagate away from the low boundary node following the drainage patterns. In contrast, adding thermal erosion as explained in Section 6.2 allows us to adjust the critical angle, and hence the shape of the slopes (Figure 10, right).

We now illustrate the versatility of our method through several applications. Fig 1 shows a landscape initialized with Perlin noise and eroded at different user-provided times to illustrate the range of possible effects, from a short post-process that carves only the main slopes (at $t\; =\; 200ky$) to the complete formation of a new mountain range ($t\; =\; 1.6My$).

User-made terrains can be authored through a combination of procedural techniques, noise and modeling tools (i.e., extrusion, smoothing). In Figure 11, we asked a user to produce a coarse terrain ($\backslash delta\; x\; =\; 30\backslash text\{m\}$, top left) and use our method (without uplift) to add erosion details. We observe the different patterns produced by varying the erosion time ($100\backslash text\{ky\}$, $200\backslash text\{ky\}$, and $300\backslash text\{ky\}$), and notice in particular that both the incision depth and the shape of the eroded cavities are affected – shorter time caves numerous small gullies, which deepen and merge into fewer large valleys at a longer erosion time.

Uplift emerges from tectonic activity and is responsible for the formation of mountain ranges. It is therefore used as another tool to control the generation of large-scale landscapes. In Figure 7, we show that the different uplift-based controls developed in previous work [CCB*17, SPF*23] readily extend to our method. These methods use various strategies to fill an uplift map, and provide this map as input to a simulation of the stream power law. We use a user-painted uplift map (Figure 7, left), which could alternatively be generated by one of the aforementioned methods, and show how the map controls the formation of the $4.6\backslash text\{My\}$ old mountain at different scales (from left to right a $5\backslash text\{km\}$, $15\backslash text\{km\}$, and $25\backslash text\{km\}$ large mountain). We notice in particular that variations in the uplift map mainly dictate the trajectory of the main rivers at smaller scales, and also influence the local ridge-line elevations for larger ranges.

Escarpments are steep slopes that separate flat areas of different elevations. The sudden change in elevation yields a strong erosive response, that illustrates the need for the advective component of the analytical solutions. Indeed, a simpler solution that would only model a progressive reduction of the slopes would only cause a local smoothing of the cliffs, while the stream power law predicts a retreat of the cliff along the drainage pattern at a speed that depends on the river discharge [SS20]. We illustrate this behavior in

Figure 12, where we start from a procedurally generated cliff that separates two areas with uniform elevation (left). Then we apply our method without uplift and observe the cliff retreating after $200$ and $500\backslash text\{ky\}$. Note that without uplift, without deposition, and with flat boundary conditions, the retreating cliff leaves open a flat area around the main rivers.

Canyons, especially the well-known examples such as the Grand Canyon, do not exhibit the large flat-bottomed areas between the cliffs observed in Figure 12. The difference comes from the boundary conditions: our results in Figure 12 assumed a flat boundary condition, while a canyon is typically tributary to another river downstream. To implement this, we assume a constant $a\; =\; a(0)$ to extend $D\_\{x,t\}$ negatively beyond the bound and we assume a constant slope to deduce $z\_0(D\_\{x,t\})$. We illustrate the impact of this new boundary condition in Figure 13 where we set a single boundary node to which we assign this small slope boundary condition. We additionally set a source point on the other side of the terrain where we impose a strong drainage area $A$ to enforce the formation of a main river connecting these two points. Erosion at time $t\; =\; 700ky$ shows a fully formed canyon with multiple tributary rivers branching out. Note that the new boundary condition resulted in the desired V-shaped walls surrounding the main river.

The main limitation of our work is that our algorithm does not enforce the time-continuity of the results: in some cases and especially for very long geological times, changing the time can modify the placement of the valleys. This is due to the fact that we predict a possible river network together with the elevations, without considering the past evolution of this network. While this issue is already documented in the fixed point case [Ste21], our multigrid implementation can sometimes degrade it, because a change in the river network at a coarser scale can affect a large region of the terrain at the finest scale.

Another limitation is inherent to analytical formulation, which restricts the possible combination with other erosion laws. We proposed new models for hillslope and thermal erosion which have a critical impact when fluvial erosion is low but including other processes, sediment deposition, or simply a non-linear dependency to the slope in the stream power law (the exponent $n\; \backslash neq\; 1$) might prove challenging.