We propose and analyze finite volume Godunov type methods based on discontinuous flux for a system of non-linear partial differential equations proposed by Hadeler and Kuttler to model the dynamics of growing sandpiles generated by a vertical source on a flat bounded rectangular table. The scheme is made well-balanced by modifying the flux function locally by including source term as a part of the convection term. Its extension to multi-dimensions is not straightforward for which an approach has been introduced here based on Transport Rays. This approach is compared with another approach for inclusion of source term which uses the idea of inverting the divergence operator relying on the Curl-free component of the Helmholtz decomposition of the source term. Numerical experiments are presented to illustrate the efficiency of the proposed scheme for both unsteady and steady state calculations and to make comparisons with the previously studied finite difference and semi-Lagrangian approaches by Falcone and Finzi Vita.