We investigate the evolution of stellar population gradients from $z=2$ to $z=0$ in massive galaxies at large radii ($r > 2R_{mathrm{eff}}$) using ten cosmological zoom simulations of halos with $6 times 10^{12} M_{odot} < M_{mathrm{halo}} < 2 times 10^{13}M_{odot}$. The simulations follow metal cooling and enrichment from SNII, SNIa and AGB winds. We explore the differential impact of an empirical model for galactic winds that reproduces the mass-metallicity relation and its evolution with redshift. At larger radii the galaxies, for both models, become more dominated by stars accreted from satellite galaxies in major and minor mergers. In the wind model, fewer stars are accreted, but they are significantly more metal poor resulting in steep global metallicity ($langle abla Z_{mathrm{stars}} rangle= -0.35$ dex/dex) and color (e.g. $langle abla g-r rangle = -0.13$ dex/dex) gradients in agreement with observations. In contrast, colour and metallicity gradients of the models without winds are inconsistent with observations. Age gradients are in general mildly positive at $z=0$ ($langle abla Age_{mathrm{stars}} rangle= 0.04$ dex/dex) with significant differences between the models at higher redshift. We demonstrate that for the wind model, stellar accretion is steepening existing in-situ metallicity gradients by about 0.2 dex by the present day and helps to match observed gradients of massive early-type galaxies at large radii. Colour and metallicity gradients are significantly steeper for systems which have accreted stars in minor mergers, while galaxies with major mergers have relatively flat gradients, confirming previous results. This study highlights the importance of stellar accretion for stellar population properties of massive galaxies at large radii, which can provide important constraints for formation models.