A phase-field model is proposed to study the inhomogeneous growth of zirconia induced by the stress–oxidation interaction, which captures the complex interplay among diffusion, oxidation kinetics, interfacial morphology evolution, and stress variation in an oxidation process. Through this numerical model, many experimentally observed but insufficiently understood phenomena can be well explained. Specifically, the numerical simulations reveal quantitatively the causes of interface roughening or smoothening during the inward oxide growth, the roughness-dependent oxide growth rate, and the nucleation sites of premature cracking. These numerical findings can be used as the theoretical references for the improving the durability of oxide scale and prolonging the service life of zirconium-based alloy cladding used in the nuclear power plant.