Shallow flow or thin liquid film models are used for a wide range of physical and engineering problems. They are applicable to free surface flows, when the flow thickness is distinctly smaller than the lateral dimension. Respective situations exist in automotive and aircraft engineering (soiling, icing, exterior water management), industrial processes (coating, fire suppression), environmental flows (river and lake hydrodynamics, atmospheric flows), and natural hazards (floods, avalanches, debris flow). Shallow flow models allow capturing the free surface of the fluid with little effort and reducing the three-dimensional problem to a quasi two-dimensional problem through depth-integrating the flow fields. Despite remarkable progress of such models in the last decade, accurate description of complex topography remains a challenge. Interaction with topography is particularly critical for granular flows, because their rheology requires modeling of the pressure field, which is strongly linked to surface curvature and associated centrifugal forces. Shallow granular flow models are usually set up in surface-aligned curvilinear coordinates, and velocity is represented as a two-dimensional surfacealigned vector field. The transformation from Cartesian to curvilinear coordinates introduces fictitious forces, however, which result in complex governing equations. In this paper, we set up the shallow flow model in three-dimensional Cartesian coordinates and preserve three-dimensional velocity in the depthintegrated model. Topography is taken into account with a constraint on velocity. This approach is commonly applied by the thin liquid film community. The advantage is a curvature-free mathematical description that is convenient for complex topographies. The constraint on velocity yields a solution for the pressure field, which is required for the pressure-dependent rheology of granular materials. The model is therefore well-suited for granular flows on three-dimensional terrain, e.g., avalanches. The mathematical model is solved with a second-order accurate, implicit finite area scheme, based on the open source software OpenFOAM. We conduct numerical simulations for various cases to verify the numerical routine based on comparisons with analytical results and a grid refinement study. We establish connections to former solutions and discuss advantages and drawbacks of the proposed method from both mechanical and numerical perspectives.