Machine learning is increasingly used to predict pollutant removal in adsorption studies, but models are often fitted to small, replicated experimental datasets and evaluated with a single random train/test split. We examined how the validation design affects model selection and surrogate-based optimization for removing phenolic compounds (quantified as tannins) from sugarcane vinasse using citric-acid-modified sugarcane bagasse, based on 315 measurements (105 conditions in triplicate) spanning temperature, contact time, and vinasse dilution. Eight regressors, from linear and quadratic models to kernel methods, tree ensembles, Gaussian processes, and neural networks, were tuned and evaluated by repeated nested cross-validation that kept all replicates of a condition together, and their errors were compared with the replicate (pure-error) standard deviation of 3.56% removal. Under this replicate-aware validation, all non-linear models performed similarly (RMSE 3.89–4.10%, R2 0.83–0.85), whereas the linear model was clearly worse (RMSE 5.06%). A naive random split increased R2 by only about 0.01 (mean ΔR2 0.005–0.013), because, with replicate scatter of this size, memorizing replicates gives no advantage over a smooth model. Holding out entire dilution levels, however, separated the models sharply: the quadratic model remained the most accurate (RMSE 4.58%), whereas random forest and gradient boosting failed (R2 ≤ 0.13). Permutation and SHAP analyses identified dilution as the dominant variable. Used as the optimization surrogate, the quadratic model placed maximum removal at 10% v/v dilution in a flat near-optimal region spanning about 37–40 °C and 74–130 min (maximum 97.7%); a confirmatory experiment at 40 °C and 100 min, inside this region, gave 97.4 ± 0.4% (n = 5) against 97.5% predicted. For small replicated adsorption datasets, the held-out unit should match the intended use of the model, and predictions away from the sampled levels require models that extrapolate smoothly.