We present a numerical model of trace-element oscillatory zoning patterns formed when zircon crystallizes from silicate melt, which is also appropriate for other accessory phases with known partition and diffusion coefficients and saturation conditions. The model accounts for diffusion-controlled accessory mineral growth and the equilibrium crystallization of major mineral phases. Consideration of recent, experimentally determined dependencies of partition coefficients on temperature, we find that thermal changes provide the simplest explanations of oscillatory zoning in accessory minerals. Numerical experiments with different cooling rates explore different crystallization scenarios with and without the precipitation of other phases and/or the interface reaction of phosphorus (P) and yttrium (Y) to form xenotime. However, these processes are monotonically related to growth rate and do not cause oscillations. Minor 3–10 °C variations in temperature do not result in zircon dissolution, but strongly influence zircon growth and lead to variations in coeval Y, Hf, and rare earth element (REE) concentrations of up to a factor of two, comparable to those observed in nature. Such temperature variations may be very common in any igneous body in response to external factors such as replenishment by hotter magmas or convection. More significant temperature fluctuations may result in initial minor dissolution at higher temperatures during a mafic recharge event but with continuous growth afterward. At high temperature (>∼850 °C) the amplitude of oscillations is relatively small that confirms observations of both less common oscillatory zoning in hot and dry volcanic rhyolites and abundant oscillations in plutonic zircons and in zircons in cold and wet crystal-rich mushes. Additional oscillations in zircon are modeled in response to oscillations of pressure on the order of ±35–50 bars, causing water concentration fluctuations of ±0.1 wt% in water-saturated melt cells with a gas bubble. These variations cause variations of Zr diffusion and zircon growth rates. Such fluctuations could result from pressure oscillations due to recharge and convection in the magma chamber. All simulated runs generate smoothed oscillatory growth zoning; similar patterns found in nature may not necessarily require post-growth intracrystalline diffusion.