The main goal of this article is to analyze a three-dimensional model for stress and velocity fields in grounded glaciers and ice sheets including the role of normal deviatoric stress gradients. This model leads to a nonlinear system of stationary partial differential equations for the velocity with a viscosity depending on the stress-tensor but which is not explicitly depending on the velocity. The existence and uniqueness of a weak solution corresponding to this model is established by using the calculus of variations. The approximation of this model is made by a finite element method with piecewise polynomial functions of degree 1 on a tetrahedral mesh and error analysis is performed. Numerical solutions show that the theoretical results we have obtained are almost optimal.