# On a dual formulation for the growing sandpile problem

Serge Dumont, Noureddine Igbida

► To cite this version:

Serge Dumont, Noureddine Igbida. On a dual formulation for the growing sandpile problem. European Journal of Applied Mathematics, 2009, 20 (2), pp.169-185. 10.1017/S0956792508007754. hal-01113763

HAL Id: hal-01113763

https://hal.science/hal-01113763

Submitted on 13 Mar 2023

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
# On a dual formulation for the growing sandpile problem

S. Dumont and N. Igbida

LAMFA, Université de Picardie Jules Verne, CNRS UMR 6140, 33, rue Saint-Leu, 80 039 Amiens, France

email: [serge.dumont,noureddine.igbida]@u-picardie.fr

In this paper, we are interested in the mathematical and numerical study of the Prigozhin model for a growing sandpile. Based on implicit Euler discretization in time, we give a simple improvement of theoretical and numerical analyses of the dual formulation for the problem. By using this model, we also give some application to the Monge-Kantorovich problem for optimal mass transportation.

## 1 Introduction

We are interested in the mathematical and numerical study of a sandpile which builds up when sand is dropped on a table. Sand might fall off the sides of the table. If a sand grain falls to a position where the slope of the pile is steep, it will slide down, and eventually cause sand grains of the existing pile to slide down as well. Since the work of Prigozhin (cf. [20]), this has been well known (see also [2, 3, 16] and the references therein) that the evolution of the surface of the sandpile when the angle of stability is equal to $\pi/4$, can be described by the following PDE:

$$ egin{cases} u_t - \nabla \cdot (m \nabla u) = f \\ m \ge 0, \quad |\nabla u| \le 1, \quad m(|\nabla u| - 1) = 0 \\ u = 0 \\ u(0) = u_0 \end{cases} \quad \begin{cases} \text{in } Q := (0, T) \times \Omega, \\ \text{on } \Sigma := (0, T) \times \partial\Omega, \end{cases} \quad (P) $$

where $\Omega \subset \mathbb{R}^D$ ($D$ equal to 1 or 2 in practical situations) is a bounded open domain describing the table, $\partial\Omega$ is its boundary from which the sand may fall out. Solution $u$ is the height of the surface that grows up (resp. grows down) under sand addition (resp. sand removal) by a source called $f$. Actually, this equation is obtained by assuming that

(1) the flow of the granular material is confined in a thin boundary layer moving down the slopes of a growing pile and

(2) the density of the material is constant.

So, we can use conservation law

$$ \partial_t u + \nabla \cdot q = f \quad (1.1) $$
where $q$ is the horizontal projection of the material flux of the moving layer. Then, by neglecting the inertia and supposing that the surface flow is directed towards the steepest descent, we have

$$q = -m\nabla u, \qquad (1.2)$$

where $m = m(x, t) \ge 0$ is an unknown scalar function. So, equation (1.1) implies

$$\partial_t u - \nabla \cdot (m \nabla u) = f. \qquad (1.3)$$

The initial free surface is given by

$$u(x, 0) = u_0.$$

The surface slope cannot exceed repose angle $\alpha$ of the material, i.e.

$$|\nabla u(t, x)| \le \gamma = \tan(\alpha).$$

Moreover, there is no pouring over the parts of angle less than $\alpha$; so that

$$|\nabla u(x, t)| < \gamma \implies m(x, t) = 0. \qquad (1.4)$$

To close the problem, one needs to add boundary condition. In this paper, we consider the open table problem; i.e. Dirichlet boundary condition

$$u(t, x) = 0 \quad \text{for } (t, x) \in (0, T) \times \partial\Omega.$$

At last, assuming that $\alpha = \pi/4$, then $\gamma = 1$ and the evolution problem (P) follows.

It is clear that the stationary problem associated with (P) is

$$\left\{
\begin{array}{ll}
-\nabla \cdot (m \nabla u) = f & \text{in } \Omega, \\
m \ge 0, \quad |\nabla u| \le 1, \quad m(|\nabla u| - 1) = 0 & \text{on } \partial\Omega,
\end{array}
\right.
\qquad (S_1)$$

where $f$ is assumed to be not depending on $t$. It is well known by now, that $(S_1)$ appears in the study of the optimal mass transport problem of the Monge-Kantorovich type with the Euclidean distance as a cost function (cf. [1, 14, 15] and the references therein). Indeed, functions $u$ solving $(S_1)$ together with $m$ are meaningful in the context of the transport problem. More precisely, $(S_1)$ is a dual problem associated with the relaxed variational formulation by Kantorovich for optimal mass transportation of $f^+$ into $f^-$, where $f^+$ and $f^-$ are the positive and negative parts of $f$, respectively. Then, $m$ is the transport density, vector $-\nabla u$ is the direction of the optimal transportation and quantity $-m \nabla u$ gives the transport flux. For numerical simulation concerning $(S_1)$, i.e. numerical approximation of the transport flux that contains all the information about the optimal transportation, we use in this paper the evolution problem (P) and its large-time behaviour.

The theoretical study of (P) may be handled by the nonlinear semi-group theory of the evolution problem governed by a sub-differential operator, so that unknown function $m$
may be not considered. However, in contrast with ($S_1$), for the numerical analysis, function $m$ and its regularity are crucial both for problems ($P$) and ($S_1$). Roughly speaking, the theoretical study of ($P$) derived from the study of a minimizing problem with a gradient constraint. For the numerical study, one can use an associated dual formulation where the flux $\sigma = m \nabla u$ is minus flux, and the analysis of the problem depends on the space of dual variables (which would be in connection with the regularity of $\sigma = m \nabla u$).

Remember that, in general $m$ is a Radon measure and the gradient of $u$ needs to be taken in an unusual sense (for more details in this direction one can see papers [6–8] and the references therein). However, under additional assumptions on $f$ and $\Omega$, $m$ may be an $L^q$ function; for more details in this direction, one can see papers [10–12] and the references therein.

In [3], the authors use a dual problem in a space of vector valued Radon measures both to prove the existence and uniqueness of a solution for ($P$) and for the numerical analysis of the problem. Our approach is different and may be simpler. For the theoretical analysis and large-time behaviour we use nonlinear semi-group theory. As for the numerical analysis, we show that it is enough to study an associated dual formulation in $H_{\text{div}}(\Omega)$, the space of $L^2$ vector functions with $L^2$ divergence.

Let us give the main lines of our approach. First, it is not difficult to see that solution $u$ of ($S_1$) is also a solution of a minimizing problem; so that ($P$) is a particular case of the evolution problem governed by a sub-differential operator stated in $L^2(\Omega)$. Applying standard results of nonlinear semi-group theory in a Hilbert space (cf. [9]), we deduce that ($P$) has a unique variational solution (see Definition 2.2), and, as $t \to \infty$, this solution converges to a solution of stationary problem ($S_1$). For the numerical analysis, thanks to the nonlinear semi-group theory, we use the fact that the variational solutions may be obtained by the Euler implicit discretization in time. So, it is enough to study the stationary problem of type

$$ egin{cases} 
v - \nabla \cdot (m \nabla v) = g \\ 
m \ge 0, \quad |\nabla v| \le 1, \quad m(|\nabla v| - 1) = 0 \\ 
v = 0 \quad \text{on } \partial\Omega, 
\end{cases} 
\quad \text{in } \Omega 
\quad (S_2) $$

where $g \in L^2(\Omega)$. It is not difficult to see that solution $v$ of ($S_2$) corresponds to the projection, with respect to the $L^2$ norm, of $g$, onto the set of $1$-Lipschitz continuous functions vanishing on the boundary. For numerical analysis of the projections, we use duality arguments. This last argument was extensively used in previous papers (see [3, 20]) in the space of vector Radon measure valued. The main novelty of our paper is that we use it simply in $H_{\text{div}}(\Omega)$.

The paper is organized as follows: in the next section, we give a formal nonlinear semi-group approach of ($P$) to show the existence and uniqueness of a variational solution, the convergence of the approximate solution by Euler implicit discretization in time and the convergence, as $t \to \infty$, of the solutions of ($P$) to the solutions of ($S_1$). In section 3, we develop the duality argument for the projection, i.e. the solution of ($S_2$), and we show how one can obtain the solution of ($S_2$) by a dual variable in $H_{\text{div}}(\Omega)$. We prove the convergence of numerical scheme that we use for numerical approximation of the projection and for
the solution of (P). At last, in section 4, results of numerical simulations for (P) and (S₁) are given.

## 2 Semi-group approach: Existence, uniqueness, ε-approximate solution and asymptotic behaviour

For ε > 0, we say that $(t_i, f_i)_{i=1,...,n}$ is an ε-discretization for the problem, if $t_0 = 0 < t_1 < \dots < t_{n-1} < T = t_n$ with $t_i - t_{i-1} \le \varepsilon$, $f_1, \dots, f_n \in L^2(\Omega)$, such that

$$ \sum_{i=1}^{n} \int_{t_{i-1}}^{t_i} \|f(t) - f_i\|_{L^2(\Omega)} \le \varepsilon. $$

For any ε > 0, we say that $u_ε$ is an ε-approximate solution of (P), if there exists $(t_i, f_i)_{i=1,...,n}$ an ε-discretization for the problem such that

$$ u_{\varepsilon}(t) = \begin{cases} u_0 & \text{for } t \in ]0, t_1], \\ u_i & \text{for } t \in ]t_{i-1}, t_i], \quad i = 1, \dots, n \end{cases} \qquad (2.1) $$

and $u_i$ solves the Euler implicit time discretization of (P)

$$ \left\{ \begin{array}{l} |\nabla u_i| \le 1, \exists m_i \ge 0, m_i (|\nabla u_i| - 1) = 0 \quad \text{in } \Omega \\ u_i - \nabla \cdot (m_i \nabla u_i) = \varepsilon f_i + u_{i-1} \\ u_i = 0 \quad \text{on } \partial\Omega. \end{array} \right. \quad \text{in } \Omega \qquad (2.2) $$

It is clear that problem (2.2) is a particular case of stationary problem (S₂). To deal with this problem, let us consider

$$ K = \{z \in W^{1,\infty}(\Omega) \cap W_0^{1,2}(\Omega); |\nabla z| \le 1 \quad \text{a.e. in } \Omega\} $$

and $\Pi_K$ the indicator function of K defined by

$$ \Pi_K(z) = \begin{cases} 0 & \text{if } z \in K, \\ +\infty & \text{otherwise.} \end{cases} $$

We denote by $\mathbb{P}_K$ the projection onto convex K, with respect to the $L^2(\Omega)$ norm, i.e. $v = \mathbb{P}_K g$ if and only if

$$ \|v - g\|_{L^2(\Omega)} = \min_{z \in K} \|z - g\|_{L^2(\Omega)}. $$

The sub-differential of $\Pi_K$ in $L^2(\Omega)$ is given by $v \in \partial \Pi_K(g)$ if and only if

$$ \int_{\Omega} v(z - g) \le 0 \quad \text{for any } z \in K. $$

It is not difficult to see that, if $v \in K$ is a solution of (S₂) in the sense that there exists a measurable function $m$ such that $m \nabla v \in (L^1(\Omega))^D$, $m (|\nabla v| - 1) = 0$ a.e. in $\Omega$ and $v - \nabla \cdot (m \nabla v) = g$ in $\mathcal{D}'(\Omega)$, then $v = \Pi_K(g)$. This gives in particular the concept of variational solutions for problems (S₁), (S₂) and (P). More precisely, we get
**Definition 2.1** For a given $f \in L^2(\Omega)$ (resp. $g \in L^2(\Omega)$), we say that $v$ is a variational solution of $(S_1)$ (resp. $(S_2)$) if $v \in K$ and $\int_{\Omega} f (v - z) \ge 0$ (resp. $\int_{\Omega} (g - v) (v - z) \ge 0$) for any $z \in K$.

**Definition 2.2** For a given $f \in L^2_{\text{loc}}(0, T; L^2(\Omega))$ and $u_0 \in K$, we say that $u$ (resp. $u_\epsilon$) is a variational solution (resp. $\epsilon$-approximate variational solution) of (P) if $u \in W^{1,1}(0, T; L^2(\Omega)) \cap L^2(0, T; W_0^{1,\infty}(\Omega))$, $u(0) = u_0$ and, for any $t \in (0, T)$, $u(t) \in K$ and $\int_{\Omega} (f(t) - u_t(t)) (u(t) - z) \ge 0$ for any $z \in K$ (resp. $u_\epsilon$ is given by (2.1) and $u_i$ is a variational solution of (2.2)).

Thanks to [9], we know that $\partial \mathbb{I}_K$ is a maximal monotone graph in $L^2(\Omega)$, i.e. for any $g \in L^2(\Omega)$, there exists a unique $v$ solution of

$$v + \partial \mathbb{I}_K(v) \ni g \quad (2.3)$$

and, if $v_i$ is the solution corresponding to $g_i$ for $i = 1, 2$, then

$$\|v_1 - v_2\|_2 \le \|g_1 - g_2\|_2.$$

Moreover, if $v_i$ is the solution corresponding to $g_i$ for $i = 1, 2$, then

$$\|(v_1 - v_2)^{+2}\|_2 \le \|(g_1 - g_2)^{+2}\|_2.$$

So, by using the nonlinear semi-group theory in Hilbert space for evolution problems governed by a sub-differential operator (cf. [5, 9, 22]), we have the following result.

**Theorem 2.3** Let $u_0 \in K$, $T > 0$ and $f \in L^2_{\text{loc}}(0, T; L^2(\Omega))$. Then,

(1) For any $\epsilon > 0$ and any $\epsilon$-discretization of (P), there exists a unique $\epsilon$-approximate variational solution of (P).

(2) There exists a unique $u \in \mathcal{C}([0, T]; L^2(\Omega))$ such that $u(0) = u_0$, and, as $\epsilon \to 0$,

$$u_{\epsilon} \to u \quad \text{in} \quad \mathcal{C}([0, T]; L^2(\Omega)).$$

(3) The $u$ function given by (2) is the unique variational solution of (P).

Moreover, if for $i = 1, 2$ $u_i$ is the solution corresponding to $f_i$, then

$$\frac{d}{dt} \int_{\Omega} (u_1 - u_2)^+ \le \int_{\Omega} (f_1 - f_2)^+ \quad \text{in} \quad \mathcal{D}'(0, T)$$

In particular, if $f \ge 0$, then $u \ge 0$ a.e. in $\Omega$.

**Proof** The first part of the theorem is a simple consequence of the fact that variational solution $u_i$ of (2.2) is equal to $\mathbb{I}_K(u_{i-1} + \epsilon f_i)$. Since $\partial \mathbb{I}_K$ is a maximal monotone graph in $L^2(\Omega)$, the second part of the theorem is a consequence of the classical nonlinear semi-group theory (cf. [9], see also [5, 22]). The third part of the theorem can be proved
using the regularity results for semi-group solutions of evolution equations governed by sub-differential operators (cf. Theorem 3.6 of [9]). □

Remark Thanks to [5], we can derive precise estimates on error $||u - u_{\varepsilon}||_{L^2(\Omega)}$. Actually, by using Theorem 4.2 and Theorem 13.3 of [5], we have

$$ \sup_{t \in [0, T)} ||u - u_{\varepsilon}||_{L^2(\Omega)} \le C (T \varepsilon^{\frac{1}{2}} + \varepsilon). $$

where $C$ is a constant depending on $||u_0||_{L^2(\Omega)}$ and $||f||_{L^2(\Omega)}$ $u_{\varepsilon}$ is an $\varepsilon$-approximate variational solution of (P).

Another main interest of evolution equation (P) is its large-time behaviour. Actually, by using Theorem 3.11 of [9], we have

**Theorem 2.4** Let $f \in L_{loc}^2(0, \infty; L^2(\Omega))$, $u_0 \in K$ and $u$ be the variational solution of (P). If there exists $f_{\infty} \in L^2(\Omega)$ such that $f - f_{\infty} \in L^1(0, \infty; L^2(\Omega))$, then there exists $u_{\infty} \in K$ such that $u_{\infty}$ is a variational solution of (S$_1$), and, as $t \to \infty$, $u(t) \to u_{\infty}$ in $L^2(\Omega)$.

**Proof** Since, for any $C > 0$, the set $\{z \in L^2(\Omega); \mathbb{I}_K(z) + ||z||_{L^2(\Omega)} \le C\} = \{z \in K; ||z||_{L^2(\Omega)} \le C\}$ is a compact set in $L^2(\Omega)$, then the theorem follows using Theorem 3.11 of [9]. □

The main application we have in mind for this theorem is the optimal mass transport problem of the Monge-Kantorovich type. However, one sees that this theorem gives only the approximations of the potential $u$, i.e. the direction of the optimal transport. We do not know about the transport densities; in this direction one can see recent work [4].

### 3 Dual formulation and numerical approximation of the projection

Now, in order to give a numerical approximation of the solution of (P), we will use part 2. of Theorem 2.3. So, we focus our attention on the projection $\mathbb{P}_{Kg}$, for a given $g \in L^2(\Omega)$. Remember that $v = \mathbb{P}_K(g)$ if and only if $v \in K$ and

$$ J(v) = \frac{1}{2} ||v - g||_{L^2(\Omega)}^2 = \min_{z \in K} J(z). \quad (3.1) $$

To treat this minimization problem we use a dual formulation. By using standard duality argument (cf. [13]), the dual problem associated with (3.1) is given by the following functional

$$ G(w) = \frac{1}{2} \int_{\Omega} (\mathrm{div}(w))^2 + \int_{\Omega} g \ \mathrm{div}(w) + \int_{\Omega} |w|. $$

Indeed, (3.1) is equivalent to

$$ \min \{F(z) + H(\Lambda z); z \in \mathcal{C}_0^1(\Omega)\}, \quad (3.2) $$
where $\Lambda z := \nabla z$ is linear operator from $\mathcal{C}_0^1(\Omega)$ to $\mathcal{C}(\Omega)^N$ and $F : \mathcal{C}_0^1(\Omega) \to \mathbb{R}^+$ and $H : \mathcal{C}(\Omega)^N \to \overline{\mathbb{R}}$ are convex functions defined by

$$F(z) = \frac{1}{2} \int_{\Omega} |z - g|^2 \quad \text{and} \quad H(\sigma) = \begin{cases} 0 & \text{if } |\sigma(x)| \le 1 \ \forall x \in \Omega, \\ +\infty & \text{otherwise.} \end{cases}$$

Thanks to [13], the dual problem associated with (3.2) is given by

$$\sup\{-F^*(\Lambda^*\sigma) - H^*(-\sigma); \sigma \in (\mathcal{C}(\Omega)^N)^*\} \quad (3.3)$$

which corresponds to

$$\sup\{-G(\sigma); \sigma \in (\mathcal{C}(\Omega)^N)^*\}. \quad (3.4)$$

So, one sees that the natural space for the study of (3.2) by duality argument is the set of vector valued Radon measures such that the divergence is a $L^2$. This kind of argument was already used in previous papers (see [3, 20]). Our aim is to simplify the analysis by using simply the space $H_{\text{div}}(\Omega)$;

$$H_{\text{div}}(\Omega) = \{w \in (L^2(\Omega))^D; \text{div}(w) \in L^2(\Omega)\}.$$

First, we have

**Lemma 1** For any $g \in L^2(\Omega)$, $w \in H_{\text{div}}(\Omega)$ and $z \in K$, we have

$$-G(w) \le J(z).$$

**Proof** Let $w \in H_{\text{div}}(\Omega)$ and $z \in K$ be fixed. Writing $\frac{1}{2}(\text{div}(w) - (z - g))^2 \ge 0$, it is clear that

$$-\frac{1}{2} \int_{\Omega} (\text{div}(w))^2 - \int_{\Omega} g \, \text{div}(w) + \int_{\Omega} z \, \text{div}(w) \le \frac{1}{2} \int_{\Omega} (z - g)^2$$

and, since $|w| - \nabla z \cdot w \ge 0$, then

$$-\int_{\Omega} |w| - \int_{\Omega} z \, \text{div}(w) = -\int_{\Omega} |w| + \int_{\Omega} \nabla z \cdot w \le 0.$$

Adding the two inequalities, the result follows. $\square$

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

In general, it is not clear whether the extremality relation is fulfilled for $(w, v) \in H_{\text{div}}(\Omega) \times K$ or not. In other words, it is not clear if there exists $(w, v) \in H_{\text{div}}(\Omega) \times K$ such that $-G(w) = J(v)$. It is known that this kind of relation is very important in order to connect the primal and dual problems. Standard analysis leads to extremality relation in $\mathcal{M}_b(\Omega)^N \times K$ and one needs to define the gradient in a non-standard way (tangential gradient with respect to Radon measure). Here, we prove that the study of $G$ in $H_{\text{div}}(\Omega)$ gives some kind of approximation for extremality relation. This approximation gives a connection between primal and dual problems and enables us to compute $v$ the projection on $K$ of $g$ by computing $\sup_{\sigma \in H_{\text{div}}(\Omega)} -G(\sigma)$. This is the aim of the following theorems (Theorem 3.1 and Theorem 3.2).
**Theorem 3.1** Let $g \in L^2(\Omega)$ and $v = P_K(g)$. Then, there exists a sequence $(w_\varepsilon)_{\varepsilon>0}$ in $H_{div}(\Omega)$, such that, as $\varepsilon \to 0$,

$$ \int_{\Omega} |w_{\varepsilon}| \to \int_{\Omega} v(g-v), \qquad (3.5) $$

$$ \operatorname{div}(w_{\varepsilon}) \to v - g \quad \text{in } L^2(\Omega) \qquad (3.6) $$

and

$$ \begin{aligned} \lim_{\varepsilon \to 0} G(w_\varepsilon) &= \inf_{w \in H_{div}(\Omega)} G(w) \\ &= -\min_{z \in K} J(z) \qquad (3.7) \\ &= -\frac{1}{2} \int_{\Omega} |g - v|^2. \end{aligned} $$

To prove this result, let us consider the following elliptic equation

$$ \left\{ \begin{array}{ll} v_{\varepsilon} - \nabla \cdot w_{\varepsilon} = g & \text{in } \Omega, \\ w_{\varepsilon} = \phi_{\varepsilon}(\nabla v_{\varepsilon}) & \text{on } \partial\Omega. \\ v_{\varepsilon} = 0 & \end{array} \right. \qquad (S_{\varepsilon}) $$

where, for any $\varepsilon > 0$, $\phi_\varepsilon : \mathbb{R}^D \to \mathbb{R}^D$ is given by

$$ \phi_{\varepsilon}(r) = \frac{1}{\varepsilon}(|r| - 1)^{+} \frac{r}{|r|}, \quad \text{for any } r \in \mathbb{R}^{D}. $$

It is not difficult to see that $\phi_\varepsilon$ satisfies the following properties:

(i) for any $r_1, r_2 \in \mathbb{R}^D$, $(\phi_\varepsilon(r_1) - \phi_\varepsilon(r_2)) \cdot (r_1 - r_2) \ge 0$

(ii) there exists $\varepsilon_0 > 0$ and $A > 1$ such that $\phi_\varepsilon(r) \cdot r \ge |r|^2$ for any $|r| \ge A$ and $\varepsilon < \varepsilon_0$

(iii) for any $\varepsilon > 0$ and $r \in \mathbb{R}$, $|\phi_\varepsilon(r)| \le \phi_\varepsilon(r) \cdot r$

So, for any $g \in L^2(\Omega)$, $(S_\varepsilon)$ has a unique solution $v_\varepsilon$, in the sense that $v_\varepsilon \in H_0^1(\Omega)$, $w_\varepsilon := \phi_\varepsilon(\nabla v_\varepsilon) \in L^2(\Omega)^D$ and $v_\varepsilon - \nabla \cdot w_\varepsilon = g$ in $\mathcal{D}'(\Omega)$. We are interested in the study of the limit of $(v_\varepsilon, w_\varepsilon)$, as $\varepsilon \to 0$.

**Lemma 2** We have

(1) $(v_\varepsilon)_{0<\varepsilon<\varepsilon_0}$ is bounded in $H_0^1(\Omega)$.

(2) For any Borel set $B \subset \Omega$,

$$ \liminf_{\varepsilon \to 0} \int_B |\nabla v_\varepsilon| \le |B|. $$

**Proof** Taking $v_\varepsilon$ as a test function in $(S_\varepsilon)$, we get

$$ \int_{\Omega} v_{\varepsilon}^{2} + \frac{1}{\varepsilon} \int_{\Omega} (|\nabla v_{\varepsilon}| - 1)^{+} |\nabla v_{\varepsilon}| = \int_{\Omega} f v_{\varepsilon}. \qquad (3.8) $$

So,

$$ \|v_{\varepsilon}\|_{L^2} \le \|f\|_{L^2(\Omega)} \qquad (3.9) $$
and

$$ \frac{1}{\epsilon} \int_{\Omega} (|\nabla v_{\epsilon}| - 1)^{+} |\nabla v_{\epsilon}| \le \|f\|_{L^2(\Omega)}^2. \quad (3.10) $$

Using property (ii) of $\phi_{\epsilon}$, for any $0 < \epsilon < \epsilon_0$, we have

$$ \begin{aligned} \int_{\Omega} |\nabla v_{\epsilon}|^2 &= \int_{[|\nabla v_{\epsilon}| \le A]} |\nabla v_{\epsilon}|^2 + \int_{[|\nabla v_{\epsilon}| > A]} |\nabla v_{\epsilon}|^2 \\ &\le \int_{[|\nabla v_{\epsilon}| \le A]} |\nabla v_{\epsilon}|^2 + \frac{1}{\epsilon} \int_{(|\nabla v_{\epsilon}| - 1)^{+} |\nabla v_{\epsilon}|} \\ &\le |A|^2 |\Omega| + \|f\|_{L^2(\Omega)}^2. \end{aligned} $$

Thus $v_{\epsilon}$ is bounded in $H_0^1(\Omega)$. Now, let $B \subseteq \Omega$ be a fixed Borel set. We have,

$$ \begin{aligned} \int_B |\nabla v_\epsilon| &\le \int_B (|\nabla v_\epsilon| - 1)^+ + |B| \\ &\le \int_B (|\nabla v_\epsilon| - 1)^+ |\nabla v_\epsilon| + |B| \\ &\le \epsilon \int_\Omega (f - v_\epsilon) v_\epsilon + |B|. \end{aligned} $$

Letting $\epsilon \to 0$, and using the fact that $v_{\epsilon}$ is bounded in $L^2(\Omega)$, the second part of the lemma follows. $\square$

**Proof of the theorem** As a consequence of Lemma 2, there exists $\tilde{v} \in H_0^1(\Omega)$ and a subsequence that we denote again by $\epsilon$, such that

$$ v_{\epsilon} \to \tilde{v} \quad \text{in } H_0^1(\Omega) - \text{weak and in } L^2(\Omega), $$

which implies that

$$ \operatorname{div}(w_{\epsilon}) \to \tilde{v} - g \quad \text{in } L^2(\Omega). $$

To prove that $\tilde{v} \in K$, let us consider $A_{\delta} = [|\nabla \tilde{v}| \ge 1 + \delta]$, with arbitrary $\delta > 0$. Since, as $\epsilon \to 0$, $\nabla v_{\epsilon} \to \nabla \tilde{v}$ in $L^1(\Omega)^p$-weak, then

$$ \begin{aligned} |A_{\delta}| &\le \frac{1}{1+\delta} \int_{A_{\delta}} |\nabla \tilde{v}| \\ &\le \frac{1}{1+\delta} \liminf_{\epsilon \to 0} \int_{A_{\delta}} |\nabla v_{\epsilon}|, \end{aligned} $$

so that, by using the second part of Lemma 2, we deduce that

$$ |A_{\delta}| \le \frac{1}{1+\delta} |A_{\delta}|, $$
which implies that $|A_\delta| = 0$. Since $\delta > 0$ is arbitrary, then, we deduce that $|\nabla \tilde{v}| \le 1$ a.e. in $\Omega$. Now, let us prove that $\tilde{v} = \mathbb{P}_K(g)$. For any $z \in K$, we have

$$
\begin{align*}
\int_{\Omega} (g - \tilde{v}) (\tilde{v} - z) &= \lim_{\epsilon \to 0} \int_{\Omega} (g - v_{\epsilon}) (\tilde{v} - z) \\
&= \lim_{\epsilon \to 0} \int_{\Omega} \phi_{\epsilon}(\nabla v_{\epsilon}) \cdot \nabla (\tilde{v} - z) \\
&= \lim_{\epsilon \to 0} \int_{\Omega} (\phi_{\epsilon}(\nabla v_{\epsilon}) - \phi_{\epsilon}(\nabla z)) \cdot \nabla (\tilde{v} - z) \\
&\ge 0.
\end{align*}
$$

To prove (3.5), we see, first, that by using property (iii) of $\phi_\epsilon$, we have

$$
\begin{align*}
\limsup_{\epsilon \to 0} \int_{\Omega} |w_{\epsilon}| &\le \limsup_{\epsilon \to 0} \int_{\Omega} \phi_{\epsilon}(\nabla v_{\epsilon}) \cdot \nabla v_{\epsilon} \\
&\le \limsup_{\epsilon \to 0} \int_{\Omega} (g - v_{\epsilon}) v_{\epsilon} \\
&\le \int_{\Omega} (g - v) v.
\end{align*}
$$

Moreover, we have

$$
\begin{align*}
\int_{\Omega} v (g - v) &= \lim_{\epsilon \to 0} \int_{\Omega} v (g - v_{\epsilon}) \\
&= \lim_{\epsilon \to 0} \int_{\Omega} w_{\epsilon} \cdot \nabla v \\
&\le \liminf_{\epsilon \to 0} \int_{\Omega} |w_{\epsilon}|,
\end{align*}
$$

which ends up the proof of (3.5). At last, since by Lemma 1 we know that $-G(w_\epsilon) \le \min_{z \in K} J(z)$, then by using (3.6) and (3.7), we deduce that, as $\epsilon \to 0$, $G(w_\epsilon) \to -\frac{1}{2} \int_\Omega |v-g|^2 = -\min_{z \in K} J(z)$. This completes the proof of the theorem. $\square$

*Remark* Theorem 3.1 also gives an extremality relation. Indeed, if $J(v) = \min_{z \in K} J(z)$; i.e. $v = \mathbb{P}_K(g)$, then, thanks to Theorem 3.1, by taking $w \in (M_b(\Omega))^D$ the weak limit of $w_\epsilon$, in $(M_b(\Omega))^D$, we deduce that $v$ is given by

$$
\left\{
\begin{array}{ll}
v - \operatorname{div}(w) = g & \text{in } \mathcal{D}'(\Omega), \\
|w|(\Omega) = \int_{\Omega} v(g-v). &
\end{array}
\right.
\qquad (3.11)
$$

This gives a new characterization of the projection of $g$ onto $K$ (see [3, 18, 19] for more details in this direction).

We are ending this section by giving results for the numerical approximation of dual problem

$$
\inf_{w \in H_{\text{div}}(\Omega)} G(w).
$$
Since we want to use a finite element method, we suppose in the following that:

* Domain $\Omega$ is a bounded, open, polyhedral subset of $\mathbb{R}^D$ ($D$ equal to 1 or 2).

* $T_h$ will be a regular partitioning (triangulation or quadrangulation) of $\bar{\Omega}$ by $n$ disjoint open simplices $\tau$ of diameter no greater than a given real $h$, with $\bar{\Omega} = \cup_{\tau \in T_h} \bar{T}_\tau$.

Let $V_h \subset V := H_{\text{div}}(\Omega)$ be the space of lowest-order Raviart-Thomas finite elements (cf. [21]):

$$ V_h = \{q_h \in (L^2(\Omega))^D : q_h|_{\tau} = a_{\tau} + b_{\tau}x, a \in \mathbb{R}^2, b \in \mathbb{R}, \forall \tau \in T_h, \text{ and } q_h \cdot v \text{ is continuous across simplex boundaries}\}. $$

Space $V_h$ is a finite dimensional subspace of $V$ with a dimension equal to $N = N(h)$. Let us denote by $r_h$ the interpolation operator onto $V_h$ given in Theorem 6.1 of [21]. By using the bounds for error interpolation (cf. Theorem 6.3 of [21]) and the density of $\mathcal{D}(\Omega)$ in $H_{\text{div}}(\Omega)$, it is not difficult to prove that, for any $w \in H_{\text{div}}(\Omega)$, as $h \to 0$, we have

$$ \begin{array}{l} r_h(w) \to w \quad \text{in } (L^2(\Omega))^N \\ \text{and} \\ \text{div}(r_h(w)) \to \text{div}(w) \quad \text{in } L^2(\Omega). \end{array} \qquad (3.12) $$

Besides, it is clear that convex optimizing problem

$$ \inf\{G(q_h); q_h \in V_h\} $$

has a solution. Let us denote a solution by $w_h$, i.e.

$$ G(w_h) = \inf\{G(q_h); q_h \in V_h\}. \qquad (3.13) $$

**Theorem 3.2** Let $g \in L^2(\Omega)$, $v = P_k g$ and $w_h$ a solution of (3.13). Then, as $h \to 0$,

$$ \text{div}(w_h) \to v - g \quad \text{in } L^2(\Omega) \qquad (3.14) $$

and

$$ -G(w_h) \to \min_{z \in L^2(\Omega)} J(z) = \frac{1}{2} \|v - g\|_{L^2(\Omega)}. \qquad (3.15) $$

**Proof** Thanks to (3.13) and the fact that $r_h(w_\epsilon) \in V_h$, we have $G(w_h) \le G(r_h(w_\epsilon))$, where $w_\epsilon$ is given by Theorem 3.1. In other words,

$$ \frac{1}{2} \int_{\Omega} (\text{div}(w_h))^2 + \int_{\Omega} g \text{ div}(w_h) + \int_{\Omega} |w_h| \le \frac{1}{2} \int_{\Omega} (\text{div}(r_h(w_\epsilon)))^2 + \int_{\Omega} g \text{ div}(r_h(w_\epsilon)) + \int_{\Omega} |r_h(w_\epsilon)|, $$
which implies that

$$ \frac{1}{2} \int_{\Omega} (\mathrm{div}(w_h) - v + g)^2 = \frac{1}{2} \int_{\Omega} (\mathrm{div}(w_h))^2 + \frac{1}{2} \int_{\Omega} (v-g)^2 - \int_{\Omega} \mathrm{div}(w_h) (v-g) \quad (3.16) $$

$$ \le \frac{1}{2} \int_{\Omega} (\mathrm{div}(r_h(w_{\varepsilon})))^2 + \int_{\Omega} g \ \mathrm{div}(r_h(w_{\varepsilon})) + \int_{\Omega} |r_h(w_{\varepsilon})| \\ - \int_{\Omega} |w_h| - \int_{\Omega} v \ \mathrm{div}(w_h) + \frac{1}{2} \int_{\Omega} (v-g)^2 \quad (3.17) $$

Since $v \in K$, then by integrating by part we get

$$ - \int_{\Omega} v \ \mathrm{div}(w_h) - \int_{\Omega} |w_h| = \int_{\Omega} \nabla v \cdot w_h - \int_{\Omega} |w_h| \\ \le 0; $$

so that (3.18) implies that

$$ \frac{1}{2} \int_{\Omega} (\mathrm{div}(w_h) - v + g)^2 \le \frac{1}{2} \int_{\Omega} (\mathrm{div}(r_h(w_{\varepsilon})))^2 \\ + \int_{\Omega} g \ \mathrm{div}(r_h(w_{\varepsilon})) + \int_{\Omega} |r_h(w_{\varepsilon})| + \frac{1}{2} \int_{\Omega} (v-g)^2 . $$

Letting $h \to 0$ and using (3.12), we obtain

$$ \limsup_{h \to 0} \frac{1}{2} \int_{\Omega} (\mathrm{div}(w_h) - v + g)^2 \le \frac{1}{2} \int_{\Omega} (\mathrm{div}(w_{\varepsilon}))^2 + \int_{\Omega} g \ \mathrm{div}(w_{\varepsilon}) + \int_{\Omega} |w_{\varepsilon}| \\ + \frac{1}{2} \int_{\Omega} (v-g)^2. $$

At last, thanks to Theorem 3.1, by letting $\varepsilon \to 0$, we get

$$ \limsup_{h \to 0} \frac{1}{2} \int_{\Omega} (\mathrm{div}(w_h) - v + g)^2 \le 0 $$

and the result follows.

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

## 4 Numerical simulations

At this stage, it is necessary to introduce an approximation of term $\int_{\tau} |w_h|$ dx for each simplex $\tau$ of the partitioning of $T_h$. We have chosen here to write

$$ \int_{\tau} |w_h| \mathrm{d}x \approx |\tau||w_h|(P_{\tau}), $$

where $|\tau|$ represents the area of simplex $\tau$, and $P_{\tau}$ is one of the vertices of $\tau$.
Using this approximation, the solution of (3.13) at each time $i\Delta t$ ($i \in \mathbb{H}$) is a minimizer of non-differential functional:

$$ 
\begin{aligned}
G_h: \quad & \mathbb{R}^n \to \mathbb{R} \\
w_h & \mapsto G_h(w_h) = \frac{1}{2} \|\mathrm{div}(w_h)\|_{L^2(\Omega)} + (t_h^i, \mathrm{div}(w_h)) + \gamma \sum_{\tau \in T_h} |\tau||w_h(P_\tau)| \\
& := \frac{1}{2}(Aw_h, w_h) + (t_h^i, \mathrm{div}w_h) + \gamma \sum_{\tau \in T_h} |\tau||w_h(P_\tau)|
\end{aligned} 
\quad (4.1) $$

where $n$ is the dimension of $V_h$, $A$ is an $n \times n$ positive semi-definite matrix, $t_h^i$ is a vector of $\mathbb{R}^n$ defined by $t_h^i = \Delta t f_h^i + u_{i-1}$ (see formula (2.2) with $\varepsilon = \Delta t$), and $\gamma = \tan(\alpha)$ is the sand angle of repose (chosen equal to 1 in the theoretical study). In the following, we denote $w_h$ a minimizer of this functional.

The minimization of this functional is performed by a Gauss Seidel type algorithm that can be summarized as follows (see for instance [17]):

* Initiate the algorithm with a vector $q_0 \in \mathbb{R}^n$, and, for $k \ge 0$ until convergence, chose a canonical direction $e_j$ in $\mathbb{R}^n$ and find $\rho_{jk}$ minimizing

$$ 
\begin{aligned}
\varphi_{jk}: \quad & \mathbb{R} \to \mathbb{R} \\
\rho & \mapsto G_h(q_k + \rho e_j).
\end{aligned} 
$$

* Take $q_{k+1} = q_k + \rho_{jk}\omega e_j$, where $\omega > 0$ is an over-relaxation parameter.

* When $\varphi_{jk}$ is differentiable, a Newton algorithm is used to find $\rho_{jk}$. Otherwise, $\rho_{jk}$ can be computed directly (because is this case, $\varphi_{jk}$ is the sum of a polynomial of degree two and an absolute value).

* This algorithm is performed until $\|q_{k+1} - q_k\|_{\ell^2(\mathbb{R}^n)} \le \varepsilon$, for a given convergence criterion $\varepsilon$. Afterwards, take $w_h = q_k$.

**Remark** The advantage of this technique over a regularisation procedure coupled with the resolution of a nonlinear PDE is that, at each step $k$, the minimization procedure can be made only on indices $j$ where sand may fall, i.e. where the function $\varphi_{jk}$ is not independent of $\rho$.

In all the simulation below, we have chosen relaxation parameter $\omega = 1.2$, and convergence criterion $\varepsilon = 10^{-5}$.

Domain $\Omega$ is square $\Omega = (-1, 1)^2$. For the discretization of the problem, we use the Raviart-Thomas elements of the lowest order [21] on a regular square grid. In this case, for a given integer $N \in \mathbb{H}$, the step of discretization $h$ is equal to $\frac{2}{N+1}$ and the simplices are squares of the form $(j_1h, (j_1+1)h) \times (j_2h, (j_2+1)h)$ for $0 \le j_1, j_2 \le N$.

Then, knowing a minimizer $w_h$ of (3.13), solution $u_i$ of (2.2) is computed using extremality relation (3.11) in a weak sense with piecewise finite elements $P0$.

In the first numerical simulation of a growing sandpile, the density of the source of sand $f(x, t)$ is a constant equal to 1, for all time and for all $x \in \Omega$ where $\Omega = (-1, 1)^2$. In this example, the exact free surface $u(x, t)$ of the sandpile can be easily computed. The numerical simulation has been done with a regular quadrangular $60 \times 60$ grid for the
<figure><img src="image_4.png" /><figcaption>FIGURE 1. Sandpile surface at time $t = 0.45$ and $t = 1.415$ for $f = 1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 2. Sand flux at time $t = 0.45$ and $t = 1.415$ for $f = 1$.</figcaption></figure>

discretization of $\Omega$, with a discretization in time $\Delta t = 0.01$, and with the maximal angle of stability $\alpha$ such that $\tan \alpha = \sqrt{2}$. We found an accuracy of primal variable $u$ equal to $\{\max_i |u(\cdot, t_i) - u_h(\cdot, t_i)|_{0,\infty,\Omega} / |u|_{0,\infty,\Omega} \times (0, T) \le 0.033$. Then, this relative error is closed to the space step of discretization $h = 2/N$. Figure 1 shows the sandpile at time $t = 0.45$, and at time $t = 1.415$ when solution $u$ becomes stationary and equal to the euclidean distance to the boundary.

Figure 2 shows the flux on the sandpile surface. We can see that the flux is parallel to the gradient of the surface and vanishes on the diagonal of square $\Omega$.

In the second example, the initial free surface of the sandpiles is equal to 1, and the source of sand is constant and is a negative measure around point $x_0 = (0, 0.4)$. The Neumann boundary conditions $q \cdot v = m \frac{\partial u}{\partial v} = 0$ are considered, which simulated the existence of a wall on the boundary. Even if the theoretical background of this problem is not studied in the article, the results shown in Figure 3 show that the method is accurate. The error on the conservative volume law $| \int_0^T \int_{\Omega} f(x,t) dx dt - \int_{\Omega} u(T,x) dx |$ is always less than $10^{-14}$ for all $T \le 4$.

We can also see in Figure 5 that the flux is concentrated around $x_0$ and has a $\frac{1}{|x-x_0|}$ singularity. Moreover, the flux vanishes on the boundary even at $t = 3$ when the height of the sandpile does not vanish on the boundary.

In the third example, the density of the source of sand $f$ is independent of $t$ and changes signs. The source is strictly positive on the disc centered on $x_1 = (-0.5, 0)$ with a radius equal to 0.2 and on the disc centered on $x_2 = (0, -0.6)$ with a radius equal to $\frac{1}{\sqrt{15}}$, and it is strictly negative on the disc centered on $x_3 = (0, 0)$ with a radius equal to 0.2 and vanishes elsewhere. On each disc, the distribution of sand is quadratic, with a maximum
<figure><img src="image_6.png" /><figcaption>FIGURE 3. Sandpile surface and sand flux at time $t = 0.14$, $t = 3.00$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 4. Sandpile surface at time $t = 0.24$, $t = 0.43$, $t = 0.50$ and $t = 3$.</figcaption></figure>

on the center of the disc. In other words, the sand is taken from $x_3$ and put over points $x_1$ and $x_2$, such that $\int_{\Omega} f(x, t) dx = 0$ for any time $t$.

Figure 4 shows the sandpile at different times, and Figure 5 shows the sand flux for the same times. The support of the sand source is included into the circles plotted on
<figure><img src="image_8.png" /><figcaption>FIGURE 5. Sand flux at time $t = 0.24$, $t = 0.43$, $t = 0.50$ and $t = 3$.</figcaption></figure>

the figure. We can see that for any time $t \le t_1 = 0.24$, the two sandpiles and the hole do not interact. After this time, the sand falls from the sandpile centering around $x_1$ into the hole. At time $t = 0.5$, sand falls from the two sandpiles to the hole. After all, at time $t = 3$, we obtain a stationary state: all the sand deposited over point $x_1$ falls into the hole, and a part of the sand deposited over $x_2$ falls into the hole, the other part creeps through the boundary.

## Acknowledgements

This work was supported by the French A.N.R. Grant JC05-41831. We would like to thank Olivier Goubet and Véronique Martin for many interesting discussions on this subject. Thanks to L. Guyot for reviewing the English for the paper.

## References

[1] AMBROSIO, L. (2003) Lecture notes on optimal transport Problems, in Mathematical Aspects of Evolving Interfaces. In: *Lecture Notes in Mathematics*, LNM 1812, Springer, Berlin.

[2] ARONSON, G., EVANS, L. C. & WU, Y. (1996) Fast/Slow diffusion and growing sandpiles. *J. Differential Equations* **131**, 304–335.

[3] BARRETT, J. W. & PRIGOZHIN, L. (2006) Dual formulation in critical state problems. *Interface Free Bound.* **8**, 349–370.

[4] BARRETT, J. W. & PRIGOZHIN, L. (2007) A mixed formulation of the Monge–Kantorovich equations. *M2AN*, **41**(6), 1041–1060.

[5] BENILAN, PH., CRANDALL, M. G. & PAZY, A. *Evolution Equations Governed by Accretive Operators*, Book to appear.
[6] BOUCHITTÉ, G. & BUTTAZZO, G. (2001) Characterization of optimal shapes and masses through Monge-Kantorovich equation. *J. Eur. Math. Soc.* **3**(2), 139–168.

[7] BOUCHITTÉ, G., BUTTAZZO, G. & SEPPECHER, P. (1997) Energies with respect to a measure and applications to low dimensional structures. *Calc. Var. Partial Differential Equations* **5**, 37–54.

[8] BOUCHITTÉ, G., BUTTAZZO, G. & SEPPECHER, P. (1997) Shape optimization solutions via Monge-Kantorovich. *C.R. Acad. Sci. Paris, t.324, Série I*, 1185–1191.

[9] BRÉZIS, H. (1973) Opérateurs maximaux monotones et semigroups de contractions dans les espaces de Hilbert. (French). *North-Holland Mathematics Studies, No. 5. Notas de Matemàtica (50).* North-Holland Publishing Co., Amsterdam-London; American Elsevier Publishing Co., Inc., New York, 1973.

[10] DE PASCALE, L. & PRATELLI, A. (2002) Regularity properties for Monge transport density and for solutions of some shape optimization problem. *Calc. Var. Partial Differential Equations* **14**(3), 249–274.

[11] DE PASCALE, L. & PRATELLI, A. (2004) Sharp summability for Monge transport density via interpolation. *ESAIM Contrôle Optim. Calc. Var.* **10**, 549–552.

[12] DE PASCALE, L. EVANS, L. C. & PRATELLI, A. (2004) Integral estimates for transport densities. *Bull. London Math. Soc.*, **36**(3), 383–395.

[13] EKELAND, I. & TÉMAM, R. (1999) Convex analysis and variational problems. In *Classics in Applied Mathematics*, 28. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA.

[14] EVANS, L. C. (1997) Partial differential equations and Monge-Kantorovich mass transfer. In *Current Developments in Mathematics*, International Press, Boston, MA, pp. 65–126.

[15] EVANS, L. C. & GANGBO, W. (1999) Differential equations methods for the Monge-Kantorovich mass transfer problem. *Mem. Am. Math. Soc.*, **137**, 653.

[16] EVANS, L. C. & REZAKHANLOU, F. (1998) A stochastic model for sandpiles and its continuum limit. *Comm. Math. Phys.* **197**(2), 325–345.

[17] GLOWINSKI, R., LIONS & TRÉMOLIÈRES (1976) Analyse numérique des inéquations variationnelles. *Méthodes Mathématiques de l'Informatique*, 5. Dunod, Paris.

[18] IGBIDA, N. On Equivalent Formulations for Monge-Kantorovich Equation. submitted.

[19] IGBIDA, N. Evolution Monge-Kantorovich equation. preprint.

[20] PRIGOZHIN, L. Variational model of sandpile growth. *Euro. J. Appl. Math.* **7**, 225–236.

[21] ROBERTS, J. E. & THOMAS, J.-M. (1991) Mixed and hybrid methods. In: P. G. Ciarlet & J. L. Lions (editors), *Handbook of Numerical Analysis, Finite Element Methods (Part 1)*, Vol. II, North-Holland, Amsterdam.

[22] SHOWALTER, R. E. (1997) Monotone operators in Banach space and nonlinear partial differential equations. In *Mathematical Surveys and Monographs*, **49**. American Mathematical Society, Providence, RI.