## 1 Introduction  
The modelling of a flowing granular mass over three-dimensional terrain is complex considering the interface between the mass and base terrain, as well as the continually moving free surface. The governing equations of the continuum dynamics models for avalanches, i.e., partial differential equations of mass and momentum conservation, can be efficiently solved using various numerical simulation methods, e.g., finite-element methods. However, difficulties such as instabilities have to be overcome in numerical simulations.  
The numerical stability of the scheme is mainly determined by the advective acceleration terms that dominate the motion of the avalanche. Suitable upwinding needs to be adopted in these instances.  
Several traditional numerical schemes have been adopted to solve hyperbolic-type partial differential equations, e.g., the first-order upwind scheme and second-order upwind beam-warming and Lax-Wendroff schemes. A Eulerian scheme similar to the MacCormack scheme [1] was adopted by [PERSON] and [PERSON] [2], and Lagrangian finite-difference schemes were adopted by [PERSON] et al. [3], [PERSON] [4], and [PERSON] [5]. These methods may cause dispersivity, which leads to large oscillations during calculations. Therefore, to avoid such oscillations, diffusive terms have to be artificially added to the right sides of the governing equations in these methods, which may affect the resolution quality.  
Shock wave formation in particulate flows is a common phenomenon in granular mass flows over an inclined plane with obstacles [6; 7; 8; 9; 10]. The traditional finite-difference schemes mentioned above fail to capture such phenomena. A high-resolution shock-capturing numerical scheme is needed to simulate the moving front accurately. Since [PERSON] and [PERSON] [11] proposed the second-order-accurate non-oscillatory central total-variation-diminishing (TVD) scheme, research on central TVD schemes has become intensive, with modifications made in successive studies [12; 13; 14], and the scheme has been recognized as an effective and commonly used approach. The central TVD schemes have already been applied to simulate the movement of dense and shallow granular \"liquid\". [PERSON] et al. [15] adopted the TVD Lax-Friedrichs scheme to simulate some typical phenomena when granular flow strikes an obstacle, such as shock waves, vacuums, and dead areas. Following [PERSON] et al. [15], a non-oscillatory central difference scheme integrated with the TVD limited or weighted essentially non-oscillatory cell formation approach was introduced by [PERSON] et al. [16] to model shock waves and obtain a smooth solution. To model the movement of landslides on natural terrain, [PERSON] et al. [17] developed an adaptive mesh and parallel Godunov solver. [PERSON] and [PERSON] [18] applied the Riemann solver algorithm to model granule-laden liquid flows. However, the above-mentioned numerical methods restrain the oscillatory behavior in the solution at the expense of increased algorithm complexity. In contrast, the TVD-MacCormak explicit scheme is a simple and efficient finite-difference method that has been applied to simulate the motion of Newtonian fluid.  
To solve traditional depth-averaged continuum dynamics models for avalanches, the present paper chooses the two-step (predictor and corrector) MacCormack scheme as the basic framework. A TVD term is introduced to eliminate steep-gradient oscillations. The validity of the numerical method is tested in an experiment on uniform granular flume flow and granular flow in a dam breaking scenario.