An existence result for the sandpile problem on flat tables with walls

We derive an existence result for solutions of a differential system which characterizes the equilibria of a particular model in granular matter theory, the so-called partially open table problem for growing sandpiles. Such result generalizes a recent theorem of Cannarsa and Cardaliaguet established for the totally open table problem. Here, due to the presence of walls at the boundary, the surface flow density at the equilibrium may result no more continuous nor bounded, and its explicit mathematical characterization is obtained by domain decomposition techniques. At the same time we show how these solutions can be numerically computed as stationary solutions of a dynamical two-layer model for growing sandpiles and we present the results of some simulations.


Introduction
In the last years an increasing attention has been devoted towards the study of differential models in granular matter theory (see, e.g., [2] for an overview of different theoretical approaches and models). This field of research, which is of course of strong relevance in the applications, has also been the source of many new and challenging problems in the theory of partial differential equations (see, e.g., [3,6,11,16]).
In this paper we deal with the rather simple phenomenon of the evolution of a sandpile created by pouring dry matter on a flat bounded table. In such a model the table is represented by a bounded domain Ω ⊂ R 2 , and the time-independent vertical matter source by a nonnegative function f ∈ L 1 (Ω). Recently, Hadeler and Kuttler [15] proposed a new model, extending the ones studied in [4] and [5], where the description of the heap evolution is based on the observation that granular matter forms heaps and slopes (the so-called standing layer), while small amounts of matter move down along the slopes, forming the so-called rolling layer. We also mention [16], where Prigozhin has studied, both from the theoretical and numerical points of view, a degenerate parabolic problem and its equivalent formulation as a variational inequality, and [3], where a similar approach has been used for growing sandpiles on the whole plane. It is worth to remark that the two different dynamical models of [15] and [16] (see for example [17] for a comparison between them) have theoretically the same set of admissible equilibria.
Let us denote by u(t, x) and v(t, x), t ≥ 0, x ∈ Ω, respectively the heights of the standing and rolling layers. Neglecting wind effects, the local dynamics depends on the local slope |Du| and on the local density of rolling matter v (that we shall call transport density). For stability reasons, at any equilibrium configuration the local slope |Du| cannot exceed a fixed constant (the critical slope), that we normalize to the value 1, and it must be maximal where transport occurs (that is, where v > 0). Moreover, we assume that the matter falls down the table when the base of the heap touches a portion Γ of the boundary of Ω. From the model point of view we are thus assuming that on ∂Ω \ Γ we have a (arbitrarily high) vertical wall, while on Γ the table is "open". In the following, we shall refer to the open table problem in the case Γ = ∂Ω, whereas in the partially open table problem, Γ will be a nonempty closed subset of ∂Ω.
The dynamical model proposed in [15] deals with the open table problem, but it can be extended to the partially open table problem in the following way: where ∂u ∂ν denotes the normal derivative of u. The nonlinear term which appears in the previous equations with opposite signs represents the exchange term between the two layers during the growth process. Before the equilibrium has been reached, a partial surface flow is allowed also at sub-critical slopes. As far as we know, a rigorous theory for this model is not known, and its equilibrium configurations have not been characterized in the general case. For the open table case, a finite difference scheme has been proposed in [13], which offers at the same time a tool for the numerical description of stationary solutions.
From the physical considerations above, an equilibrium configuration (u, v) for (1), with u, v nonnegative functions in Ω, must satisfy In the case of the open table problem (i.e., Γ = ∂Ω), solutions of (2) have been completely characterized by Cannarsa and Cardaliaguet [6] (see also [7] for an extension to higher dimensions). More precisely, denoting by d ∂Ω : Ω → R the distance function from the boundary of Ω, they proved that there exists a nonnegative continuous function v f (with an explicit integral representation) such that (d ∂Ω , v f ) is a solution of (2), and any other Aim of this paper is to extend these results to the case of the partially open table problem. As we shall see in the sequel, the presence of vertical walls has a relevant influence on the regularity of stationary solutions. Namely, we cannot expect the transport density v to be a continuous function as in the case of the open table problem. This fact has several consequences.
First of all, we cannot give a pointwise meaning to the boundary conditions in (2). Our choice here it to set v in L 1 (Ω), and to consider a weak formulation of the problem (see (9) below). Another possible choice, which we do not pursue here, could be to set v in the class of functions with bounded variation, so that the trace of v on ∂Ω is well defined. Our main result (see Theorem 2.2 below) states that there exists a nonnegative function v f ∈ L 1 (Ω), of which we give an explicit integral representation, such that (d Γ , v f ) is a solution to the weak formulation of the problem, where d Γ denotes the distance function from Γ (see (5) for its precise definition).
The second major consequence of the lack of continuity of v concerns the uniqueness of solutions. Namely, the uniqueness of the transport density v for the open table problem was proved in [6] using a blow-up argument, first introduced in [11] for the analysis of mass transport problems in the framework of the Monge-Kantorovich theory, that relies on the continuity of v. In our opinion this argument cannot be adapted to the case v ∈ L 1 . Nevertheless, it could be possible to prove a uniqueness result in a restricted class of more regular functions F ⊂ L 1 (Ω), provided that one can show that v f ∈ F.
A precise formulation of the existence result mentioned above requires some notation. The table Ω ⊂ R 2 will be a Lipschitz domain, i.e. an open bounded connected set with Lipschitz boundary, and the open boundary Γ ⊂ ∂Ω will be a nonempty closed subset of ∂Ω. Let d : Ω × Ω → [0, ∞) denote the path distance in Ω, defined by It is easily seen that, if Γ = ∂Ω, then d Γ is the Euclidean distance from the boundary of Ω. Moreover, d Γ belongs to the space of functions It is also known that d Γ is the maximal function among all functions in Lip 1 Γ (Ω). Thus d Γ is the maximal size of the standing layer.
Given a nonnegative function f ∈ L 1 (Ω), let u f : Ω → [0, ∞) denote the function defined by where, for every x, z ∈ Ω, and supp(f ) denotes the (essential) support of f , that is, the complement in Ω of the union of all relatively open subsets A ⊆ Ω such that f = 0 a.e. in A.
It is readily seen that the graph of the positive part of G z is, in the path metric, the maximal cone of unitary slope, with apex in z, whose base is contained in Ω and touches Γ at some point. It is clear that 0 ≤ u f ≤ d Γ , and that u f ∈ Lip 1 (Ω). Moreover, since u f is the sup-envelope of all the cones G z with z ∈ supp(f ), it is plain that u f represents the minimal standing layer for an equilibrium configuration with respect to the given support of the source.
We can now state the weak formulation of problem (2): Find nonnegative functions u, v : Ω → [0, +∞) such that In this weak formulation the boundary conditions on u and v are embedded in the choice of the test functions space C ∞ c (R 2 \ Γ) and in the condition u ∈ Lip 1 Γ (Ω). It is readily seen that, if u and v are smooth enough, then problem (9) is equivalent to (2).
The plan of the paper is the following. In Section 2 we state the main existence result, constructing explicitly a transport density v f (see formula (16) below) and showing that a pair (u, v f ) is a solution to problem (9) for every u ∈ Lip 1 Γ (Ω) satisfying u f ≤ u ≤ d Γ . Moreover, we show that no other function u ∈ Lip 1 Γ (Ω) can be the standing layer for the problem. Section 3 contains the proofs of these results, that are mainly based on a Change of Variables formula of some independent interest (see Theorem 3.3 below). In Section 4 we compute the explicit solution in a simple case, which will be compared in Section 5 to the numerical equilibrium solution of the dynamical model (1) obtained via some finite difference schemes. Indeed, in this last section the specific difficulties of such a numerical approximation are discussed in details.

Existence of a solution
Throughout this paper, Ω ⊂ R 2 and Γ ⊂ ∂Ω satisfy the following assumptions.

a nonempty open bounded connected set with
Lipschitz boundary. (H2) Γ = N i=1 Γ i is a nonempty closed subset of ∂Ω, with Γ 1 , . . . , Γ N connected arcs of ∂Ω, pairwise disjoint (up to the endpoints) and of class C 2 . We denote by A i , B i the endpoints of the arc Γ i , i = 1, . . . , N , and by Γ e the collection of all these endpoints.
Remark 2.1. Observe that (H2) does not prevent the intersection of two arcs at the endpoints. This is the case, for example, if Ω is a square and Γ i , i = 1, . . . , 4, its sides. Nevertheless, in order to simplify the proofs, in the following we assume that the arcs Γ i do not intersect at the endpoints. The general case can be treated with minor modifications.
Condition (H3) says that, for every x ∈ Ω, there is a point y ∈ Γ such that the closed segment [[x, y]] with endpoints x and y is a path of minimal length joining x to Γ.
There are three relevant cases where these conditions are satisfied: (a) Γ = ∂Ω, and Ω is a domain of class C 2 .
Case (a) and (b) refer to the open table problem. They were considered respectively in [6] (see also [7] for an extension to Ω ⊂ R n ) and [14] in the case of piecewise C 2,1 boundary with outer (i.e. "convex") corners.
For every x ∈ Ω we denote by the set of all projections of x on Γ. By (H3) it is clear that Open boundary Γ (dotted), extended ridge R (bold) and transport rays For every x ∈ Ω and every y ∈ Π Γ (x), let us define where ]]x, y[[ denotes the segment joining x to y without the endpoints. The set R x is called a distance ray (or transport ray) through x and, in general, it depends on the projection point y. On the other hand, it is not difficult to prove that m(x) and l(x) do not depend on the choice of y ∈ Π Γ (x). It is readily seen that l(x) is the length of the transport ray R x . The set will be called the extended ridge of Ω (see Figure 1). It coincides with the usual definition of ridge when Γ = ∂Ω (see [12]). Finally let τ : Ω → [0, ∞) denote the normal distance to the ridge, defined by This function is clearly bounded from above by max{d Γ (x); x ∈ Ω}.
Let Ω * ⊂ Ω be the set of regular points of d Γ , that is the set of all x ∈ Ω such that the projection Π Γ (x) is a singleton. It is well-known that L 2 (Ω \ Ω * ) = 0. Let x ∈ Ω * , let Π Γ (x) = {y}, and define for every t ∈ [0, τ (x)] Let us define the function We In other words, the transport density v f vanishes at the end of each transport ray.
The main theoretical contribution of this paper is the following existence result. Let Ω and Γ satisfy (H1)-(H3). Then the pair (u, v f ) is a solution of (9) for every u ∈ Lip 1 It is apparent that the result above characterizes all the possible standing layers as the On the other hand, the uniqueness of the transport density v f remains an open problem.

Proof of Theorem 2.2
In what follows we shall always use the following decomposition of the set Ω * of regular points of d Γ . For every i = 1, . . . , N , let us define the sets Then (see Fig. 2 for an example) Next, we need an approximation of this decomposition in terms of sets that can be easily parameterized (see Theorem 3.2 below). For every > 0 let us define the sets From Remark 2.1 and (H3), there exists r > 0 with the following property: for every ∈ (0, r], the sets G i are pairwise disjoint and of class C 1,1 . Given ∈ (0, r] and i = 1, . . . , N , for every y ∈ Γ i let ν (y) denote the outer (with respect to G i ) unit normal vector to Γ i at y. In the following, will always denote a number in (0, r]. Let us define the maps The following lemma states that Ω can be parameterized using the maps Φ i . Moreover, this parametrization is independent of the choice of . Let us define and denote by Π = Π Γ the projection operator on Γ . Γ ε ε Ω ε Γ Figure 3. The set Ω Lemma 3.1. Let x ∈ Ω and let y ∈ Π Γ (x). If y ∈ Γ i , then for every ∈ (0, r) there exists a unique y ∈ Γ i such that y ∈ Π (x) and Proof. Let y be the intersection of Γ i with the ray starting from y and passing through x. From the definition of Γ i it is straightforward to check that y satisfies all the stated properties.
Let y ∈ Γ , and let y ∈ Γ be the unique projection of y on Γ. Let us define It can be checked that κ (y ) is the curvature of Γ at every point y where the curvature is defined. We remark that, for every i = 1, . . . , N , κ is defined at all but four points of Γ i . Let x ∈ Ω * , let Π (x) = {y }, and define It is easy to check that, if Π Γ (x) = {y}, then The main tool needed for the proof of Theorem 2.2 is a Change of Variables formula (see Theorem 3.3 below). As a first step, we need the following approximation result (see Figure 3).
Proof. It is not difficult to see that the set covers Ω up to a set of Lebesgue measure N π 2 , i.e. L 2 (Ω \ O ) < N π 2 . Moreover, from Lemma 3.1 it is clear that Ω ⊂ O ⊂ Ω. The conclusion now follows observing that {Φ i (y, l(y)); y ∈ Γ i } ⊆ R , and that D has vanishing Lebesgue measure (see [9], Corollary 6.8).

Theorem 3.3 (Change of Variables). For every
Proof. From Theorem 3.2 we have that For every i = 1, . . . , N we have the decomposition defined as in (17) replacing Ω * with Ω .
Using the arguments developed in [9, Sect. 7], it can be checked that We shall now prove that (25)

Let us define the region
It is clear that A ⊂ Ω ,A i . Hence, it is enough to prove that Observe that, by the very definition of κ , we have that κ (y) = −1/ for every y ∈ Γ i ∩Ω A i . Then, the integral in brackets becomes On the other hand, the integral over A can be computed using polar coordinates (ρ, θ). We remark that Γ i ∩ Ω ,A i is an arc of circumference of radius , so that dH 1 (y) = dθ, hence formula (26) follows.
It is clear that (25) holds if Ω ,A i is replaced by Ω ,B i . The identity (23) then follows from (24) and (25).
Lemma 3.4. The function u f defined in (7) Proof. From the very definition of u f it is plain that u f (x) ≥ d Γ (x) for every x ∈ supp(f ). On the other hand, by the maximality of d Γ among all functions of Lip 1 (Ω) vanishing on Γ, we conclude that   Let (u, v) be a solution of (9). Then Proof. By a standard approximation argument we have that Let w ∈ Lip 1 Γ (Ω). Since v ≥ 0, |Dw| ≤ 1, |Du| ≤ 1, and |Du| = 1 a.e. in {v > 0}, we have that

Hence we infer that
where the last equality follows from (29) with φ = w − u. Since this inequality holds for every w ∈ Lip 1 Γ (Ω), (28) follows. Proof of Theorem 2.2. We divide the proof into three steps. In the first one we shall prove that (d Γ , v f ) is a solution to (9). Then, in the second step, we shall prove that also (u, v f ) is a solution for every u ∈ Lip 1 (Ω) satisfying u f ≤ u ≤ d Γ . Finally, in the last step we shall prove that if u ∈ Lip 1 Γ (Ω) is a nonnegative function with {u < u f } = ∅, then for every nonnegative function v ∈ L 1 (Ω) the pair (u, v) is not a solution of (9).
Step 1. We give only a sketch of the proof, since it follows the lines of the proof of Theorem 7.2 in [9].
Given φ ∈ C ∞ c (R 2 \ Γ), we have to prove that (9) holds with u = d Γ . Using the Change of Variables formula (23) we have that For every fixed y ∈ Γ , let us integrate by parts the term in brackets. Recalling that Φ (y, 0) = y − ν (y) ∈ Γ, we have that φ(Φ (y, 0)) = 0, hence the integration by parts gives

From the definition (16) of v f and (22) we deduce that
so that Observe now that Dd Γ (Φ (y, t)) = ν (y) for every y ∈ Γ and every t ∈ (0, l(y)). Substituting the last expression for I(y) in (30) and using again the Change of Variables formula (23) we finally get By the way, choosing φ = d Γ (see (29)) we have that hence the nonnegative function v f belongs to L 1 (Ω).
where the last equality follows from Step 1.
Step 3. Let u ∈ Lip 1 Γ (Ω) be a nonnegative function with {u < u f } = ∅. From [10], Proposition 4.4(iii), we deduce that also the set {x ∈ supp(f ); u(x) < d Γ (x)} is not empty and, in particular, u < d Γ on a set of positive measure where f > 0. Since u ≤ d Γ and f ≥ 0, we infer that On the other hand (u, v) is a solution to (9), so that (28) holds, in contradiction with (31).

A test example
In this section we describe a simple example which illustrates very well how the presence of vertical walls on the boundary can influence the regularity of solutions (u, v) of (9). Let Ω = (0, 1) 2 be the unit square of R 2 , and Γ = {(x, y) ∈ R 2 : 0 ≤ x ≤ 0.5 ; y = 0} the only open part of its boundary. Assume f ≡ 1 in all of Ω. From the picture in Fig. 4 we see that the sand transport rays behave differently in the two half sides of the table: in the left-hand side they lay parallel in the direction of Γ, whereas in the right-hand side they converge all together into the extremal point P = (0.5, 0), creating a singularity.
Since f = 1 in Ω, we have that u f = d Γ so that the only possible standing layer is u = d Γ . The explicit computation for the solution (d Γ , v f ) of Theorem 2.2 can be done by decomposition of the domain Ω along the segment P Q. Using polar coordinates (r, θ) S X Γ Q R P Figure 4. The domain Ω of the example: P is a singular boundary point, the solution v results discontinuous along the line P Q centered in P (with θ ∈ [0, π 2 ]) in the right hand side, from (15) and (16) we get in particular where l(θ) denotes the length of the transport ray from P to the ridge on the wall boundary along the θ direction (see (10)). It results that v f ∈ L 1 (Ω) is unbounded near P , it is discontinuous along the segment P Q, and its gradient is discontinuous along the segment P S. The graph of the functions d Γ , v f and their level lines are shown in Fig. 5.

Numerical detection of stationary solutions
In [13] the numerical approximation of the two-layer model of [15] was studied to simulate growing sandpiles on an open flat table. Here we have considered the natural generalization (1) of such a model in the case of the partially open table problem, in order to get solutions of (9) as equilibrium solutions of a system of two evolutive partial differential equations. The extension of the finite difference scheme introduced in [13] to such a system is enough straightforward. For a given discretization step h = ∆x, we introduce in the domain Ω (for simplicity, a rectangle) a uniform grid of nodes x i,j , and we denote as usual by (u n i,j , v n i,j ) the components of the discrete solutions at time t n = n∆t. Then our fully explicit finite difference scheme can be written as where the discrete gradient vectors Du n and Dv n are computed respectively, component by component, through the maxmod and the upwind finite difference operators, and D 2 u denotes the standard five-points discretization of the Laplace operator on the grid (see [13] for the details). What is new in this scheme is the wall boundary condition (37),  Figure 5. Exact solution (d Γ , v f ) of the test example with level lines whose implementation requires some comments. The standard way is the following: after (33) and (34) have been applied, we look for the sign of (Du n+1 · ν) at the wall nodes. If it is strictly positive (as it happens for nodes which are in the extended ridge R, that is which are starting points of a transport ray to Γ) then v n+1 i,j is set to zero. If this is not the case, one should modify u n+1 i,j on the boundary in order to fulfill (37). This is not the best strategy. In fact this situation corresponds to the pathological case of nodes belonging to boundary transport rays, that is when there exist straight portions of the wall boundary, as in the test example of Section 4. Referring to Fig. 4, on the west side of the square the sand flow is parallel to the boundary (and also to the mesh in this particular case) and there is no need to impose any boundary condition: the discrete solution u n naturally satisfies a no flux condition at those points. Also the south portion of the wall in the example coincides with a transport ray, but in that case the normal derivative of u is naturally negative in the boundary nodes, and it becomes zero only asymptotically in time (at the equilibrium). Then, by continuity arguments, the best choice seems to us simply to impose a no flux boundary condition for v at those points.
The direct application of scheme (33)-(37) is anyway not so efficient, due to the numerical difficulty of handling unbounded discontinuous solutions. In Fig. 6, the computed stationary solutions for (1) and their level lines are shown in the test example of Section 4 (compare with Fig. 5).
Despite the fact that the real sand flow is completely separated in the left and the right subregions of Ω, at the numerical level the flow travels through the grid points and  Figure 6. Numerical stationary solutions u and v of system (1) in the test example and their level lines then it can cross the separation line. More precisely, the transport path for sand from a point X in the right hand side should be the segment XP in Fig. 4; on the contrary, the algorithm splits this flow along vertical and horizontal segments connecting nodes and then part of this sand reaches the segment P Q even far from P (and from there eventually the left-hand side of the table). That is why the simple use at the discrete level of the same decomposition strategy adopted to characterize the stationary solutions is not able to reduce this phenomenon. As a test we applied in fact on the same uniform grid the scheme (33)-(37) separately in the two subregions of Ω, with suitable wall boundary conditions on the cut (the P Q segment). The results (see Fig. 7) show an evident improvement of the solutions only in the left (that is the regular) subregion.
Better results can be expected by coupling decomposition with suitable grid strategies. Keeping the uniform grid, the use of semi-lagrangian type schemes along characteristics should give a better trace of the correct transport directions. On the other hand, a different idea could be to employ unstructured grids (and mesh refinements near the singularity regions) in order to improve the accuracy. The discussion of these approaches will be the goal of a forthcoming paper. The main difficulty is, anyway, that a sharp domain decomposition requires the a priori knowledge of the ridge set, which is not in general an easy task. For example, if we slightly modify the table in Fig. 4 by simply opening a symmetric portion of the boundary on its northern side ({(x, y) : y = 1, 0.5 ≤ x ≤ 1}), the situation becomes completely different: a curved internal ridge appears, and with the help of the normal directions to the singular boundary points it subdivides the table into