We analyze and model the dynamics of a complex hydro-geochemical system associated with a variety of physico-chemical processes. These take place in a soil-water environment mimicked through a laboratory-scale column hosting a redox zonation and subject to injection of a solution rich in the pharmaceutical diclofenac. Experimental evidences suggest that sorption is a key driver to diclofenac fate, while dissolution-precipitation and biotic redox reactions control the evolution of speciation patterns in pore water. Resting on these evidences, we develop an original reactive transport model to interpret dynamics of diclofenac in the considered scenario. Our model formulation is characterized by a high parameterization degree. We frame the analysis in a stochastic context and perform model calibration against available data through a maximum likelihood approach driven by global sensitivity analysis. We further develop a set of reduced-complexity models to explore the potential of diverse formulations to $$(a)$$ understand dominant processes (and parameters) driving system dynamics and $$(b)$$ effectively assist stochastic model calibration and ensuing uncertainty quantification. We assess the relative skill of each model to data interpretation through the Kashyap discrimination criterion. We then quantify predictive uncertainty emanating from model inputs to outputs. Our results reveal that the original system model is generally outperformed by its reduced-complexity counterparts. Our findings also demonstrate that relying on a modeling framework yielding a robust uncertainty quantification associated with estimates of geochemical heterogeneity patterns is key to provide a sound description of the fate of emerging contaminants such as diclofenac as they migrate through the subsurface.