  Motivated by the need for efficient and accurate simulation of the dynamics
of the polar ice sheets, we design high-order finite element discretizations
and scalable solvers for the solution of nonlinear incompressible Stokes
equations. We focus on power-law, shear thinning rheologies used in modeling
ice dynamics and other geophysical flows. We use nonconforming hexahedral
meshes and the conforming inf-sup stable finite element velocity-pressure
pairings $\mathbb{Q}_k\times \mathbb{Q}^\text{disc}_{k-2}$ or $\mathbb{Q}_k
\times \mathbb{P}^\text{disc}_{k-1}$. To solve the nonlinear equations, we
propose a Newton-Krylov method with a block upper triangular preconditioner for
the linearized Stokes systems. The diagonal blocks of this preconditioner are
sparse approximations of the (1,1)-block and of its Schur complement. The
(1,1)-block is approximated using linear finite elements based on the nodes of
the high-order discretization, and the application of its inverse is
approximated using algebraic multigrid with an incomplete factorization
smoother. This preconditioner is designed to be efficient on anisotropic
meshes, which are necessary to match the high aspect ratio domains typical for
ice sheets. We develop and make available extensions to two libraries---a
hybrid meshing scheme for the p4est parallel AMR library, and a modified
smoothed aggregation scheme for PETSc---to improve their support for solving
PDEs in high aspect ratio domains. In a numerical study, we find that our
solver yields fast convergence that is independent of the element aspect ratio,
the occurrence of nonconforming interfaces, and of mesh refinement, and that
depends only weakly on the polynomial finite element order. We simulate the ice
flow in a realistic description of the Antarctic ice sheet derived from field
data, and study the parallel scalability of our solver for problems with up to
383M unknowns.
