This work measures the evolution of the iron content in galaxy clusters by a rigorous analysis of the data of 130 clusters at 0.1 < z < 1.3. This task is made difficult by a) the low signal-to-noise ratio of abundance measurements and the upper limits; b) possible selection effects; c) boundaries in the parameter space; d) non-Gaussian errors; e) the intrinsic variety of the objects studied; and f) abundance systematics. We introduce a Bayesian model to address all these issues at the same time, thus allowing cross-talk (covariance). On simulated data, the Bayesian fit recovers the input enrichment history, unlike in standard analysis. After accounting for a possible dependence on X-ray temperature, for metal abundance systematics, and for the intrinsic variety of studied objects, we found that the present-day metal content is not reached either at high or at low redshifts, but gradually over time: iron abundance increases by a factor 1.5 in the 7 Gyr sampled by the data. Therefore, feedback in metal abundance does not end at high redshift. Evolution is established with a moderate amount of evidence, 19 to 1 odds against faster or slower metal enrichment histories. We quantify, for the first time, the intrinsic spread in metal abundance, 18 ± 3%, after correcting for the effect of evolution, X-ray temperature, and metal abundance systematics. Finally, we also present an analytic approximation of the X-ray temperature and metal abundance likelihood functions, which are useful for other regression fitting involving these parameters. The data for the 130 clusters and code used for the stochastic computation are provided with the paper.