Pharmaceutical products (PPs) are emerging water pollutants with adverse environmental and health-related impacts, owing to their toxic, persistent, and undetectable microscopic nature. Globally, increasing scientific knowledge and advanced technologies have allowed researchers to study PP-associated problems and their removal for water reuse. Experimental modeling methods require laborious, lengthy, expensive, and environmentally hazardous lab-work to optimize the process. On the other hand, predictive machine learning (ML) models can trace the complex input-output relationship of a process using available datasets. In this study, ensemble ML techniques, including decision tree (DT), random forest (RF), and Xtreme gradient boost (XGB), were used to explore PP (diclofenac, iopromide, propranolol, and trimethoprim) removal by a managed aquifer recharge (MAR) system. The model input parameters included characteristics of reclaimed water and soil used in the columns, pH, dissolved organic carbon, operating time, nitrogen dioxide, sulfate, nitrate, electrical conductivity, manganese, and iron. The selected PP removal was the model output. Datasets were collected through a one-year experimental study of continuous MAR system operation to predict the removal of PPs. DT, RF, and XGB models were then developed for one of the selected compounds and tested for the others to check the reliability of the ML model results. The developed models were assessed using statistical performance matrices. The experimental results showed >80% removal of propranolol and trimethoprim; however, removal of diclofenac and iopromide was only ≈50% by the MAR system. The proposed DT and RF models presented higher coefficients of determination (R2 ≥ 0.92) for diclofenac, propranolol, and trimethoprim than for iopromide (R2 ≤ 0.63). In contrast, the XGB model showed better results for diclofenac, iopromide, propranolol, and trimethoprim, with R2 values of 0.92, 0.72, 0.96, and 0.97, respectively. Therefore, XGB could be the best predictive model to provide insight into the adaptation of ML models to predict PP removal by the MAR system, thereby minimizing experimental work.