Techniques for Solving the Stokes Equations of Ice Sheet Dynamics at Continental Scale

Ice sheets exhibit incompressible creeping flow with shear-thinning rheology. On a continental scale, the flow is characterized by localized regions of fast flow that are separated from vast slow regions by thin transition zones. To address these issues, we wish to use a parallel, adaptive mesh, higher-order finite element discretization to model in 3D the equations for ice sheet dynamics on a continental scale. Our approach to solving the resulting nonlinear equations is a fast Newton method with a Krylov subspace iteration to solve the linearized system of equations. We enforce incompressibility with element-discontinuous pressure spaces with uniformly bounded inf-sup constants for discrete stability. The Krylov iterations are preconditioned using block preconditioners, with algebraic multigrid (AMG) for the viscous stress block and an approximate commutator approach to preconditioning the Schur complement of the system. The matrices for higher-order discretizations are relatively dense and poorly preconditioned by AMG, so sparser low order discretizations are used to construct the AMG preconditioner. An important difference between ice sheet flows and similar processes is the friction at the base of the ice sheet, which can be significant. When this friction is modeled by a Robin-type boundary condition, terms are added to the viscous stress block of the Stokes operator. These changes have implications for preconditioning both this block and the Schur complement of the system, which we will address. The length scales of a continental ice sheet necessitate highly anisotropic discretizations. Anisotropy drastically impairs the efficiency of standard multilevel preconditioners. We will discuss techniques for overcoming this when preconditioning the viscous stress block, such as semi-coarsening and robust smoothers, and how they can be implemented using popular solver libraries.