# Arenite: A Physics-based Sandstone Simulator

Zhanyu Yang, Aryamaan Jain, Guillaume Cordonnier, Marie-Paule Cani, Zhaopeng Wang, Bedrich Benes

► To cite this version:

Zhanyu Yang, Aryamaan Jain, Guillaume Cordonnier, Marie-Paule Cani, Zhaopeng Wang, et al. "Arenite: A Physics-based Sandstone Simulator." ACM Transactions on Graphics, 2025, 44 (4), 10.1145/3731201. hal-05168310

HAL Id: hal-05168310

https://hal.science/hal-05168310v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Arenite: A Physics-based Sandstone Simulator

ZHANYU YANG, Purdue University, USA

ARYAMAAN JAIN, Inria, Université Côte d'Azur, France

GUILLAUME CORDONNIER, Inria, Université Côte d'Azur, France

MARIE-PAULE CANI, Ecole Polytechnique, CNRS (LIX), IP Paris, France

ZHAOPENG WANG, Purdue University, USA

BEDRICH BENES, Purdue University, USA

<figure><img src="image_2.png" /><figcaption>Fig. 1. A large sandstone structure generated by Arenite. The user defines simple initial conditions of each object (layer hardness and vegetation), and the physics-based simulation generates the results.</figcaption></figure>

We introduce Arenite, a novel physics-based approach for modeling sandstone structures. The key insight of our work is that simulating a combination of stress and multi-factor erosion enables the generation of a wide variety of sandstone structures observed in nature. We isolate the key shape-forming phenomena: multi-physics fabric interlocking, wind and fluvial erosion, and particle-based deposition processes. Complex 3D structures such as arches, alcoves, hoodoos, or buttes can be achieved by creating simple 3D structures with user-painted erodable areas and vegetation and running the simulation. We demonstrate the algorithm on a wide variety of structures, and our GPU-based implementation achieves the simulation in less than 5 minutes on a desktop computer for our most complex example.

**CCS Concepts: • Computing methodologies → Physical simulation.**

**Additional Key Words and Phrases:** Physical simulation of natural phenomena

Authors' Contact Information: Zhanyu Yang, Purdue University, West Lafayette, USA, yang2334@purdue.edu; Aryamaan Jain, Inria, Université Côte d'Azur, Sophia-Antipolis, France, aryamaan.jain@inria.fr; Guillaume Cordonnier, Inria, Université Côte d'Azur, Sophia-Antipolis, France, guillaume.cordonnier@inria.fr; Marie-Paule Cani, Ecole Polytechnique, CNRS (LIX), IP Paris, France, marie-paule.cani@polytechnique.edu; Zhaopeng Wang, Purdue University, West Lafayette, USA, wang4433@purdue.edu; Bedrich Benes, Purdue University, West Lafayette, USA, bbenes@purdue.edu.

**ACM Reference Format:**

Zhanyu Yang, Aryamaan Jain, Guillaume Cordonnier, Marie-Paule Cani, Zhaopeng Wang, and Bedrich Benes. 2025. Arenite: A Physics-based Sandstone Simulator. *ACM Trans. Graph.* 44, 4 (August 2025), 13 pages. https://doi.org/10.1145/3731201

## 1 INTRODUCTION

Sandstone landscapes are breathtaking structures with dramatic cliffs, mesas, natural bridges, canyons, smaller-scale caves, arches, and pinnacles. Their striking visual beauty and somewhat surreal appearance have made them a popular setting for feature films due to their ability to create a sense of grandeur, isolation, and adventure. Iconic in Western movies (e.g., Monument Valley and Arches Canyon), they have also been frequently used in science fiction and fantasy movies to represent alien or otherworldly planets, creating an environment that feels both ancient and mysterious. Quite surprisingly, sandstone formations are, however, almost absent from digitally-generated landscapes. This is due to the difficulty, so far, in simulating the specific type of layered erosion that modeled them under the millennial action of wind, water, and climatic variations.

Sandstone is a sedimentary rock made of sand grains cemented together by a matrix of minerals such as quartz, feldspar, or calcite [Pettijohn et al. 2012]. The nature of this cement, and therefore its resistance to different types of erosion, can vary from one layer to
another [Garzanti 2019; Okada 1971]. In addition, in desert or semi-arid regions, the wind carries abrasive particles that wear down the surface of the sandstone. As was recently discovered, the impact of these eroding factors is reduced in the case of local stress, thanks to a unique process called *fabric interlocking* that makes compressed sandstone material much stronger. Differential erosion, where softer layers erode faster than harder ones, and irregular stress distribution create unique features such as pinnacles, natural bridges, and arches. Over millions of years, water flows, dissolving cemented layers, and digs steep canyons and gorges.

The key insight of this work is to use a combination of stress simulation and multi-factor erosion mechanisms to address the specificity of sandstone erosion. We carefully isolate the key shape-forming phenomena: multi-physics fabric interlocking, wind, fluvial erosion, and particle-based deposition processes. We show how their interplay can be used to simulate a wide variety of observable phenomena [Härtel et al. 2007]. One common challenge of physics-based algorithms is controlling the simulation to allow the authoring of the resulting landscapes. Our interactive authoring system allows users to guide the evolution of physically plausible landscapes by interactively painting key parameters such as erosion rates, material distribution, and stress feedback that directly affect crack formation. An example in Fig. 1 shows a 3D sandstone structure generated by Arenite. The user positions several 3D blocks, paints vegetation and erodibility on the surfaces, and the physics simulation then generates the results per each object. The scene definition takes less than a minute, and the overall simulation time was less than five minutes per sandstone structure.

Our contributions include: (1) a numerical method based on the material point method scheme, which enables accurate physical simulation of sandstone evolution, grounded in the negative feedback mechanism between stress and erosion discovered in recent geomorphology studies. (2) A novel approach to simulate flow and depression routing over 3D shapes represented by particles, surpassing the limitations of existing methods that only support 2.5D terrains. In particular, our new deposition method calculates how eroded particles are deposited in 3D environments, incorporating key findings in geomorphology. (3) A novel wind erosion mechanism that has a significant impact on sandstone erosion.

## 2 RELATED WORK

Terrain modeling (see also the surveys [Galin et al. 2019; Smelik et al. 2014]) has been introduced to computer graphics by early works that exploited the strong visual similarity of fractals to terrains [Fournier et al. 1982; Mandelbrot 1983]. It was noticed [Musgrave et al. 1989] that this similarity holds for geologically young terrains, which was later confirmed by large-scale user studies and perceptual metrics on terrains [Rajasekaran et al. 2022; Scott and Dodgson 2022]. However, most real terrains undergo substantial morphological changes due to environmental causes. These changes follow a conceptually simple pattern: some part of the terrain is captured by an agent, transported, and deposited. Many terrain modeling algorithms have been introduced, depending on the morphological agent. The key terrain-forming agents are wind, changes in temperature (thermal

weathering), and water (fluvial or hydraulic erosion). However, these agents may operate at different temporal and spatial scales.

Early works focused on two main phenomena: hydraulic erosion and thermal weathering. Musgrave et al. [1989]'s pioneer work based on cellular automata introduced hydraulic and thermal erosion on regular height fields and was extended to work on layered terrains to capture the underlying strata of materials [Benes and Forsbach 2001]. The diffusion model for hydraulic erosion was extended by considering more complex physical properties, such as a force exerted on the terrain [Anh et al. 2007; Chiba et al. 1998; Neidhold et al. 2005] and was implemented even at interactive fram erates [Štava et al. 2008]. Large-scale hydraulic and thermal erosion simulations were achieved by utilizing terrain tiling on the GPU [Vaněk et al. 2011]. Rivers are an essential terrain-forming phenomenon [Génevaux et al. 2013]. The related erosion phenomena can be captured through the concept of Fluvial erosion from geomorphology [Braun and Willett 2013; Whipple and Tucker 1999], which is more exact than hydraulic erosion, as it considers the drainage areas and more complex water flow. Special attention has been paid to the mutual effect of continental uplift and fluvial erosion that work in tandem on similar time scales [Cordonnier et al. 2016, 2018; Schott et al. 2023]. A recent work introduced “unerosion” algorithms that invert the thermal and fluvial erosion in time [Yang et al. 2024].

Other phenomena that affect terrain formation considered in computer graphics were glaciers [Cordonnier et al. 2023], flora and fauna [Cordonnier et al. 2017; Ecormier-Nocca et al. 2021].

Heightfields or layers are a suitable representation for large-scale terrains, but cannot fully capture small and medium-scale (tens or hundreds of meters) 3D features. Close to our approach are Arches [Peytavie et al. 2009], which used a layered representation [Benes and Forsbach 2001] in an interactive way to carve complex 3D geomorphological objects. While their approach is based on user interaction, ours uses a different data representation and physics. Medium-scale terrain erosion has been implemented as Eulerian [Benes et al. 2006] and Lagrangian [Krištof et al. 2009] approaches, and a recent work addressed the problem of debris captured by the moving fluid [Jain et al. 2024a]. Paris et al. [2019a] used construction trees and open shape grammars to model complex volumetric terrain structures, and Becher et al. [2019] introduced Feature Curves to efficiently generate large-scale 3D features. A full 3D physics simulation enabled dynamic medium-scale features, but it is challenging to control. Close to our work is the spherical erosion introduced by Jones et al. [2010], the wind erosion [Krs et al. 2020; Paris et al. 2019b], and the particle-based erosion in [Hartley et al. 2024]. Our approach also uses particles but uses complex stress simulation and flow routing in 3D manifolds.

Several recent approaches are based on deep learning, and they attempt to generate terrains by using models trained on real-world data [Guérin et al. 2017; Jain et al. 2024c] or, similar to superresolution, enhance sparse data with details [Argudo et al. 2018; Grenier et al. 2024; Schott et al. 2024]. The style of the terrain has been disentangled via deep neural models, and the user is provided a conditional generator and sketch-based tools for control [Huang and Yuan 2023; Liu and Benes 2025; Lochner et al. 2023]. The lack of specific, accurate 3D data makes it difficult for these models to capture the unique geological shapes observed in sandstone. This
work uses a simulation approach to model complex geomorphic phenomena on sandstone.

Sandstones have been extensively studied in geology [Adamovic et al. 2006; Cilek et al. 2007; Turkington and Paradise 2005] but with a prime focus on the qualitative understanding of the processes responsible for sandstone features and formation. In particular, recent advances highlighted the interplay between stress and erosion [Bruthans et al. 2014], stress and arcade formation [Safonov et al. 2020], the importance of wind erosion [Loope and Mason 2009; Ostanin et al. 2017] (see [Migoń 2021] for a recent survey).

We present a novel and practical solution for modeling sandstone formations for computer graphics applications. Our method relies on the 3D simulation of physical phenomena and is fast enough for interactive shape control.

# 3 OVERVIEW

Sandstone formation in CG is an understudied topic due to several specific challenges: such mid-scale objects (tens to hundreds of meters) must be fully represented in 3D, which adds complexity to the physics simulation. Moreover, sandstones undergo a complex interplay of erosion and internal stresses, requiring novel algorithms that need to operate on dedicated data representations.

## 3.1 Fabric Interlocking Controls Sandstone Erosion

Sandstones are sedimentary rocks composed of sand grains (quartz or feldspar) glued by a cementing material (e.g., silica) [Okada 1971; Pettijohn et al. 2012]. These rocks are formed by the accumulation and compaction of sand deposits at the outlet of rivers, beaches, deserts, or shallow continental shelves. After millions of years, changes in geological forces exhume the sandstones and subject them to erosion, which leads to intricate patterns, especially in arid environments where chemical weathering is minimal and vegetation does not conceal the underlying rocks.

Iconic structures stemming from sandstone erosion, such as mesas, canyons, arches, alcoves, buttes, or caves, have been widely studied in geomorphology. Yet, only recently, a consensus emerged on the key process in their formation, namely a negative feedback between erosion and internal stress [Bruthans et al. 2014]. Stress impacts erosion through a process called fabric interlocking, where high loads favor the chemical dissolution and re-precipitation of sand grains. The rearranged internal structure is then characterized by a lower porosity and higher grain contact area, which increases resistance to erosion [Migoń 2021]. A typical example is load-bearing areas such as pillars, where the local stress increases when weathering reduces the cross-section. When a critical threshold is reached, fabric interlocking transforms granular sediment into a strong, rock-like material that can resist further erosion.

## 3.2 Arenite

Our sandstone simulator Arenite (see Fig. 2) tackles the efficient simulation of the fabric interlocking process, enabling the emergence of most features observed in sandstone landscapes. Our solution relies on a hybrid data structure that captures the interplay between stress and erosion, as well as a unique way to simulate the progressive

<figure><img src="image_3.png" /><figcaption>Fig. 2. Arenite overview. The rough input shape has some areas painted by the user (viability, vegetation) and is converted into particles in a grid. Stress distribution is calculated, and erodibility is determined. The wind and fluvial erosion then modify the shape, and some material is deposited, which results in the final shape.</figcaption></figure>

change of rock state through the notions of erodibility (capacity to erode) and viability (current state), both locally updated over time.

**Stress and Erosion.** We model the interplay between stress and erosion as follows. Considering a small volume of rock at the surface of the sandstone, we estimate its viability b that describes the expectancy that it has not yet been eroded. The viability is initialized as $b(t = 0) = 1$, and evolves as:

$$ \frac{\partial b}{\partial t} = -\mathcal{E}(\sigma_c) (W + \mathcal{F}) \quad (1) $$

where W and F are scalars representing local erosion effects caused by the wind and fluvial processes.

When b is small, the local eroded rock is transported and deposited. This decreases the rock volume, which in turn increases its Cauchy stress tensor $\sigma_c$. When stress exceeds a threshold I, erodibility E switches to the fabric interlocking state:

$$ \mathcal{E}(\sigma_c) = \begin{cases} k_s & \text{if } \operatorname{trace}(\sigma_c) > I \\ k_w, & \text{otherwise.} \end{cases} \quad (2) $$

The parameters $k_w$ and $k_s$, where $k_w \gg k_s$ stand for weak versus strong material, and respectively model the erodibility of the local rock layer before and after the interlocking phenomenon.

**Hybrid Data Structure.** Modeling such negative feedback between stress and erosion and the evolution of the 3D rock shape requires a data structure capable of capturing the stress distribution within a 3D volume and its temporal evolution. We adopt a hybrid representation that combines a regular grid with particles. In contrast to similar hybrid structures in fluid simulation and continuum dynamics [Brackbill et al. 1988; Sulsky et al. 1995], we change the purpose of particles: Instead of tracking dynamic motion -which is not needed to model static geomorphological features- they capture the local evolution of weathering, and by being moved, model both rock erosion and deposition. The benefit is that topological changes can be allowed at low costs while the grid is used for efficiently computing the internal stress forces within the remaining rock structure.

## 3.3 Erosion Algorithm

The input to our erosion simulation algorithm is the initial shape of the rock structure and a set of local erodibility parameters. The
former can be built using any geometric representation compatible with subsequent 3D sampling with particles, such as voxels, manifold meshes, or implicit surfaces. In our implementation, we use a constructive solid geometry (CSG) tree made of the union and intersection of simple cubes and spheres. Note that rough initial shapes have little influence on the visual quality of results, as even a cube may erode into a realistic sandstone pillar (see Sect. 8). Due to the sedimentary nature of sandstone, we specify the local erodibility parameters through a stack of altitude values defining the limits of sediment layers and a set of values ($k_w, k_s$) in each of these layers. This vertical stacking enables us to mimic successive layers of sandstone sediments, whose composition and, therefore, resistance to erosion depend on the climatic and geological conditions during their formation.

**ALGORITHM 1:** Sandstone simulation

**Data:** A constructive solid geometry representing the initial shape.

**Result:** Output sandstone shapes.

1 Initialization

2 Sample particles;

3 Extract grid representation;

4 **forall** Frames **do**

5 Update stress distribution using MSL-MPM (Sect.4.2);

6 Update normal vectors (Sect. 5.1);

7 **Wind Process**

8 Compute wind deflation using Eq. (12);

9 Compute wind abrasion using Eq. (8);

10 **Fluvial Process**

11 Compute receivers (Sect. 6.2);

12 Reroute depressions (Sect. 6.2);

13 Compute water discharge (Sect. 6.2);

14 Compute fluvial erosion using Eq. (14);

15 Compute overall erosion using Eq. (1);

16 **Deposition**

17 Update deposition pointers (Sect. 7);

18 Deposit eroded particles;

19 **end**

To complement this input, we sample the particles uniformly within the initial shape and assign them the erodibility parameters corresponding to their respective sediment layers. The user can paint specific erodibility values directly on certain particles for more expressive control and/or perturb them with Perlin noise (see Sect. 8.3).

We simulate Eq. (1) iteratively, splitting the equation and successively evolving each term over time with forward Euler in the following steps, as outlined in Alg. 1:

(1) **Stress computation:** We rely on the Material Point Method on the grid to efficiently compute $\sigma_c$ and deduce $\varepsilon$.

(2) **Wind erosion:** We model wind erosion $W$ through two mechanisms: (1) a background effect that simulates global, cumulative erosion, and (2) a local process that simulates localized, non-uniform erosion via SPH, capturing particle interactions between wind and sandstone.

(3) **Fluvial erosion:** We introduce a 3D flow transport algorithm to compute the fluxes of water running over the sandstone and deduce the associated fluvial erosion $\mathcal{F}$ (Eq. (13)).

(4) **Deposition:** We transport all the eroded particles to a stable location to model deposition.

The user interactively controls the simulation process, allowing parameter changes. Finally, for rendering purposes, we convert the final set of particles into a mesh by first extracting normals (see Sect. 5.1) and then using Screened Poisson Reconstruction [Kazhdan and Hoppe 2013], complemented by Least Squares Subdivision Surfaces [Boyé et al. 2010].

## 4 STRESS COMPUTATION

Internal stress has a critical impact on the response of sandstone to erosion since it alters the microscopic arrangement of sand grains through fabric interlocking. We use the Moving Least Squares formulation of the Material Point Method (MLS-MPM) [Hu et al. 2018] to evaluate the Cauchy stress tensor. MPM extends to broader continuous mechanics and captures, for instance, the detachment and fall of large rock blocks (Fig. 3). MLS-MPM also simplifies stress computation, yielding higher performance and GPU parallelization.

<figure><img src="image_4.png" /><figcaption>Fig. 3. With the dynamic nature of the MPM method, Arenite simulates broader scenarios such as rock block fall when over-eroded.</figcaption></figure>

### 4.1 Mechanical Model

Even while rock formations are visually at rest, a dynamic simulation model is needed to update the internal stress after small shape changes due to erosion. We thus model sandstone as a stiff elastic material subject to conservation of mass and momentum:

$$ \frac{D\rho}{Dt} + \rho\nabla \cdot \mathbf{v} = 0 \qquad (3) $$

$$ \rho \frac{D\mathbf{v}}{Dt} = \nabla \cdot \sigma + \rho g, \qquad (4) $$

where $D/Dt$ is the material derivative, $\rho$ the density, $\mathbf{v}$ the velocity, $\sigma$ the Cauchy stress tensor, and $g$ is gravitational acceleration.

Closing this equation requires a constitutive model that relates stress to the deformation gradient $\mathbf{F} = \partial\mathcal{Z}/\partial\mathbf{X}$, where $\mathcal{Z}(\mathbf{X}, t)$ is the deformation map and $\mathbf{X}$ the material space. We use the corotated model from [Stomakhin et al. 2012], which yields [Jiang et al. 2016]:

$$ \sigma = \frac{1}{J} (2\mu(\mathbf{F} - \mathbf{R}) + \lambda(J-1)\mathbf{F}^{-\top}) \mathbf{F}^{\top}, \qquad (5) $$

where $J = \det(\mathbf{F})$, $\mathbf{R}$ is the polar decomposition of $\mathbf{F}$, where $\mathbf{R}$ is the orthogonal matrix, and $\mathbf{S}$ is the symmetric positive-definite matrix, and $\mu$ and $\lambda$ are mechanical constants related to Young's
modulus and Poisson's ratio of the sandstone. Additionally, we assume Dirichlet boundary conditions $v = 0$ at the ground-sandstone interface.

## 4.2 MPM Discretization

We follow [Hu et al. 2018] to discretize Eqs. (4) and (5) with MLS-MPM. We use our hybrid representation as the MPM grid and particle data structure. The particles serve as material points and store properties such as mass $m_p$, position $\mathbf{x}_p$, velocity $\mathbf{v}_p$ and deformation gradient $\mathbf{F}_p$. At each timestep, particles transfer their mass and velocity to the grid. The grid velocity is then updated to account for internal and external forces. The resulting velocities are mapped back to the particles and used to advect their position.

We recall the MPM algorithm and explicit the formula used for stress computation. Let us use indices $i$ for grid nodes and $p$ for particles. At each time step (denoted here by the exponent $n$):

(1) Transfer of mass and momentum from particles to the grid, using APIC [Jiang et al. 2015]:

$$m_i^n = \sum_p m_p w_{ip}^n, \\ (mv)_i^n = \sum_p m_p w_{ip}^n (v_p^n + C_p^n (\mathbf{x}_i - \mathbf{x}_p^n)),$$

where $C$ is the affine velocity matrix, initialized as 0, and $w_{ip}$ the quadratic weighting function, centered at $\mathbf{x}_i$.

(2) Estimate the stress-related forces in the grid with the Moving Least Squares formulation [Hu et al. 2018]:

$$\mathbf{f}_i^n = - \sum_p J_p^n V_p^0 w_{ip}^n M_p^{-1} \sigma_p^n (\mathbf{x}_i - \mathbf{x}_p^n),$$

where $V_p^0$ is the initial volume of the particle and $M_p$ is a weighting factor (here $M_p = \frac{1}{4}\Delta x^2$ as we chose a quadratic weighting function).

(3) Update of the momentum on the grid with timestep $\Delta t$:

$$\mathbf{v}_i^{n+1} = ((mv)_i^n + (g + \mathbf{f}_i^n) \Delta t) / m_i^n, \quad (6)$$

(4) Transfer the grid velocities to the particles, advect, and update the affine velocity:

$$\mathbf{v}_p^{n+1} = \sum_i w_{ip}^n \mathbf{v}_i^{n+1}, \\ \mathbf{x}_p^{n+1} = \mathbf{x}_p^{n+1} + \mathbf{v}_p^{n+1} \Delta t, \\ C_p^{n+1} = \sum_i M_p^{-1} w_{ip}^n \mathbf{v}_p^{n+1} (\mathbf{x}_i - \mathbf{x}_p^{n+1})^T.$$

(5) Update the deformation gradient:

$$\mathbf{F}_p^{n+1} = (\mathbf{I} + \Delta t C_p^{n+1}) \mathbf{F}_p^n,$$

and finally compute the stress per particle with Eq. 5. The stress will be used to estimate the force at the next timestep (step 2) and to estimate the particle erodibility due to fabric interlocking (Eq. 1).

We chose a quadratic kernel for MPM due to performance concerns.

# 5 WIND EROSION

Wind erosion has often been neglected in terrain modeling, with water erosion usually being the dominant factor. However, it plays a significant part in shaping arid sandstone landscapes.

While height fields can adequately capture the formations caused by fluvial erosion, wind erosion generates complex 3D shapes through abrasion and deflation. Wind abrasion [Loope and Mason 2009] corresponds to the mechanical weathering by wind-blown particles of sand or dust. As such, abrasion is non-uniform and directional and is responsible for the uneven sandstone features [Greeley and Iversen 1985; Loope and Mason 2009]. With deflation, wind acts directly and more globally on the rock surface, from which it detaches fine-grained particles [Ostanin et al. 2017]. We separate the wind erosion term in Eq. (1) into two scalar values that respectively model wind abrasion and deflation per unit time:

$$\mathcal{W} = \mathcal{W}_a + \mathcal{W}_d. \quad (7)$$

## 5.1 Wind Abrasion

Abrasion happens when wind-blown particles impact the rock surface and dislodge sandstone grains. Geomorphology studies show [Greeley 1982; Greeley and Iversen 1985] that wind abrasion depends on the size of the particles, their kinetic energy, and the angle between the incoming velocity $v$ and the normal to the rock surface $n$. Given these observations, we use:

$$\mathcal{W}_a = k_a \|v\| (-\mathbf{n} \cdot \mathbf{v})_+, \quad (8)$$

where $k_a$ is an abrasion parameter accounting for the mass and diameter of the blown particles, $v$ is their velocity approximated to that of the local wind, and $(\cdot)_+$ is clamped to positive.

To compute $v$ at a fine-grain resolution, we simulate the wind as an incompressible fluid using Smoothed Particle Hydrodynamics (SPH), using Bender and Koschier [2015] divergence-free formulation for the favorable tradeoff between the performance and precision, and a cubic spline SPH kernel for smooth rigid and fluid couplings between wind and sandstone. We define a wind domain as twice the size of the sandstone grid. We emit incoming wind particles at the boundary of this domain and assign them the user-provided global wind direction. Particles leaving the domain are deleted. Noting that the resolution of wind velocity does not need to be as high as for the eroded sandstone particles, we used a coarser sampling of about 1/10 in our implementation. Inspired by Akinci et al. [2012], sandstone particles are treated as static obstacles for rigid/fluid coupling [Akinci et al. 2012]. Wind/rock friction is modeled using a laminar artificial viscosity model [Becker et al. 2009; Müller et al. 2004].

The normal vector $n$ to the sandstone surface from Eq. (8) is computed by projecting sandstone particles on a regular grid with a resolution close to their spacing. We identify surface particles $p$ as those in non-empty cells with at least one empty neighbor and their normal vectors from the eigenvectors of the covariance matrix of their k-nearest neighbors [Hoppe et al. 1992]. Our implementation accelerates neighbor search using a kD-Tree and $k = 20$.
## 5.2 Wind Deflation

Wind directly affects the surface by progressively detaching small particles through deflation. In contrast to abrasion caused by very strong winds during episodic storms, deflation models the long-term averaging of subtle wind effects. It is assumed to come from uniform directions, which is consistent with geological findings [Ostanin et al. 2017].

We follow Ostanin et al. [2017], which uses the Mohr-Coulomb law for dry friction to provide a condition on the drag force $F_d$ required to pull a particle out of a compact sandstone structure. In particular, the component of the force in the direction normal to the sandstone surface $F_n = F_d \cdot n$ should exceed a critical value $F_n^c$, given as:

$$F_n^c = \mu_c + \mu_f \mathbf{tr}(\sigma), \qquad (9)$$

where $\mu_c$ is a coefficient related to the internal cohesion, $\mu_f$ is a coefficient for dry friction, and $\sigma$ the Cauchy stress (Sect. 4).

We are inspired by the derivation from [Ostanin et al. 2017] to compute the effect of wind deflation on the sandstone surface. First, we average the detachment of particles over long periods:

$$W_d = k_d \int_{F_n^c}^{\infty} p(F_n) dF_n, \qquad (10)$$

where $k_d$ is a constant that accounts for the particle size to eroded surface ratio, and the characteristic time for the detachment of a single particle, and $p(F_n)$ describes the probability distribution of the normal drag force. We assume a Raleigh distribution for the normal drag [Hennessey Jr 1978]:

$$p(F_n) = \frac{F_n}{\alpha^2} e^{-|F_n|^2/(2\alpha^2)}, \qquad (11)$$

where $\alpha$ is a scaling coefficient on the wind strength. Combining Eqs. 9, 10, and 11, yields:

$$W_d = k_d e^{-(\mu_c + \mu_f \mathbf{tr}(\sigma))^2/(2\alpha^2)}. \qquad (12)$$

The sum of abrasion and deflation coefficients is used for updating through Eqs (1) and (7). The viability coefficient $b$ of surface particles, the latter being detected using the method of Sect. 5.1. Particles with $b = 0$ are considered eroded and they are removed from the still sandstone particles and deposited (see Sect. 7).

# 6 FLUVIAL EROSION ON SANDSTONES

Even in the arid landscape surrounding sandstones, water occasionally surges from episodic storms and acts as an additional process for sandstone erosion. Fluvial patterns are mainly visible as narrow gullies that converge to form dendritic networks. We discuss how to compute the fluvial erosion coefficient $\mathcal{F}$ to update a sandstone particle's viability via Eq. (1). Similar to wind erosion, fluvial erosion affects only the surface particles.

## 6.1 Choice of an Erosion Model

We borrow the shear stress erosion law from geomorphology [Howard 1994], which erodes proportionally to the shear stress of water $\tau$, when it exceeds a constant critical threshold $\tau_c$:

$$\mathcal{F} = -k_t(\tau - \tau_c)_+. \qquad (13)$$

Re-expressing $\tau$ as a function of the terrain slope $S$, flow discharge $Q$, and fluvial erosion coefficients $k_t$ and $k_f$ in [Howard 1994], yields:

$$\mathcal{F} = k_t (k_f \cdot Q^{0.6} S^{0.7} - \tau_c)_+. \qquad (14)$$

Note that while Eq. (14) is similar to the Stream Power Law used to erode large-scale landscapes [Cordonnier et al. 2016], this expression has the advantage of explicitly using the critical stress threshold. This is particularly relevant for smaller-scale geomorphological features, where discharge and stress are low, as in our case.

In practice, the critical shear stress $\tau_c$ controls the distance from higher elevations at which fluvial erosion patterns emerge, as our results show (see Fig. 17).

## 6.2 Computation of the Fluvial Erosion Terms

We now detail how we evaluate the different terms in Eq. 14.

The slope $S$ of our 3D representation of sandstone is derived from the normal vector $n$ as $S = \sqrt{n_x^2 + n_y^2}/|n_z|$. The normal was previously computed in Sect. 5.1.

Discharge $Q$ is defined as the upstream integral of precipitation and water sources [Cordonnier et al. 2016]. It can be computed for a given point $p$ as the sum of the precipitation of all the points $q$, so there is a path of strictly decreasing elevations from $q$ to $p$. While this definition was proposed for 2.5D terrains, it generalizes to 3D with the constraints that all paths are either located on the surface of an object if this surface is oriented upward ($n_z > 0$, this is not a ceiling), or traverse empty space vertically (waterfalls).

Computing $Q$ is challenging as water transfers are non-local. At the geological scales of erosion, water flows over long distances during a time step. Approximating the problem as a local one (where water is only transferred between neighboring elements) [Schott et al. 2023] stabilizes if the erosion is slow, but cannot handle *de-pressions* i.e., local minima that interrupt the flow path. Instead, we follow Fastflow [Jain et al. 2024b], which presents the estimation of discharge (flow routing) along with a solution to the problem of depression. Flow paths are updated to route the flow out of the local minima as if the depression fills and outflows.

Fastflow was designed to run at the surface of a height-field terrain. To generalize it to 3D, we use a 3D regular grid aligned to the one used for the MPM stress evaluation.

Receivers. Fastflow uses a tree-like structure for flow paths, where each element has a single receiver of strictly lower elevation. To extend this concept to 3D, we first tag grid cells as empty if they do not contain particles and as surface cells if their upper neighbor is empty. We find the altitude of the highest particle in each surface cell $c$ and assign it as a proxy elevation $z_c$. This elevation provides a total ordering of the surface cells, preventing any cycle in the flow paths.

We assign receivers to 3D cells as follows: (see Fig. 4). If the cell below $c$ is empty, we choose it as a receiver. Otherwise, we randomly select one of the four horizontal neighbors of $c$ that are empty, if any. If none of the above applies but $c$ is a surface cell, we select one of these horizontal neighbors $n$ if $n$ is a surface cell and $z_n < z_c$ (randomly with a probability proportional to $z_c - z_n$ if there is more than one choice). Otherwise, no valid receiver exists. If $c$ is a surface
<figure><img src="image_5.png" /><figcaption>Fig. 4. Choice of the receiver (blue arrow) of the yellow cells. (a) All types of cells point to the cell below if empty. (b) Otherwise, cells point to a horizontal neighbor if empty. (c) Otherwise, surface cells give to a lower surface cell.</figcaption></figure>

cell, it is tagged as a local minimum; otherwise, it is tagged as solid and masked out of the subsequent steps of the algorithm.

Note that cells on the walls and ceiling of the sandstone have a specific status: they are not surface cells but have a lower or lateral receiver. Even if water cannot reach these cells, this choice allows the evacuation of sediments eroded by wind.

*Depression Routing.* The Fastflow algorithm is applied to the new set of receivers, with modifications arising from our extension to 3D. In particular, depression routing removes local minima by connecting them to external nodes. Fastflow tags a set of cells as a basin if they flow toward the same local minima and then computes a minimum spanning tree to connect these basins. This can fail in 3D if solid surfaces entirely surround a cavity. We detect this situation during the minimum spanning tree computation (a connected component does not have a valid neighbor) and remove all the nodes of this basin from the computation.

*Flow Routing.* We now have all the components to compute the discharge $Q$. We assign to each cell at the upper boundary of the domain a constant precipitation value $p$ and set the precipitation to zero everywhere else. We use the rake-compress algorithm from Fastflow to efficiently integrate precipitation along the flow paths, which yields the discharge $Q$.

# 7 DEPOSITION

Our model simulates the gravitational settling of rock and debris particles dislodged by wind and water erosion. We neglect smaller sand grains and their transport by wind, as they are carried away from our domain size. As for fluvial erosion, the challenge of deposition lies in the multiple-cell displacement of falling debris over erosion timescales. We propose an algorithm to estimate the deposited location of the particles directly without explicitly modeling sand motion. This algorithm builds upon the 3D flow routing method detailed in Sect. 6. In particular, we reuse the 3D grid that embeds the receiver information, updated to account for depression routing.

Granular materials deposit where the surface slope falls below an angle of repose, also known as the talus slope [Musgrave et al. 1989]. We define a *stable cell* as a grid cell where the average slope of surface particles is less than this critical angle. For each surface cell $c$, we compute a *deposition pointer* $p(c)$ that indicates the next closest downward stable cell. This pointer predicts where an eroded particle originating from cell $c$ would be deposited. Initially, $p(c)$ is set to the receiver of cell $c$. To efficiently determine the final deposition location, we iteratively update these pointers using a

parallel pointer jumping technique: $p(c) = p(p(c))$. This process is repeated for a maximum of $\log(n)$ iterations, where $n$ is the total number of surface cells. Fig. 5 illustrates this algorithm, converging in three iterations in this simple 2D case.

Each cell containing eroded particles uses $p$ to increase the count of the number of particles that should be deposited in each stable cell. We deposit particles at random 2D positions within the cell and at the cell proxy elevation $z$, increasing after each deposited particle to maintain a target particle density. This approach, however, could introduce staircase artifacts. To mitigate this, we move the grid by a random offset for particles at each deposition step, ensuring particles are not grouped in the same grid cell each iteration.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Deposition Algorithm: each cell points to its deposition position after three iterations of pointer jumping.</figcaption></figure>

# 8 IMPLEMENTATION, RESULTS AND VALIDATIONS

We implemented our method using Python with PyTorch, Taichi [Hu et al. 2019] and custom CUDA kernels on the GPU. We generated all results on a desktop computer with an Intel Xeon Gold CPU clocked at 2.10GHz, 128GB RAM, and an Nvidia RTX A6000 with 48GB of memory. The meshes were textured using tri-planar mapping, vegetation was generated procedurally, and the results were rendered using Houdini (and Terragen for Fig. 6).

<figure><img src="image_7.png" /><figcaption>Fig. 6. Fluvial erosion: The cave network (left) was carved by fluvial erosion from an initial underground water path and high discharge value at the source. Setting this path on the surface results in the formation of a canyon (right).</figcaption></figure>

## 8.1 Emergent Sandstone Structures Caused by Wind

We claim that the key ingredient for the emergence of sandstone structures is the negative feedback between erosion and stress-driven fabric interlocking. We validate this assertion by showing that iconic sandstone structures emerge from our simulation from simple initial conditions. For each experiment, we show a photograph of a real rock formation, the simulation's initial conditions, the algorithm's output as a set of particles, and the final rendering.

We provide the statistics and the key parameters for each experiment in Tab. 1. We estimated the overall erosion intensity for each time step and used that estimate to choose the corresponding erodibility parameters. For the other parameters, the values in Tab. 2 can uniformly be used.
<figure><img src="image_8.png" /><figcaption>Fig. 7. A real sandstone structure (column a), the input to our simulator (column b), particle representation (column c), and realistic rendering (column d). (row 1) The natural arch from Arches National Park, Utah (Creative Commons, Wikipedia). We set a rectangular cuboid as our initial on top of a distinct support condition. The rock is composed of stacked layers with varied erodibility (darker is stronger) and is subjected to a lateral wind. (row 2) A mushroom (toadstool hoodoo) from Kanab Hoodo Trail in Utah (courtesy of Thomas Jundt, Flickr) is simulated by forming a layered block of material with a hard layer on the top (shown as a darker color). The lateral wind and our simulation then lead to the final shape. (row 3) Arcades (courtesy Spanaut, Flickr). The user brushes weaker parts on the weak layer, and the erosion simulation causes the arcades to emerge. (row 4) An alcove is a large recess in sandstone. An example from Zion National Park, Utah (courtesy James St. John, Flickr) is replicated by a large block with an erodable block at its foot. A small hole at its foot initiates the erosion process. (row 5) The butte is a large vertical block eroded by fluvial erosion from the top. The Monument Valley butte from Monument Valley, Utah (courtesy John Fowler, Flickr) is replicated using three blocks and interactively marking the top one with several instabilities. The fluvial erosion then forms the resulting object.</figcaption></figure>

**Arches:** One of the most iconic sandstone features is natural arches (natural bridges) isolated or joining neighboring hills. We reproduce a natural arch from the Arches National Park in Utah in Fig. 7 row 1. We start with a rectangular cuboid posed on two disjoint supports.

We vary the rock erodibility to match the stratigraphic appearance of the target image and add a lateral wind to achieve a similar asymmetry. The feedback between erosion and stress produces a
<table><caption>Table 1. Key parameters for the simulation of real sandstone structures. The range of parameter values corresponds to different sandstone layers with varying degrees of erosion resistance.</caption><thead><tr><th>Example</th><th>Particle Number</th><th>Fabric interlocking threshold: I</th><th>Wind abrasion: k<sub>a</sub></th><th>Wind deflation: k<sub>d</sub></th><th>Fluvial shear stress threshold: τ<sub>c</sub></th><th>Fluvial erosion: k<sub>f</sub></th><th>Deposit Viability b</th></tr></thead><tbody><tr><td>(Fig.7 row 1) Arch</td><td>1M</td><td>1.5 to 3MPa</td><td>2.5e-6 to 2.5e-5 m<sup>-2</sup>s<sup>-2</sup></td><td>5e-8 to 5e-7m<sup>-2</sup>s<sup>-2</sup></td><td>5Pa</td><td>-</td><td>0</td></tr><tr><td>(Fig.7 row 2) Mushroom</td><td>1M</td><td>3 to 6MPa</td><td>2e-6 to 2e-5m<sup>-2</sup>s<sup>-2</sup></td><td>1e-8 to 1e-7m<sup>-2</sup>s<sup>-2</sup></td><td>5Pa</td><td>-</td><td>0</td></tr><tr><td>(Fig.7 row 3) Arcade</td><td>2.5M</td><td>1 to 6MPa</td><td>-</td><td>2e-7 to 8e-7m<sup>-2</sup>s<sup>-2</sup></td><td>5Pa</td><td>-</td><td>0</td></tr><tr><td>(Fig.7 row 4) Alcove</td><td>1.5M</td><td>1.5 to 3MPa</td><td>2e-6 to 4e-5m<sup>-2</sup>s<sup>-2</sup></td><td>2e-8 to 4e-7m<sup>-2</sup>s<sup>-2</sup></td><td>5Pa</td><td>1e2Pa<sup>-1</sup></td><td>0.2</td></tr><tr><td>(Fig.7 row 5) Butte</td><td>1.5M</td><td>6MPa</td><td>1e-6m<sup>-2</sup>s<sup>-2</sup></td><td>1e-8m<sup>-2</sup>s<sup>-2</sup></td><td>1Pa</td><td>4e2 to 8e2Pa<sup>-1</sup></td><td>0.1</td></tr></tbody></table>

shape that optimizes stress to support the structure, similar to the arches found in medieval architecture.

<table><caption>Table 2. General parameters.</caption><thead><tr><th>Name</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Weak rock erodability: k<sub>w</sub></td><td>1e-4</td><td>year<sup>-1</sup></td></tr><tr><td>Strong rock erodability: k<sub>s</sub></td><td>1e-6</td><td>year<sup>-1</sup></td></tr><tr><td>Lamé parameter: μ</td><td>15</td><td>GPa</td></tr><tr><td>Shear Modulus: λ</td><td>12</td><td>GPa</td></tr><tr><td>Internal cohesion: μ<sub>c</sub></td><td>1</td><td>MPa</td></tr><tr><td>Dry friction: μ<sub>f</sub></td><td>0.75</td><td></td></tr><tr><td>Discharge multiplier: k<sub>f</sub></td><td>8e-5</td><td>Pam<sup>-3</sup>s<sup>-1</sup></td></tr><tr><td>Time step: Δt</td><td>1000</td><td>year</td></tr></tbody></table>

Mushrooms (toadstool hoodoos) are another typical sandstone structure (Fig. 7 row 2), where a large block of rock lies on top of a thin column in a seemingly unstable equilibrium. We chose a hoodoo from the Kanab Hoodo Trail in Utah. We use a vertical cuboid with varied erodibility as the input for this experiment. The upper layer is set as strong with low erodibility as the top of mushrooms is usually from a different type of rock, much more resistant to erosion [Dewing et al. 2016]. We add a lateral wind as the target image shows wind-driven asymmetry. From this initial condition, we obtain a mushroom with a shape similar to the target.

Arcades (Fig. 7 row 3) are sandstone features caused by stress and erosion that resemble architectural man-made arcades [Migoń 2021]. Previous work analyzed arcades and simulated them by using finite elements. It has been found that their size and shape are determined by the redistribution stress along planar discontinuities and the resulting stress shadows [Filippi et al. 2018; Safonov et al. 2020]. Our method simulates this phenomenon by defining erodable, weaker areas painted by the user. The stress increases between these points while they are eroded into small cavities. Eventually, high stress concentrates on thin pillars separating the newly emerged arcades.

Alcoves are large empty recess areas that form arcs in sandstones (see an example in Fig. 7 row 4). We use a cuboid as our initial shape. We disable erosion at the back, left, and right faces, and increase the hardness of the top layer. The process is initialized by a small hole at the foot of the cave, and the simulation takes care of the rest.

## 8.2 Fluvial Patterns in Sandstones.

Sandstone features from Sect. 8.1 showed primarily the interplay between stress and wind erosion, with subtle effects from fluvial erosion. At large scales, water can have a more significant impact, for instance, on the formation of buttes, canyons, and cave networks.

Buttes are isolated hills with vertical walls and flat tops (see Fig. 7 row 5). The initial pillars formed by stream erosion as the area was submerged. After it drained, fluvial erosion caused the pillar to split into vertical columns. The hard top layer -or caprock- is slowly eroded by wind and rain, causing mesa-like deposition at the butte's foot. We reproduce the butte by defining a cubical block standing on two larger ones, forming a pyramidal structure. We carve a small crevice on the caprock to create an instability, which causes fluvial erosion to converge on the right side of the butte.

<figure><img src="image_9.png" /><figcaption>Fig. 8. GPU memory usage and computational time for different particle numbers and scaled grid size.</figcaption></figure>

Canyons are primarily made of sedimentary rocks, such as limestones and sandstones [Hereford et al. 1993]. Sandstones are subject to fluvial erosion for high water discharge while capable of painting large cliffs when little or no water is present because of the critical shear stress term in Eq. 14. Fig. 6 right shows an example of a canyon simulated by fluvial erosion. We input an initial water path on a rock block caped with a gentle slope upon a hard layer of rock and impose a high discharge at the water source. Fluvial erosion carves a deep canyon, which is widened by wind erosion.

Cave networks are mainly a feature of limestone but can also occur for specific combinations of permeable and impermeable sandstones [Veress 2020]. These conditions allow the resurgence of small water streams along underground natural cracks, which amplify to form tunnels and natural cavities. We paint the initial water course inside a sandstone block with layers of rocks with different erodibility, impose a strong water discharge at the river's source, and let fluvial erosion amplify the network (Fig 6 right). We enable wind deflation with a small coefficient: while the wind is generally absent from these caves, this process accounts for any other internal disturbance due to gravity, internal flooding, and earthquakes.

## 8.3 Scale, Performance and Interactive User Control

Scale and Performance. To enable real-time simulation and interactive control, our simulations employed 1-3 million particles within a 128³ grid. On average, a single frame took around 200ms, with individual components: MPM, wind erosion and 3D flow routing (fluvial erosion + deposition) requiring 20ms, 5ms and 180ms, respectively.
By increasing the particle count and grid resolution, our simulation can generate sandstone structures with greater detail (Fig. 9). This improvement, however, incurs additional computational time and GPU memory requirements, as shown in Fig. 8.

<figure><img src="image_10.png" /><figcaption>1</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>2</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>3</figcaption></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>4</figcaption></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

Fig. 9. Increasing the particle count and grid resolution leads to structures with greater detail. For the arch case, we compare 1 million particles on a 128³ grid (row 1) and 9 million particles on a 256³ grid (row 2). The higher-resolution simulation exhibits greater detail with more pronounced pillar-like formations in the weaker layers on the top, which are characteristic of sandstone structures. For the butte, comparisons between 1.5 million particles on a 128³ grid (row 3) and 27 million particles on a 320³ grid (row 4) reveal enhanced surface irregularity and a more natural appearance.

*Visualization and Interaction.* Our method allows users to interactively visualize the current state during the simulation stage and dynamically adjust the direction and speed of the global wind. At any point in the simulation, users may also consider the current shape as part of a new input, enabling its editing with the brush tool and adding new shapes. This allows for seamless editing sessions, as illustrated in Fig. 10.

The user can paint layers of *vegetation* as shown in Fig. 11. The plants adjust the erodibility and provide a shielding effect of erosion of plants [Łukasz Pawlik et al. 2023].

<figure><img src="image_18.png" /></figure>

Fig. 10. User interaction: A lightly eroded rock cube (left) is marked at three locations by the user and eroded. Then the user added two cylinders on the top and eroded again.

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

Fig. 11. Plants: The plants hold soil and protect sandstone from wind and water erosion, playing an important role in shaping channels for sandstone. New plants grow at the deposit areas at the foot of the hill.

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

Fig. 12. Comparison to [Bruthans et al. 2014] (left). The contours of both shapes show that Arenite can generate shapes similar to those formed under natural erosion.

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

Fig. 13. Volumetric erosion from Jones et al. [2010] (left) compared to ours (right). Accounting for stress is critical to achieving natural arches (top) and arcade formations (bottom).

## 8.4 Comparisons and Ablation Studies

### 8.4.1 Comparisons.
We compare our results with lab experiments from Bruthans et al. [2014], where a weak cemented conglomerate of sand was eroded by water for several months until an arch emerges.
Fig. 12 shows that Arenite can generate a similar arch under similar conditions. We also compare in Fig 13 our method with [Jones et al. 2010], which produces results similar to sandstones with erosion based on rock surface curvature. Without considering stress, this method cannot create realistic arches and arcades.

**8.4.2 Ablation Studies.** We show the impact of our different phenomena and parameters on the formation of sandstone structures.

<figure><img src="image_26.png" /></figure>

Fig. 14. Hoodoos: are formed by fluvial and wind erosion (right). Without fluvial erosion (left), the initial sandstone cannot be separated into multiple pillars. Without the wind erosion (center), fluvial erosion cannot shape the distinctive stratification.

<figure><img src="image_27.png" /></figure>

Fig. 15. Wind erosion: If only abrasion is present, the wind does not create the typical arch (left). The sole effect of deflation is that it makes the arc regular (middle). The combination of both effects results in a more natural structure.

**Erosions.** Fig. 14 shows that both wind and fluvial erosion are critical to forming hoodoos. Their pillars are separated by fluvial erosion, but their characteristic shape is directed by wind erosion. When the fluvial erosion is missing, the pillars are not separated (left). When wind erosion is absent, the pillars do not have their typical form.

Fig. 15 shows the effect of the two types of wind erosion: deflation and abrasion. Abrasion (left) cannot create arch structures, but it adds non-homogeneity. Deflation only generates a hole in the structure, but the shape is symmetrical (middle). A combination of both effects causes an asymmetric structure with a hole.

**Fabric interlocking.** The threshold *I* in Eq. 2 critically impacts the shape of the sandstone, as highlighted in Fig. 16. Low threshold values prevent erosion even for small stresses, while high values disable the stress's influence entirely, and the sandstone is eroded uniformly.

**Fluvial erosion critical stress.** Another threshold is the $\tau_c$ variable in Eq. 12 that controls the erosion rate. Fig. 17 shows that while a high threshold keeps a significant part of the shape unaffected and separates the erosion patterns, a low threshold results in joined

<figure><img src="image_28.png" /></figure>

Fig. 16. Changing erosion threshold *I* greatly affects the sandstone shape. A low value of *I* = 0.25 stops the erosion at a very early stage. With the increasing value, the bottom part of the arch structure becomes thinner compared to the top part. When *I* = +∞, the whole sandstone would be eroded simultaneously, leading to an unnatural shape.

ridges similar to large-scale version models [Cordonnier et al. 2016].

<figure><img src="image_29.png" /></figure>

Fig. 17. Fluvial erosion critical stress: initial shape (left), $\tau_c$ = 1.5 (middle), and $\tau_c$ = 0.0. Small values lead to ridges and saddle points.

**Viability of eroded particles.** We assign an initial value for the viability of the deposit lower than the initial sandstone viability to account for the lack of cementing material, which makes them more prone to erosion. We show in Fig. 18 the impact of the choice of the deposit viability. A higher viability (left) increases the stability of the deposit slopes and self-shielding of the structure against erosion. A lower value further increases the slope's spread.

<figure><img src="image_30.png" /></figure>

Fig. 18. Viability allows user to control 1) spread of deposit 2) conformity to slope and 3) self-shielding for erosion.

## 9 CONCLUSIONS AND FUTURE WORK

We introduced Arenite, a physics-based approach for modeling sandstone geological structures. It simulates stress and multifactor erosion on a hybrid 3D representation. We showed that complex 3D structures can emerge from simple initial conditions by executing the physics simulation. We have shown the generation of a wide range of sandstone structures commonly found in nature, such as
arches, hoodoos, alcoves, arcades, and buttes. Our ablation shows the effect of individual phenomena and the simulation parameters.

Our work has limitations. While larger scenes could be theoretically calculated, our implementation is limited by computational time and memory requirements of the GPU. One prominent feature in sandstones are cracks that cause large blocks to fall. We did not model the crack propagation of sandstone because of the heterogeneity and anisotropy of the material structure. In particular, cracks propagate from pre-existing microcracks and flaws, which would require the complete modeling of the past tectonic events and the formation of sandstone. A similar challenge lies in accurately depicting moisture propagation and interactions with the sandstone, which depend on the cementing material and pre-existing defects. We used the plants as a homogeneous layer that affects the underlying material. In reality, vegetation has more complex interactions with sandstone. On one hand, plant roots could penetrate the sandstone, leading to its biological erosion. On the other hand, plant roots can help maintain small grains, such as dirt, on top of the sandstone, protecting it from other forms of erosion.

The potential avenues for future work stem directly from the limitations: the modeling of crack propagation, generalized to different types of rocks, and the effect of vegetation at varying scales. Future work also includes the study of multi-scale and distributed algorithms to enable the scaling of our method to larger scenes, and its coupling with large-scale terrain generation. Finally, Arenite uses a simple wind simulation. It would be interesting to explore further and study the impact of the wind accuracy on the realism of the results.

## Acknowledgments

This work was partially supported by NRCS #NR233A750004G044 and NIFA grants #2024-67013-42449, #2023-68012-38992, #2024-67021-42879 to Benes. The findings and conclusions should not be construed to represent any agency determination or policy. This work was supported in part by the U.S. National Science Foundation under awards 2417510 and 2412928 to Benes. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect those of the National Science Foundation. This project was sponsored by the Agence Nationale de la Recherche project Invterra ANR-22-CE33-0012-01 to Cordonnier. We want to thank Václav Cílek for sharing his insights with us on this topic.

## References

* Jirí Adamovic, Radek Mikuláš, and Václav Cílek. 2006. Sandstone districts of the Bohemian Paradise: emergence of a romantic landscape. *Geolines* 21 (2006), 5.

* Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler, and Matthias Teschner. 2012. Versatile rigid-fluid coupling for incompressible SPH. *ACM Trans. on Graph*. 31, 4 (Article 62 (July 2012), 8 pages). https://doi.org/10.1145/2185520.2185558

* Nguyen Hoang Anh, Alexei Sourin, and Parimal Aswani. 2007. Physically based hydraulic erosion simulation on graphics processingunit. In *GRAPHITE '07*. ACM, New York, NY, USA, 257–264. https://doi.org/10.1145/1321261.1321308

* Oscar Argudo, Antoni Chica, and Carlos Andujar. 2018. Terrain super-resolution through aerial imagery and fully convolutional networks. In *Computer Graphics Forum*, Vol. 37. Wiley Online Library, 101–110.

* Michael Becher, Michael Krone, Guido Reina, and Thomas Ertl. 2019. Feature-Based Volumetric Terrain Generation and Decoration. *IEEE Transactions on Visualization and Computer Graphics* 25, 2 (2019), 1283–1296. https://doi.org/10.1109/TVCG.2017.2762304

* Markus Becker, Markus Ihmsen, and Matthias Teschner. 2009. Corotated SPH for deformable solids. In *Proceedings of the EG Conf. on Natural Phenomena* (Munich, Germany) (NPH'09). Eurographics Association, Gilos, DEU, 27–34.

* Jan Bender and Dan Koschier. 2015. Divergence-free smoothed particle hydrodynamics. In *Proceedings SCA* (Los Angeles, California) (SCA '15). Association for Computing Machinery, New York, NY, USA, 147–155. https://doi.org/10.1145/2786784.2786796

* Bedrich Benes and Rafael Forsbach. 2001. Layered data representation for visual simulation of terrain erosion. In *Proceedings Spring Conference on Computer Graphics*. IEEE, 80–86. https://doi.org/10.1109/SCCG.2001.945341

* Bedrich Benes, Václav Těšinský, Jan Horný, and Sanjiv K. Bhatia. 2006. Hydraulic erosion. *Computer Animation and Virtual Worlds* 17, 2 (2006), 99–108. https://doi.org/10.1002/cav.v17:2

* S. Boyé, G. Guennebaud, and C. Schlick. 2010. Least Squares Subdivision Surfaces. *Computer Graphics Forum* 29, 7 (2010), 2021–2028. https://doi.org/10.1111/j.1467-8659.2010.01788.x

* Jeremiah U Brackbill, Douglas B Kothe, and Hans M Ruppel. 1988. FLIP: a low-dissipation, particle-in-cell method for fluid flow. *Computer Physics Communications* 48, 1 (1988), 25–38.

* Jean Braun and Sean D. Willett. 2013. A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. *Geomorphology* 180-181 (2013), 170–179.

* Jirí Bruthans, Jan Soukup, Jana Vaculikova, Michal Filippi, Jana Schweigstillová, A. L. Mayo, David Mašín, Gunther Kletetschka, and Jaroslav Říhošek. 2014. Sandstone landforms shaped by negative feedback between stress and erosion. *Nature Geoscience* 7 (2014), 597–601.

* N. Chiba, K. Muraoka, and K. Fujita. 1998. An Erosion Model Based on Velocity Fields for the Visual Simulation of Mountain Scenery. *The Journal of Visualization and Computer Animation* 9 (1998), 185–194.

* V Cílek, R Williams, A Osborne, P Migoń, and R Mikuláš. 2007. The origin and development of sandstone landforms. *Sandstone landscapes. Prague: Academia* (2007), 34–43.

* Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, Adrien Peytavie, and Eric Guérin. 2016. Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion. *Computer Graphics Forum* 35, 2, 165–175. https://doi.org/10.1111/cgf.12820

* Guillaume Cordonnier, Marie-Paule Cani, Bedrich Benes, Jean Braun, and Eric Galin. 2018. Sculpting Mountains: Interactive Terrain Modeling Based on Subsurface Geology. *IEEE TVCG* 24, 5 (May 2018), 1756–1769. https://doi.org/10.1109/TVCG.2017.2689022

* Guillaume Cordonnier, Eric Galin, James Gain, Bedrich Benes, Eric Guérin, Adrien Peytavie, and Marie-Paule Cani. 2017. Authoring Landscapes by Combining Ecosystem and Terrain Erosion Simulation. *ACM Trans. on Graph*. 36, 4 (Article 134 (July 2017), 12 pages). https://doi.org/10.1145/3072959.3073667

* Guillaume Cordonnier, Guillaume Jouvet, Adrien Peytavie, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, Eric Guérin, and James Gain. 2023. Forming terrains by glacial erosion. *ACM Trans. on Graph*. 42, 4 (2023), 1–14.

* Keith Dewing, Austin Springer, Bernard Guest, and Thomas Hadlari. 2016. Geological evolution and hydrocarbon potential of the salt-cored Hoofoo Dome, Sverdrup Basin, Arctic Canada. *Marine and Petroleum Geology* 71 (2016), 134–148.

* Pierre Ecomier-Nocca, Guillaume Cordonnier, Philippe Carrez, Anne-marie Moigne, Pooran Memari, Bedrich Benes, and Marie-Paule Cani. 2021. Authoring Consistent Landscapes with Flora and Fauna. *ACM Trans. on Graph*. 40, 4 (2021). https://doi.org/10.1145/3450626.3459952

* Michal Filippi, Jirí Bruthans, Jaroslav Říhošek, Martin Slavík, Jirí Adamovič, and David Mašín. 2018. Arcades: Products of stress-controlled and discontinuity-related weathering. *Earth-Science Reviews* 180 (2018), 159–184.

* Alain Fournier, Don Fussell, and Loren Carpenter. 1982. Computer rendering of stochastic models. *Commun. ACM* 25, 6 (1982), 371–384.

* Eric Galin, Eric Guérin, Adrien Peytavie, Guillaume Cordonnier, Marie-Paule Cani, Bedrich Benes, and James Gain. 2019. A Review of Digital Terrain Modeling. *Computer Graphics Forum* 38, 2 (2019), 553–577. https://doi.org/10.1111/cgf.13657

* Eduardo Garzanti. 2019. Petrographic classification of sand and sandstone. *Earth-science reviews* 192 (2019), 545–563.

* Jean-David Génevaux, Éric Galin, Eric Guérin, Adrien Peytavie, and Bedrich Benes. 2013. Terrain generation using procedural models based on hydrology. *ACM Trans. on Graph*. 32, 4 (2013), 1–13.

* Ronald Greeley. 1982. Aeolian Modification of Planetary Surfaces. In *The Comparative Study of the Planets, Angioletta Coradini and Marcello Fulchignoni* (Eds.). Springer Netherlands, Dordrecht, 419–434.

* Ronald Greeley and James D. Iversen. 1985. *Wind as a Geological Process: On Earth, Mars, Venus and Titan*. Cambridge University Press.

* Charline Grenier, Eric Guérin, É Galin, and Basile Sauvage. 2024. Real-time Terrain Enhancement with Controlled Procedural Patterns. In *Computer Graphics Forum*, Vol. 43. Wiley Online Library, e14992.

* Éric Guérin, Julie Digne, Éric Galin, Adrien Peytavie, Christian Wolf, Bedrich Benes, and Benoît Martinez. 2017. Interactive Example-based Terrain Authoring with
Conditional Generative Adversarial Networks. *ACM Trans. on Graph*. 36, 6, Article 228 (Nov. 2017), 13 pages. https://doi.org/10.1145/3130800.3130804

Handrj Härtel, Václav Cilek, Tomáš Herben, Andrew Jackson, and Rendel Williams. 2007. *Sandstone landscapes*. Academia Praha.

Marc Hartley, Nicolas Mellado, Christophe Fiorio, and Noura Faraj. 2024. Flexible terrain erosion. *The Visual Computer* 40, 7 (01 Jul 2024), 4593–4607. https://doi.org/10.1007/s00371-024-03444-w

Joseph P Hennessy Jr. 1978. A comparison of the Weibull and Rayleigh distributions for estimating wind power potential. *Wind Engineering* (1978), 156–164.

Richard Hereford, HC Fairlye, KS Thompson, and JR Balsom. 1993. Surficial geology, geomorphology, and erosion of archeologic sites along the Colorado River, eastern Grand Canyon, Grand Canyon National Park, Arizona. *US Geological Survey, Open-File Report* 93 (1993), 517.

Hugues Hoppe, Tony DeRose, Tom Duchamp, John McDonald, and Werner Stuetze. 1992. Surface reconstruction from unorganized points. In *Proceedings of the 19th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH '92)*. Association for Computing Machinery, New York, NY, USA, 71–78.

Alan D. Howard. 1994. A detachment-limited model of drainage basin evolution. *Water Resources Research* 30, 7 (1994), 2261–2285. https://doi.org/10.1029/94WR00757

Yuanming Hu, Yu Fang, Ziheng Ge, Ziyin Qu, Yixin Zhu, Andre Pradhana, and Chen-fanfu Jiang. 2018. A moving least squares material point method with displacement discontinuity and two-way rigid body coupling. *ACM Trans. on Graph*. 37, 4, Article 150 (July 2018), 14 pages. https://doi.org/10.1145/3197517.3201293

Yuanming Hu, Tzu-Mao Li, Luke Anderson, Jonathan Ragan-Kelley, and Frédéric Durand. 2019. Taichi: a language for high-performance computation on spatially sparse data structures. *ACM Trans. on Graph*. 38, 6, Article 201 (Nov. 2019), 16 pages. https://doi.org/10.1145/3355089.3356506

You-Liang Huang and Xue-Feng Yuan. 2023. *StyleTerrain: A novel disentangled generative model for controllable high-quality procedural terrain generation*. *Computers & Graphics* 116 (2023), 373–382.

Aryamaan Jain, Bedrich Benes, and Guillaume Cordonnier. 2024a. Efficient Debris-flow Simulation for Steep Terrain Erosion. *ACM Trans. on Graph*. 43, 4, Article 58 (July 2024), 11 pages. https://doi.org/10.1145/3658213

Aryamaan Jain, Bernhard Kerbl, James Gain, Brandon Finley, and Guillaume Cordonnier. 2024b. FastFlow: GPU Acceleration of Flow and Depression Routing for Landscape Simulation. *Computer Graphics Forum* 43, 7 (2024). https://doi.org/10.1111/cgf.15243

Aryamaan Jain, Avinash Sharma, and KS Rajan. 2024c. Learning Based Infinite Terrain Generation with Level of Detailing. In *2024 International Conference on 3D Vision (3DV)*. IEEE, 1048–1058.

Chenfanfu Jiang, Craig Schroeder, Andrew Selle, Joseph Teran, and Alexey Stomakhin. 2015. The affine particle-in-cell method. *ACM Trans. on Graph*. 34, 4, Article 51 (July 2015), 10 pages. https://doi.org/10.1145/2766996

Chenfanfu Jiang, Craig Schroeder, Joseph Teran, Alexey Stomakhin, and Andrew Selle. 2016. The material point method for simulating continuum materials. In *ACM SIGGRAPH 2016 Courses* (Anaheim, California) (*SIGGRAPH '16*). Association for Computing Machinery, New York, NY, USA, Article 24, 52 pages. https://doi.org/10.1145/2897826.2927348

Michael D. Jones, McKay Farley, Joseph Butler, and Matthew Beardall. 2010. Directable Weathering of Concave Rock Using Curvature Estimation. *IEEE TVCG* 16, 1 (2010), 81–94. https://doi.org/10.1109/TVCG.2009.39

Michael Kazhdan and Hugues Hoppe. 2013. Screened poisson surface reconstruction. *ACM Trans. on Graph*. 32, 3, Article 29 (July 2013), 13 pages. https://doi.org/10.1145/2487228.2487237

Peter Krištof, Bedrich Benes, J Křivánek, and Ondřej Stáva. 2009. Hydraulic Erosion Using Smoothed Particle Hydrodynamics. *Computer Graphics Forum* 28, 2 (2009), 219–228. https://doi.org/10.1111/j.1467-8659.2009.01361.x

Vojtech Krs, Torsten Haedrich, Dominik L. Jesch, Oliver Deussen, Soeren Pirk, and Bedrich Benes. 2020. Wind Erosion: Shape Modifications by Interactive Particle-based Erosion and Deposition. In *EG/ACM SIGGRAPH SCA*. The Eurographics Association. https://doi.org/10.2312/sca.20201216

Y. Liu and B. Benes. 2025. Single-Shot Example Terrain Sketching by Graph Neural Networks. *Computer Graphics Forum* 44, 1 (2025), e15281. https://doi.org/10.1111/cgf.15281 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.15281

Joshua Lochner, James Gain, Simon Perche, Adrien Peytavie, Eric Galin, and Eric Guérin. 2023. Interactive authoring of terrain using diffusion models. In *Computer Graphics Forum*, Vol. 42. Wiley Online Library, e14941.

David B. Looper and Joseph A. Mason. 2009. Wind erosion of the wind-deposited Navajo Sandstone, USA. *Geology Today* 25, 2 (2009), 65–70. https://doi.org/10.1111/j.1365-2451.2009.00709.x

Benoît B. Mandelbrot. 1983. *The Fractal Geometry of Nature*. W.H. Freeman and Company, San Francisco.

Piotr Migón. 2021. Sandstone geomorphology – Recent advances. *Geomorphology* 373 (2021), 107484. https://doi.org/10.1016/j geomorph.2020.107484

Forest Kenton Musgrave, Craig Kolb, and Robert Mace. 1989. The synthesis and rendering of eroded fractal terrains. ACM Press, 41–50. https://doi.org/10.1145/74333.74337

Matthias Müller, Simon Schirm, Matthias Teschner, Bruno Heidelberger, and Markus Gross. 2004. Interaction of fluids with deformable solids. *Computer Animation and Virtual Worlds* 15, 3-4 (2004), 159–171. https://doi.org/10.1002/cav.18

B. Neidhold, M. Wacker, and O. Deussen. 2005. Interactive physically based Fluid and Erosion Simulation. In *Proceedings of Eurographics Workshop on Natural Phenomena*, Vol. 1. 25–32.

Hakuyu Okada. 1971. Classification of sandstone: analysis and proposal. *The Journal of Geology* 79, 5 (1971), 509–525.

Igor Ostanin, Alexander Safonov, and Ivan Oseledets. 2017. Natural Erosion of Sandstone as Shape Optimisation. *Scientific Reports* 7 (12 2017). https://doi.org/10.1038/s41598-017-17777-1

Axel Paris, Eric Galin, Adrien Peytavie, Eric Guérin, and James Gain. 2019a. Terrain Amplification with Implicit 3D Features. *ACM Trans. Graph*. 38, 5, Article 147 (Sept. 2019), 15 pages. https://doi.org/10.1145/3342765

A. Paris, A. Peytavie, E. Guérin, O. Argudo, and E. Galin. 2019b. Desertscape Simulation. *Computer Graphics Forum* 38, 7 (2019), 47–55. https://doi.org/10.1111/cgf.13815 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.13815

Francis John Pettijohn, Paul E Potter, and Raymond Siever. 2012. *Sand and sandstone*. Springer Science & Business Media.

Adrien Peytavie, Eric Galin, Jérôme Grosjean, and Stéphane Mérillou. 2009. Arches: a framework for modeling complex terrains. In *Computer graphics forum*, Vol. 28. Wiley Online Library, 457–467.

Suren Deepak Rajasekaran, Hao Kang, Martin Čadík, Eric Galin, Eric Guérin, Adrien Peytavie, Pavel Slavík, and Bedrich Benes. 2022. PTMR: Perceived Terrain Realism Metric. *ACM Trans. Appl. Percept*. 19, 2, Article 6 (July 2022), 22 pages. https://doi.org/10.1145/3514244

Alexander Safonov, Michal Filippi, David Mašín, and Jíří Bruthans. 2020. Numerical modeling of the evolution of arcades and rock pillars. *Geomorphology* 365 (2020), 107260.

Hugo Schott, Eric Galin, Eric Guérin, Axel Paris, and Adrien Peytavie. 2024. Terrain Amplification using Multi Scale Erosion. *ACM Trans. on Graph*. 43, 4 (2024), 1–12.

Hugo Schott, Axel Paris, Lucie Fournier, Eric Guérin, and Eric Galin. 2023. Large-scale terrain authoring through interactive erosion simulation. *ACM Trans. on Graph*. 42, 5 (2023), 1–15.

Joshua J Scott and Neil A Dodgson. 2022. Evaluating realism in example-based terrain synthesis. *ACM Transactions on Applied Perceptions* (TAP) 19, 3 (2022), 1–18.

Ruben Smelik, Tim Tutenel, Rafael Bidarra, and Bedrich Benes. 2014. A Survey on Procedural Modelling for Virtual Worlds. *Computer Graphics Forum* 33, 6 (2014).

Alexey Stomakhin, Russell Howes, Craig Schroeder, and Joseph M. Teran. 2012. Energetically consistent invertible elasticity. In *Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA '12)*. Eurographics Association, Goslar, DEU, 25–32.

Deborah Sulsky, Shi-Jian Zhou, and Howard L. Schreyer. 1995. Application of a particle-in-cell method to solid mechanics. *Computer Physics Communications* 87, 1 (1995), 236–252. https://doi.org/10.1016/0010-4655(94)00170-7 Particle Simulation Methods.

Alice V. Turkington and Thomas R. Paradise. 2005. Sandstone weathering: a century of research and innovation. *Geomorphology* 67, 1 (2005), 229–253. https://doi.org/10.1016/j geomorph.2004.09.028 Weathering and landscape evolution.

Juraj Vaněk, Bedrich Benes, Adam Herout, and Ondřej Štáva. 2011. Large-Scale Physics-Based Terrain Editing Using Adaptive Tiles on the GPU. *IEEE Computer Graphics and Applications* 31, 6 (2011), 35–44. https://doi.org/10.1109/MCG.2011.66

Márton Veress. 2020. Karst types and their karstification. *Journal of Earth Science* 31, 3 (2020), 621–634.

Ondřej Štáva, Bedrich Benes, Matthew Brisbin, and Jaroslav Křivánek. 2008. Interactive Terrain Modeling Using Hydraulic Erosion. In *Proceedings of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (Dublin, Ireland) (SCA'08)*. Eurographics Association, Goslar, DEU, 201–210. https://doi.org/10.5555/1632592.1632622

Kelin X. Whipple and Gregory E. Tucker. 1999. Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs. *Journal of Geophysical Research: Solid Earth* 104, B8 (1999), 17661–17674.

Zhanyu Yang, Guillaume Cordonnier, Marie-Paule Cani, Christian Perrenoud, and Bedrich Benes. 2024. Unerosion: Simulating Terrain Evolution Back in Time. *Computer Graphics Forum* 43, 8 (2024), e15182. https://doi.org/10.1111/cgf.15182

Lukasz Pawlik, Piotr Gruba, Anna Gałązka, Anna Marzec-Grządziel, Dawid Kupka, Krzysztof Szopa, Brian Buma, and Pavel Šamonil. 2023. Weathering and soil production under trees growing on sandstones – The role of tree roots in soil formation. *Science of The Total Environment* 902 (2023), 166002. https://doi.org/10.1016/j.scitotenv.2023.166002