We use a centimeter-scale 2-D numerical model to investigate the effect of the presence of a second phase with various volume percent, shape, and orientation on strain localization in a viscoelastic matrix. In addition, the evolution of bulk rheological behavior of aggregates during uniaxial compression is analyzed. The rheological effect of dynamic recrystallization processes in the matrix is reproduced by viscous strain softening. We show that the presence of hard particles strengthens the aggregate, but also causes strain localization and the formation of ductile shear zones in the matrix. The presence of soft particles weakens the aggregate, while strain localizes within the particles and matrix between particles. The shape and the orientation of second phases control the orientation, geometry, and connectivity of ductile shear zones. We propose an analytical scaling method that translates the bulk stress measurements of our 2-D simulations to 3-D experiments. Comparing our model to the laboratory uniaxial compression experiments on ice cylinders with hard second phases allows the analysis of transient and steady-state strain distribution in ice matrix, and strain partitioning between ice and second phases through empirical calibration of viscous softening parameters. We find that the ice matrix in two-phase aggregates accommodates more strain than the applied bulk strain, while at faster strain rates some of the load is transferred into hard particles. Our study illustrates that dynamic recrystallization processes in the matrix are markedly influenced by the presence of a second phase.