Geochemical models are usually fueled by experimental data to improve their ability to predict the environmental fate of trace elements. This was typically the case for the prediction of the binding properties towards trace metals of dissolved organic matter (Dom) in soils. The empirical parameterization of the binding properties of humic susbtances was used as a surrogate to predict trace metal speciation in soil solution by considering a ubiquitous homogeneity of Dom binding properties. A decade of modelling of copper (Cu) speciation in the solution of agricultural soils using the Windermere humic aqueous model (Wham) led however to challenge the homogeneity of Dom binding properties and to suggest that rigorously supervised model optimization can be used conversely to evidence their heterogeneity. We illustrated this consideration in a recent experiment mimicking under controlled conditions the field application of 31 organic residues on an agricultural soil (plus a non-amended control). Soil treatments induced a one-unit pH change, a 4-fold Dom concentration change, and a 10-fold change in total, available, and free Cu concentrations in the soil. When considering homogeneous binding properties of Dom, the measured speciation of Cu in soil solutions was very poorly predicted (R2 < 0.01; RMSR = 1.0). In contrast, the optimization of Dom binding properties for each single datum point within physically meaningful values (i.e. the range of binding site density and Cu affinity constant was restricted according to the literature review) and by accounting for the input uncertainty fitted neatly Cu2+ activity in soil solutions (R2 = 0.92, RMSR = 0.1). The optimized Dom binding properties consequently differed by 10-fold between the soil treatments with a high level of statistical significance. Previously achieved datasets showed however that Wham-optimized Dom binding properties were only occasionally correlated to Dom properties that can be routinely characterized by e.