AbstractWe investigate the origin of stellar metallicity gradients in massive galaxies at large radii (r > 2 R_eff) using ten cosmological zoom simulations of halos with 6 × 10¹²M_⊙ < M_halo < 2 × 10¹³M_⊙. 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 evolution of the mass-metallicity relation. At larger radii, the galaxies 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 (〈 ▽ Z_stars 〉= -0.35 dex/dex) gradients in agreement with observations. Metallicity gradients of models without winds are inconsistent with observations. For the wind model, stellar accretion is steepening existing in-situ metallicity gradients by about 0.2 dex by the present day and is required to match observed gradients. Metallicity gradients are significantly steeper for systems, which have accreted stars in minor mergers. In contrast, galaxies with major mergers have relatively flat gradients, confirming previous results. We highlight the importance of stellar accretion for stellar population properties of massive galaxies at large radii, which provide important constraints for formation models.