To establish a theoretical basis for predicting and interpreting the behavior of rapid mass movements on Earth's surface, we develop and test a new computational model for gravity-driven motion of granular avalanches across irregular, three-dimensional (3-D) terrain. The principles embodied in the model are simple and few: continuum mass and momentum conservation and intergranular stress generation governed by Coulomb friction. However, significant challenges result from the necessity of satisfying these principles when deforming avalanches interact with steep and highly variable 3-D terrain. We address these challenges in four ways. (1) We formulate depth-averaged governing equations that are referenced to a rectangular Cartesian coordinate system (with z vertical) and that account explicitly for the effect of nonzero vertical accelerations on depth-averaged mass and momentum fluxes and stress states. (2) We compute fluxes of mass and momentum across vertical cell boundaries using a high-resolution finite volume method and Roe-type Riemann solver. Our algorithm incorporates flux difference splitting, an entropy correction for the flux, and eigenvector decomposition to embed the effects of driving and resisting forces in Riemann solutions. (3) We use a finite element method and avalanche displacements predicted by Riemann solutions to compute Coulomb stresses conjugate to the displacements in 3-D stress space. (4) We test the model output against analytical solutions, a sand cone conceptual experiment, and (in a companion paper) data from detailed laboratory experiments. Model results illustrate a complex interplay of basal traction and internal stress, and they successfully predict not only the gross behavior but also many details of avalanche motion from initiation to deposition.