We propose a mesh-free method to solve the full Stokes equation for modeling the glacier dynamics with nonlinear rheology. Inspired by the Deep-Ritz method proposed in [13], we first formulate the solution to the non-Newtonian Stokes equation as the minimizer of a variational problem with boundary constraints. Then, we approximate its solution space by a deep neural network. The loss function for training the neural network is a relaxed version of the variational form, in which penalty terms are used to present soft constraints due to mixed boundary conditions. Instead of introducing mesh grids or basis functions to evaluate the loss function, our method only requires uniform sampling from the physical domain and boundaries. Furthermore, we introduce a re-normalization technique in the neural network to address the significant variation in the scaling of real-world problems. Finally, we illustrate the performance of our method by several numerical experiments, including a 2D model with the analytical solution, the Arolla glacier model with realistic scaling and a 3D model with periodic boundary conditions. Numerical results show that our proposed method is efficient in solving the non-Newtonian mechanics arising from glacier modeling with nonlinear rheology.