An operator-splitting algorithm is presented for the solution of a partial differential equation arising in the modeling of deposition processes in sand mechanics. Sand piles evolution is modeled by an advection-diffusion equation, with a non-smooth diffusion operator that contains a point-wise constraint on the gradient of the solution. Piecewise linear finite elements are used for the discretization in space. The advection operator is treated with a stabilized SUPG finite element method. An augmented Lagrangian method is proposed for the discretization of the fast/slow non-smooth diffusion operator. A penalization approach, together with a Newton method, is used for the treatment of inequality constraints. Numerical results are presented for the simulation of sand piles on flat and non-flat surfaces, and for extensions to water flows.