The ability of softening strain gradient plasticity to simulate intragranular plastic slip localization modes called slip and kink bands, at is investigated at incipient plasticity. A strain gradient crystal plasticity model based on a quadratic energy contribution of the Nye tensor has been implemented within the massively parallel FFT-based solver AMITEX_FFTP. It is shown that a classical evaluation of the double curl operator does not properly account for the backstress induced by GNDs in slip and kink bands. Consequently, a better suited differentiation operator is proposed for this purpose. Besides, four practical implementations of interface conditions at grain boundaries are implemented and discussed. It is demonstrated that these conditions have a strong influence on GND induced hardening but weakly affect the formation of slip localization patterns.A theoretical and numerical study of an idealized kink band shows that this modeling framework, contrary to classical crystal plasticity models, is able to capture their essential physical features. A quantitative analysis of high resolution 3D polycrystalline simulations, based on imaged processing and a peak detection on plastic slip profiles, is applied to show that gradient effects strongly influence the mechanisms of selection between slip and kink banding, to the favor of the first. In particular, some kink bands are shown to decompose into an homogeneous distribution of parallel slip bands of finite length. Finally, the competition between structural effects and energetic cost of lattice curvature in the formation of localized slip bands networks in polycrystals is discussed. The present results demonstrate the relevance of softening crystal plasticity modeling including an energetic contribution of GND density to accurately capture slip localization in polycrystals.