Existence of a solution to the quasi-variational inequality problem arising in a modelfor sand surface evolution has been an open problem for a long time. Another long-standingopen problem concerns determining the dual variable, the flux of sand pouring down theevolving sand surface, which is also of practical interest in a variety of applications ofthis model. Previously, these problems were solved for the special case in which theinequality is simply variational. Here, we introduce a regularized mixed formulationinvolving both the primal (sand surface) and dual (sand flux) variables. We derive,analyse and compare two methods for the approximation, and numerical solution, of thismixed problem. We prove subsequence convergence of both approximations, as the meshdiscretization parameters tend to zero; and hence prove existence of a solution to thismixed model and the associated regularized quasi-variational inequality problem. One ofthese numerical approximations, in which the flux is approximated by thedivergence-conforming lowest order Raviart–Thomas element, leads to an efficient algorithmto compute not only the evolving pile surface, but also the flux of pouring sand. Resultsof our numerical experiments confirm the validity of the regularization employed.