Magnesium is distinguished by its highly anisotropic inelastic deformation involving a profuse activity of deformation twinning. Instrumented micro/nano-indentation technique has been widely applied to characterize the mechanical properties of magnesium, typically through the analysis of the indentation load–depth response, surface topography, and less commonly, the post-mortem microstructure within the bulk material. However, experimental limitations prevent the real-time observation of the evolving microstructure. To bridge this gap, we employ a recently-developed finite-strain model that couples the phase-field method and conventional crystal plasticity to simulate the evolution of the indentation-induced twin microstructure and its interaction with plastic slip in a magnesium single-crystal. Particular emphasis is placed on two aspects: orientation-dependent inelastic deformation and indentation size effects. Several outcomes of our 2D computational study are consistent with prior experimental observations. Chief among them is the intricate morphology of twin microstructure obtained at large spatial scales, which, to our knowledge, represents a level of detail that has not been captured in previous modeling studies. To further elucidate on size effects, we extend the model by incorporating gradient-enhanced crystal plasticity, and re-examine the notion of ‘smaller is stronger’. The corresponding results underscore the dominant influence of gradient plasticity over the interfacial energy of twin boundaries in governing the size-dependent mechanical response.