A model based on the elastic theory of continuously distributed dislocations that accounts for the transport of geometrically necessary dislocations (GND) on slip systems is developed. It allows keeping the crystallographic nature of glide by allocating velocities specific to slip systems to GND. At grain boundaries, the dislocation transport equation is resolved between a specific system in a grain and a specific system in the adjacent grain. It is used to simulate a compression creep test followed by unloading of a multiple slip deforming multi-crystal of ice during which kink band formation, grain boundary migration and localized grain nucleation are observed. The model predictions are compared to 2D strain fields obtained by digital image correlation and show a good agreement. Besides, the kink band position corresponds very well with an area of strong lattice misorientation predicted by the model and is also bounded by opposite densities of edge dislocations, in agreement with kink banding theory and characterization. Furthermore, the grain boundary migration is observed to happen from predicted low dislocation density area towards high dislocation ones, also in agreement with the theory. Lastly, the triple junctions where nucleation is observed are also characterized by high GND density and especially strong gradient of elastic energy density. These different features show the relevance of using a continuum theory of polarized dislocations per slip system to study the onset of relaxation mechanisms like kink banding, grain boundary migration and grain nucleation and possibly to propose nucleation and migration criteria.