In this chapter, the flow of grounded ice sheets on geophysical scales is investigated. Two ice flow configurations are considered: plane and radially-symmetric. Assuming that ice viscosities depend on local temperature, strain-rate and current strength of anisotropy of the material, computational models have been developed to solve the system of equations governing the flow of a large, gravity-driven, polythermal polar ice sheet. The plane flow problem is solved by a finite-element method, whereas the radially-symmetric problem is solved by applying a method of asymptotic expansions in a small parameter defining the ratio of an ice sheet’s thickness to its lateral span. The results of numerical simulations illustrate the effect of ice anisotropy on both the free-surface profile and the velocity field in a polar ice sheet. In addition, the influence of the bed topography features on the overall flow of an ice sheet is examined. The chapter is complemented with the presentation of results showing the effects of the dynamic recrystallization process on the flow of a polar ice sheet.