Summary: This paper develops two first-order system – in this context, first-order refers to the order of the PDE not to the model – least-squares, fluidity-based formulations of a nonlinear Stokes flow model for ice sheets that attempt to overcome the difficulties introduced by unbounded viscosity. One commonly used way to define viscosity, Glen’s law, allows viscosity to become unbounded as the strain rates approach zero. Often, numerical approaches overcome these singularities by modifying viscosity to limit its maximum. The formulations in this paper, however, reframe the problem to avoid viscosity altogether by defining the system in terms of the inverse of viscosity, which is known as fluidity. This results in a quantity that approaches zero as viscosity approaches infinity. Additionally, a set of equations that represent the curl of the velocity gradient is added to help approximate the solution in a space closer to \(H^1\), which improves algebraic multigrid convergence. Previous research revealed that the first-order system least-squares formulation has difficulties in maintaining optimal discretization convergence on more complex domains. This paper discovers that this problem is linked to how the curl equations are scaled and that stronger scalings result in better solver performance but worse discretization convergence. Determining if there is an optimal scaling that balances performance and convergence is still an open question. Additionally, the fluidity-based formulations are tested with three 2D benchmark problems. Two of these benchmark problems involve basal sliding and one involves a time-dependent free surface. The fluidity-based solutions are consistent with the standard Galerkin method using Taylor-hood elements while better resolving viscosity.