Geophysical hazards such as avalanches and landslides necessitate the development of efficient predictive models for early warning systems. Depth-averaged models employing $\mu(I)$-rheology have become a standard approach due to their computational efficiency compared to complex 3D Navier-Stokes solvers. However, several significant scientific challenges persist, including the numerical preservation of "lake-at-rest" steady states, the incorporation of non-hydrostatic effects, and the proper implementation of local coordinates. Regarding the well-balanced property, granular stationary states differ fundamentally from standard fluids; they are characterized by an inequality where friction forces balance gravity, resulting in an infinite family of equilibria with non-flat free surfaces. Current state-of-the-art methods often rely on first-order approximations or regularization techniques that introduce artificial residual velocities. To address these limitations, we propose a novel reconstruction procedure enabling the definition of high-order well-balanced schemes. Numerical validation demonstrates that the proposed scheme achieves second and third-order accuracy in both space and time, while remaining exactly well-balanced for a broad class of stationary solutions. Furthermore, non-hydrostatic effects may be incorporated and a new geometric formulation may be introduced . This approach provides a more general and robust shallow model for dry avalanches, ensuring high-fidelity simulations across complex topographies and providing a reliable tool for hazard assessment.