Windblown sand around obstacles - simulation and validation of deposition patterns

Sand dunes are iconic landmarks of deserts, but can also put human infrastructures at risk, for instance by forming near buildings or roads. We present a simulator of sand erosion and deposition to predict how dunes form around and behind obstacles under wind. Inspired by both computer graphics and geo-sciences, our algorithm couples a fast wind flow simulation with physical laws of sand saltation and avalanching, which suffices to reproduce characteristic patterns of sand deposition. In particular, we validate our approach via a qualitative comparison of the erosion and deposition patterns produced by our simulator against real-world patterns measured by prior work under controlled conditions.


INTRODUCTION
Around 1 15 ℎ of the Earth's lands is covered by sand [L.Raffaele 2019], giving rise to fascinating dune landscapes, but also to unintended deposition of sand on infrastructures (buildings, roads, railways).Simulating the formation of dunes has attracted attention in computer graphics to create plausible virtual worlds [Paris et al. 2019], as well as in geo-sciences to anticipate and mitigate real-world sand deposition [Raffaele et al. 2022].On the one hand, methods from computer graphics offer high speed and user control, but make simplifying assumptions that are seldom validated against real-world experiments [Taylor and Keyser 2023].On the other hand, methods from geo-sciences achieve high accuracy by accounting for multiple physical effects, but these effects are costly to simulate [Andrea Lo Giudice and Preziosi 2020].
In this work, we draw inspiration from both communities to propose a simulator of sand erosion and deposition that focuses on saltation -the phenomenon of sand particles hoping over the terrain under the action of the wind, which dominates sand deposition patterns [Kok et al. 2012].Our method relies on an efficient fluid simulator to model the 3D+t wind field over the terrain, coupled with a sand transport model derived from geo-science laws that only acts at the surface of the terrain.While simpler than advanced simulators from the engineering literature, our approach better accounts for wind flow around obstacles than prior work in computer graphics (Figure 1).
We validate our approach by simulating windblown sand around controlled configurations of obstacles, demonstrating a qualitative agreement with reference deposition patterns that have been recently measured in the real world [Poppema, Wijnberg, et al. 2022].

RELATED WORK
Sand erosion and deposition.Sand dunes shape through complex interactions between the sand bed and the air flowing above it, referred to as aeolian processes [Nichols 2009].Intuitively, wind erodes the sand bed where its intensity is high, and it deposits the eroded particles where its intensity is low.Several modes of transportation of the eroded particles can be considered [Bagnold 1941], including reptation and creep erosion where heavy particles hope or roll on the ground over a few centimeters, suspension where light particles are carried by the wind over kilometers, and saltation where medium particles hope on the ground over distances of a few meters.Note however that there is no strict separation between these modes, and that some studies adopt a definition of saltation that includes reptation and creeping [Kok et al. 2012].Once deposited, sand particles are also subject to avalanching, a process that typically triggers when the sand bed surface reaches a critical angle.Our method focuses on saltation since it is responsible for the majority of the sand mass transported by wind over short distances.Similar phenomena of erosion and deposition also happen for snow transport [G. E. Liston 1998;Schneiderbauer and Prokop 2011;Schneiderbauer, Tschachler, et al. 2008;Vionnet et al. 2014].
Aeolian sand transport simulation.The phenomena listed above can be simulated by coupling a fluid simulator (for the wind) and a transport simulator (for the sand).Andrea Lo Giudice, Nuca, et al. [2019], Tominaga [2018] and Zhou and T. Zhang [2023] survey methods to simulate sand and snow transport based on computational fluid dynamics (CFD).Closest to our problem setting is the work by Andrea Lo Giudice and Preziosi [2020], who propose a multiphase model to accurately simulate multiple phenomena induced by the mutual interactions of sand, air, and obstacles.However, accounting for all these phenomena within a volumetric simulation results in prohibitive computation time, up to four days on a 10-processor cluster for a single simulation.In contrast, we adopt here an efficient fluid simulator and a simplified sand transport model to achieve qualitatively similar deposition patterns in a few minutes (Sec.4).
Instead of simulating macroscopic physical laws, methods based on cellular automata reproduce physical phenomena by implementing a small set of stochastic rules that dictate how sand is transported across neighboring cells of a grid discretization of the sand bed [O.Rozier 2013].While early methods only consider sand dunes [Narteau et al. 2009;Werner 1995], recent works propose additional rules to model the effect of simple obstacles [Poppema, Baas, et al. 2022], or adapt the rules to model snow transport [Kochanski et al. 2019].Since these methods model the impact of obstacles on wind and sand transport using a phenomenological point of view, each new effect requires developing dedicated rules.In contrast, we developed our method around physical laws, such that sand patterns emerge from a simple coupling between sand transport and wind dynamics.
The computer graphics community has proposed several approaches to simulate sand or snow transport with fast, albeit approximate algorithms.Early methods focus on the formation of ripples, as well as simple deposition and shadowing around obstacles [Benes and Roa 2004;Onoue and Nishita 2000;Wang and Hu 2012].More recently, Paris et al. [2019] propose a fast aeolian simulation tailored to the interactive modeling of desert landscape.Key to their approach is a simplified model of fluid flow, which consists in warping a uniform wind field to align with the local slope of the terrain, along with erosion and deposition rules akin to cellular automata.This model targets the fast generation of large desertscapes (several hundreds of meters wide), but it does not reproduce the smaller-scale wind patterns we target, such as vortices that appear behind obstacles and that are responsible for characteristic sand erosion and deposition patterns.Taylor and Keyser [2023] extended this method to model interactions with obstacles and propose a GPU implementation for real-time simulation, to which we compare in Sec. 4. We differ by adopting a generalist and more realistic, albeit more expensive, fluid simulator.Our approach also relates to the work of Feldman and O'Brien [2002], who simulate wind flow and snow transport in a 3D grid surrounding obstacles.Similarly, we simulate a volumetric wind field to capture intricate flow patterns around obstacles, but we rely on a 2D grid to model sand saltation at the surface of the terrain.
Importantly, while most prior works in computer graphics were validated by comparing the overall shape of simulated dunes with photographs of real dunes [Paris et al. 2019], we conduct a more systematic experiment where we qualitatively compare our simulations on specific obstacle configurations for which ground truth deposition patterns have been measured [Poppema, Wijnberg, et al. 2022].

SIMULATING WINDBLOWN SAND
We now describe our model to simulate sand transport and deposition under the effect of wind.While we seek for a simple model, we build it on laws derived from geo-science experiments.

Sand transport
As mentioned in Sec. 2, sand transport under wind is dominated by saltation and avalanches [Kok et al. 2012].Most previous works in geo-sciences [O.Rozier 2013] and computer graphics [Paris et al. 2019;Taylor and Keyser 2023] model sand as a height field and rely on proxies to simplify the wind dynamics at the sand surface.As an exception, Andrea Lo Giudice and Preziosi [2020] considers a fully coupled 3D transport of both the wind and the sand.While we take inspiration from this coupled model by simulating 3D fluid dynamics, we note that saltation and avalanches occur at the surface of the terrain, making these phenomena amenable to a 2.5D representation that tracks the evolution of the sand thickness ℎ over the terrain.In what follows, we allow for a non-flat ground of altitude , and define the surface altitude as  = ℎ + .
In computer graphics, transport of sand [Paris et al. 2019] or snow [Cordonnier et al. 2018] is simplified as directly linked to the wind velocity: sand is eroded by strong wind, and is deposited when the wind weakens or is shadowed by obstacles.While this approach is easily controllable and fits well with procedural wind fields, it fails to capture the complex erosion and deposition patterns induced by physically-based 3D wind flows.Instead, we propose a formulation based on a sand transport rate Q (unit: kg/m/s), which readily expresses the evolution of the sand thickness through volume conservation: ℎ where  is the sand density.Intuitively, if we consider a small volume of sand, the sand transport rate Q expresses the mass of sand that flows through the bound of the volume, and its divergence gives the net balance in sand mass per time unit.
We split the sand transport rate into two terms: where Q  accounts for the sand moved by saltation, while Q  models the sand moved by avalanching.G. Strypsteen [2021] compares several formulations of saltation rate Q  and suggests that the following law by Zingg [1952] is a good fit to measured data: where   depends on the grain size; we use   = 0.13 . 2  −4 which corresponds to a median grain size of 1.3, and u * is the shear velocity of the wind at the sand surface.The shear velocity is used in fluid mechanics to characterize the flow velocity profile at the wall boundary, and therefore represents here the capacity of the flow to entrain sediments.We note that the wind velocity is null at the ground altitude, and therefore, linearizing the velocity profile at the vicinity of the ground gives the following relationship: where  = 10 −1  is the distance to the flow surface, and we obtained the coefficient 1/50 by comparing our velocities with the velocities and shear velocities reported by G. Strypsteen [2021] from real measurements.Sand is a granular material and as such flows down when the surface slope exceeds a critical Coulomb angle   .We express this avalanching process in the avalanche rate Q  .We set Q  = 0 if ℎ = 0, otherwise we follow a rotation-invariant expression from geosciences [Lo Giudice et al. 2019]: The constant   depends on the adhesive properties of the sand but is poorly constrained [Lo Giudice et al. 2019;Andrea Lo Giudice and Preziosi 2020].Compared to the time scale of saltation, avalanches are much faster and can be considered at equilibrium.Therefore, we assign a specific timestep to this phenomenon and set     = 0.24 2 -a compromise between numerical stability and fast convergence to equilibrium.Discretization.We use a 2D staggered grid embedded in the heightmap (Fig. 2).We store the elevations at the centers of the cells (e.g., ℎ(, ) for a cell centered at position , ), while we store the transport rates and shear velocity at the interfaces between cells (e.g.,   ( + 1/2, ) or   (,  + 1/2)).The surface gradient ∇ arising in the avalanche term naturally complies with this discretization, but the norms ∥u * ∥ and ∥∇ ∥ require interpolating the vector components to obtain their values at the flux locations.Note that when the surface is not flat, the tangential velocity (and hence the shear velocity) has a vertical component.This component is discarded in the definition of the transport rate but is still required for the computation of the norm of the shear velocity and therefore stored in the center of the cells.
The cases Q  = 0 and Q  = 0 when ℎ = 0 require checking the presence of sand at the cells interface, while ℎ is stored at the center.Furthermore, a largely divergent transport rate where ℎ is small can lead to a negative sand thickness.To solve both issues, we consider independently each component of the transport rate, stored at a cell interface, e.g.,   ( + 1/2, ).We want to limit the amount of material transported by the amount of material available in the incoming cell, called the upwind cell.We determine the cell upwind to an interface from the sign of the transport rate component stored on that interface (e.g., (, ) if   ( + 1/2, ) > 0), and denote by the subscript *  a quantity stored in the upwind cell.Finally, we adjust the transport rate to the maximum amount that can be transported from the cell upwind during a time :   = ℎ       where  ∈ {1, . . ., 4} is, for each cell, the number of interfaces where the sign of the transport rates indicates that the sand leaves the cell.The final transport rate is then given component-wise ( ∈ {, }) by:

Wind flow
We assume that the wind is a Newtonian fluid and follows the incompressible Navier-Stokes Equation: u where  = 1.293  −3 is the air density and  = 1.5 10 −5  2  −1 is the air viscosity.We add the incompressibility condition We do not account for turbulence phenomena, which could be present in real life experiments, due to their complexity of implementation.Nevertheless, turbulence mainly acts as additional diffusion effects, which are inherently included in numerical diffusion.
We set the boundary conditions u = 0 at the sand-air interface and we force a constant air inflow u =  inflow on the front, left, right and top sides of the domain.We leave an open boundary condition at the back side, with Neumann boundary conditions for the velocity and Dirichlet conditions ( = 0) for the pressure.
We use Stable Fluid's [Stam 1999] Helmholtz decomposition: we first advect the velocity with an RK3 semi-lagrangian scheme, then we compute a pressure field that projects back the residual velocities into the divergence-free space.The Poisson problem arising from the pressure projections requires additional Neumann boundary conditions at the interface with the sand and the inflow domain bounds: n • ∇ = 0, where n is the normal to the obstacle or to the fluid domain.We use a sparse preconditioned conjugate gradient method for the pressure projection.Specifically, we run it with a threshold of  = 10 −6 and a maximum number of 6000 iterations as our convergence parameters.
Sand-wind interactions.We store the wind velocity and pressure in a 3D staggered grid (Fig. 2), with each three components of the velocity stored on the faces of the cell and pressure at the center of the cells.We align the 2D staggered grid of transport rates with the 3D staggered grid of wind velocities.Indeed, the transport rate for saltation depends on the shear velocity, itself depending on the wind velocity tangent to the sand surface.To compute the tangent velocity, we use a vertical linear interpolation of the wind velocity to fetch it close to the sand surface (u( + )), and we deduce the tangential velocity with u  = u − n (u • n).Note that non-flat sand surfaces lead to a  and Preziosi [2020]) and simulated by our method (b).Note that we produced this figure by running a 3D simulation and extracting a 2D slice of the sand bed and wind field, yielding more complex vortex patterns than the 2D idealized case.Our method (b) aligns well with the idealized scenario (a) by eroding sand on the left side of the dune, which faces the wind, and by depositing sand on the right side of the dune, which faces the recirculation area where the wind forms a vortex.vertical tangential velocity component, which is stored at the center of the 2D staggered grid.It might be confusing to handle this vertical component because the saltation transport rate Q  is inherently 2D and only embeds the two horizontal components, but the computation of Q  in Eq. 2 includes the norm of the 3D shear velocity which does include a vertical component.
The sand surface itself acts as a boundary condition for the wind.Encoded in a heightmap, the altitude values freely vary vertically and do not conform to a specific cell or intercell altitude.More generally, handling non-regular domains is a common problem in computer graphics, and we follow [Ng et al. 2009] who estimate the proportion  of the cell's face that intersects with the free air volume.Following a finite-volume approach, they scale the pressure gradients and the velocity by these proportions, effectively changing the equation for pressure projection to: where  is the free fluid proportion and ū is the intermediate velocity after advection.We also set the velocity to u = 0 when the associated face is entirely below the surface ( = 0).

Summary of our algorithm
In summary, our algorithm starts with an initial heightmap with ground altitude  and sand thickness ℎ, and velocities initialized at u = 0 except for the prescribed inflow values.Then, we iterate over the following steps: • Compute the proportions of the cell interface  intersecting with the fluid domain above the sand surface  = ℎ + .• Project ū in the divergence-free space by solving Eq. 8.This results in the wind velocity u.
• Fetch u close to the sand surface  (altitude  + ), deduce the tangential velocity u  and then the shear velocity u * from Eq. 3.
• Compute the saltation transport rate (Q  , Eq. 2) and the avalanche transport rate (Q  , Eq. 4), compute the component-wise maximum flux   and deduce the final transport rate Q from Eq. 5. • Evaluate erosion/deposition from mass conservation (Eq. 1) and update the sand thickness ℎ.
In practice, we let the fluid flow simulation run for 100 iterations before starting the erosion and deposition processes.

RESULTS
We validate our algorithm by reproducing sand deposition patterns documented in prior work.We start with idealized patterns for which we have reference diagrams or simulations, followed by real-world patterns for which we have reference measurements.

Idealized deposition patterns
Transverse dune.Prior work describes how transverse dunes exhibit a characteristic triangle profile that moves in the direction of the wind as sand gets eroded on the side facing the wind, and gets deposited on the side facing away from the wind, where wind recirculation occurs [Bruno and Fransos 2015;Andrea Lo Giudice and Preziosi 2020].Figure 3a illustrates this configuration schematically.Figure 3b shows that running our simulator on a triangle-shaped dune for a single time step produces a similar behavior.
Gable roof.Predicting sand deposition over human-made structures is of critical importance to preserve such structures.[L.Raffaele 2019] report typical sand deposition patterns around common structures, including gable roofs (Figure 4(a)).Our method produces similar patterns on a simple house, with sand deposition against the front and back wall, as well as over the back side of the roof (Figure 4(b)).Barchan dune.Barchans are crescent-shaped dunes that commonly appear in sand deserts all over the world [Courrech du Pont 2015].Andrea Lo Giudice and Preziosi [2020] reproduced such dunes by initializing their high-end simulator with a Gaussian-shaped pile of sand under a strong uni-directional wind.After a few time steps, the crescent shape emerges.Figure 5 reports their visualizations, along with the results of our simulator.While approximate, our fast method reproduces well the characteristic shape of the barchan.Figure 6 illustrates the result of our simulation when we let it run longer, along with the picture of a real barchan dune for visual comparison.
Star dune.[D.Zhang et al. 2012] used cellular automata to simulate the formation of star dunes.
In their experiments, star dunes emerge from a truncated conical sand pile when wind alternates over several directions.They observed that star arms form in the presence of an odd number of wind directions, the number of arms being equal to the number of directions.We reproduced their experiment for wind alternating over three and five directions, and obtained similar star-shaped patterns, as shown in Figure 7.However, because this experiment amounts to simulate wind flow over multiple, alternating expositions, it requires more iterations and computation time than our other experiments (Table 1).

Real-world patterns
We now compare our results to sand deposition patterns measured in real-world experiments.
In particular, we build on the extensive series of experiments conducted by Poppema, Wijnberg, et al. [2022], who measured sand deposition patterns around cuboid shapes placed under diverse configurations (spacing, orientation with respect to the wind) on a windy sand beach in the Netherlands.
We reproduced these configurations with our simulator and provide a visual comparison to their measurements in Figure 8.We adjusted the wind speed and viscosity of our simulator to best match their measures.While qualitative, this evaluation shows that our simulator reproduces well the erosion and deposition pattern around obstacles, as a function of the spacing between the obstacles and the inclination of the wind.
We also use measurements from Poppema, Wijnberg, et al. [2022] as a reference to compare our method to the real-time simulator of Taylor and Keyser [2023], which extends the work of Paris et al. [2019] to handle obstacles (see Figure 9).This simulator is based on stochastic rules that erode sand in front of the obstacle and deposit sand behind the obstacle.As such, this simulator cannot reproduce the complex erosion and deposition patterns observed in the real-world experiment,  where the wind flow wraps around the obstacle, creating arc-shaped sand ripples in front and on the side of the cuboid.On the contrary, our method is able to reproduce such patterns.

Impact of resolution
The spatial accuracy of our simulator depends on the resolution of the underlying heightmap and volumetric grid.

Time performance
We implemented our algorithm in Python, using PyTorch for easy GPU integration.We provide our source code at https://gitlab.inria.fr/nrosset/windblown_sand.All timings are reported from a desktop computer equipped with a single NVidia RTX 3070 GPU with 8 GB dedicated memory.
We ran experiments on domains with varying discretizations -more complex scenes requiring higher resolutions (Table 1).All simulations shown in this paper took between a few seconds to a few minutes, depending on the size of the grid and the number of iterations, which is orders of magnitude faster than the complex simulators used in engineering [Andrea Lo Giudice and Preziosi 2020].Since all the erosion/deposition updates are made on a 2D grid, the bottleneck for our time performances is the 3D fluid simulation, and in particular, the pressure projection that requires the inversion of a matrix.The number of iterations here specified is the number of times we performed this projection during a given simulation.

LIMITATIONS
In addition to missing long-range sand transport modes (such as suspension in higher altitudes), the main limitations of our simulator stem from our representation of the ground as a heightfield defined over a grid.This representation is better suited to represent axis-aligned shapes, such as the cuboid in Figure 9, than curved shapes that suffer from aliasing artifacts.Furthermore, heightfields cannot represent obstacles with overhangs, like arches and fences, while such obstacles can form effective shielding structures [Bruno, Coste, et al. 2018].

CONCLUSION AND FUTURE WORK
We introduce a simulator that combines the physical phenomena responsible for windblown sand deposition patterns.Specifically, a fluid flow simulator captures the wind dynamics over a heightfield, and drives the transportation of sand through a process called saltation.We derive the erosion and deposition of sand from the balance of saltation and of another transport mechanism, avalanching, that dominates when the slope of the sand surface exceeds a critical angle.Our experiments demonstrate that our method can reproduce reference patterns from alternative methods and field studies, and capture both the transient emergence of dunes and the static buildup of typical sand patterns around obstacles.While the simplicity of our simulator makes it faster than the ones used in engineering, it does not yet reach the speed necessary for the interactive design of obstacles for sand mitigation measures [Raffaele et al. 2022].Future work could rely on our simulator to train fast data-driven surrogate models suitable for inverse design applications [Rosset et al. 2023].

Fig. 3 .
Fig. 3. Evolution of a transverse dune, schematically (a, from Andrea Lo Giudiceand Preziosi [2020]) and simulated by our method (b).Note that we produced this figure by running a 3D simulation and extracting a 2D slice of the sand bed and wind field, yielding more complex vortex patterns than the 2D idealized case.Our method (b) aligns well with the idealized scenario (a) by eroding sand on the left side of the dune, which faces the wind, and by depositing sand on the right side of the dune, which faces the recirculation area where the wind forms a vortex.

Fig. 4 .
Fig. 4. Windfield and sand deposition over a gable roof, schematically (a, from L. Raffaele [2019]) and simulated by our method (b).Note that we produced this figure by running a 3D simulation and extracting a 2D slice of the sand bed and wind field.
Fig. 5. Formation of a barchan dune by blowing wind on a Gaussian-shaped pile of sand.Reference simulation by Andrea Lo Giudice and Preziosi [2020](a) and our simulation (b).Our method produces similar results yet employs a faster fluid solver and only considers saltation and avalanching.
Fig. 6.Visual comparison between our simulation (a) and the picture of a real barchan dune (b, Mesquite Sand Dunes in the Death Valley by Daniel Mayer, Wikimedia Commons).

Fig. 7 .
Fig. 7. Formation of star dunes by blowing wind on a truncated cone of sand, alternating over three (left) and five (right) wind directions.We refer to Fig2a) and Fig2c) from D. Zhang et al. [2012] for a comparison.

Fig. 8 .
Fig. 8.Comparison to real-world measurements byPoppema, Wijnberg, et al. [2022]  for different configurations of cuboids.Blue denotes erosion and red denotes deposition.Our method faithfully reproduces sand erosion and deposition patterns for cuboids with varying spacing under a vertical wind (a,b) and inclined wind (c,d).Note in particular how sand is eroded in front and around the cuboids, and deposited behind and between the cuboids as they are more spaced.
Fig 10 (top)  visualizes barchan dunes obtained by simulating unidirectional wind over a Gaussian pile of sand, on spatial domains of increasing resolution.The characteristic crescent shape of the dune emerges in all cases, but the ridge of the dune becomes sharper at higher resolution.Fig 10 (bottom)  shows the result of simulating a thin bed of sand under unidirectional wind around a pyramid, at different resolutions.At higher resolution, the sand deposition in front of the obstacle and in the wake of the flow is more detailed.But high resolution comes at an increased computational cost, as reported in

Fig. 9 .
Fig. 9. Comparison to the real-time method byTaylor and Keyser [2023].The measurements byPoppema,  Wijnberg, et al. [2022]  provide a reference of how sand erodes and deposits around a cuboid (a).Our simulator reproduces well this pattern (b,c), where sand level increases progressively in front of the object, then drops quickly and rises again to form a thin secondary dune, before dropping again to form a crescent-shaped cavity that wraps around the cube.In contrast, the procedural rules ofTaylor and Keyser [2023] only produce a straight dune and cavity parallel to the front face of the cube (d, reproduced fromTaylor and Keyser [2023]).Note that we used a lower initial sand thickness in (c) to better match the conditions in (d).

Fig. 10 .
Fig. 10.Rendering of the simulations at different resolutions of a barchan dune emerging from a gaussian stack of sand (a, b, c), and of a thin layer of sand blown of a pyramid (d, e, f).While the different resolutions yield a similar pattern of depositions overall, higher resolution allows to capture sharper details.

Table 1 .
Timings for the main results in this paper.