This paper highlights the use of incorporating strain gradient into flow stress to study localization behavior in materials. Pioneered by Zbib and Aifantis in the late 1980s, the formulation enabled incorporation of length scales into continuum formulations naturally. The formulation has also evolved into being able to study the effects of microstructure and heterogeneity on localization in granular materials. A multi-slip Mohr–Coulomb type plasticity model with the flow stress in the constitutive equation modified with a higher order gradient term of the effective plastic strain is used for this purpose. The possibility of abrupt changes of mobilized friction caused by intense shearing rate often leads to particle breakage. Its effects on localization are accounted for by modifying the material properties such as mobilized friction using a scaling parameter averaged over a representative elementary area. The change of shearing rate in the integration points was monitored through quasi-statistically measure parameter called inertia number. The inertia number was set to be less than l.0 × 10−3 all the time to consider quasi-static. The formulation was implemented into a finite element code and used to simulate plane strain compression tests on dry sand. The model highlights the effects of confining pressure, anisotropic microstructure, and the non-coaxial angle between the direction of principal stress and principal plastic strain rate directions on shear band characteristics.