A FINITE-DIFFERENCE APPROXIMATION OF A TWO-LAYERS SYSTEM FOR GROWING SANDPILES * MAURIZIO FALCONE † AND STEFANO FINZI VITA ‡ Abstract. We study from a numerical point of view a system of partial differential equations recently proposed in [13] to model the dynamics of growing sandpiles generated by a vertical source on a flat bounded table. In such a system an eikonal equation for the standing layer of the pile is coupled to an advection equation for the rolling layer. We analyze a suitable explicit finite dif- ference scheme in one dimension and discuss its properties. A comparison is made with respect to the variational approach of [14], particularly in terms of steady-state solutions. Since the effective equilibrium configuration reached by the growing process is not completely clear for this model, we show experiments in 1D and 2D which characterize the steady-state solutions for different source terms. Key words. growing piles, granular matter, hyperbolic systems, finite difference schemes AMS subject classifications. 65M06, 65M12, 70F45 1. Introduction. In the last years several papers have been devoted to the mathematical study of granular matter, and different models have been proposed in the framework of kinetic [3, 4], cellular automata [2, 9, 17] or partial differential equations theories [6, 7, 1, 10, 14, 11, 13]. A typical prototype for granular matter is sand and one of the simplest phe- nomenon to describe is the following. Let us consider an initially empty table (a bounded open domain Ω of IR 2 ). The sand is poured on a subset of Ω according to the value of a nonnegative function f (x, t) (the source). The pile grows in height with a slope which is always lower than a characteristic value (strongly depending on the physical properties of the granular matter). This means that the sand is added to the pile at every point where f is strictly positive while the slope of the pile profile u remains below the critical value α. When the slope becomes steeper than α the sand rolls down enlarging the heap base. Physical experiments on real sandpiles show that, for f independent of t in the steady avalanching regime, the profile reaches a steady state in a finite time. The height of the pile stops growing up when at some point its base reaches the boundary of the table. After that time the sand will fall down from the table at that point (the tip moves slightly in the opposite direction), so that the pile profile remains unchanged. We will focus our attention on the numerical simulations based on differential models and in particular on a recent model proposed by Hadeler and Kuttler [13]. Their model is an improvement of the one proposed by de Gennes [10] as a simplifi- cation of the one dimensional BCRE model (see [6, 7]). All the above models use two interacting variables: the standing layer u and the rolling layer v. The interactions between u and v are described in [13] by a system of nonlinear partial differential equations for which it is quite delicate to prove existence and uniqueness results. However, a partial characterization of its stationary solutions has been recently given * Work partially supported by MIUR Project 2003 ”Metologie Numeriche Avanzate per il Calcolo Scientifico”. † Dipartimento di Matematica, Universit` a di Roma ”La Sapienza”, P.le Aldo Moro, 5 - 00185 Roma, Italy (falcone@mat.uniroma1.it). ‡ Dipartimento di Matematica, Universit` a di Roma ”La Sapienza”, P.le Aldo Moro, 5 - 00185 Roma, Italy (finzi@mat.uniroma1.it). 1