In ice sheet and glacier modelling, the Finite Element Method is rapidly gaining popularity. However, constructing and updating meshes for ice sheets and glaciers is a non-trivial and computationally demanding task due to their thin, irregular, and time dependent geometry. In this paper we introduce a novel approach to ice dynamics computations based on the unfitted Finite Element Method CutFEM, which lets the domain boundary cut through elements. By employing CutFEM, complex meshing and remeshing is avoided as the glacier can be immersed in a simple background mesh without loss of accuracy. The ice is modelled as a non-Newtonian, shear-thinning fluid obeying the p-Stokes (full Stokes) equations with the ice atmosphere interface as a moving free surface. A Navier slip boundary condition applies at the glacier base allowing both bedrock and subglacial lakes to be represented. Within the CutFEM framework we develop a strategy for handling non-linear viscosities and thin domains and show how glacier deformation can be modelled using a level set function. In numerical experiments we show that the expected order of accuracy is achieved and that the method is robust with respect to penalty parameters. As an application we compute the velocity field of the Swiss mountain glacier Haut Glacier d'Arolla in 2D with and without an underlying subglacial lake, and simulate the glacier deformation from year 1930 to 1932, with and without surface accumulation and basal melt.