The groundwater (GW) resource plays a central role in securing water supply in the coastal region of Bangladesh and therefore the future sustainability of this valuable resource is crucial for the area. However, there is limited research on the driving factors and prediction of phosphate concentration in groundwater. In this work, geostatistical modeling, self-organizing maps (SOM) and data-driven algorithms were combined to determine the driving factors and predict GW phosphate content in coastal multi-aquifers in southern Bangladesh. The SOM analysis identified three distinct spatial patterns: K+Na+pH, Ca2+Mg2+NO₃−, and HCO₃−SO₄2−PO43−F−. Four data-driven algorithms, including CatBoost, Gradient Boosting Machine (GBM), Long Short-Term Memory (LSTM), and Support Vector Regression (SVR) were used to predict phosphate concentration in GW using 380 samples and 15 prediction parameters. Forecasting accuracy was evaluated using RMSE, R2, RAE, CC, and MAE. Phosphate dissolution and saltwater intrusion, along with phosphorus fertilizers, increase PO43− content in GW. Using input parameters selected by multicollinearity and SOM, the CatBoost model showed exceptional performance in both training (RMSE = 0.002, MAE = 0.001, R2 = 0.999, RAE = 0.057, CC = 1.00) and testing (RMSE = 0.001, MAE = 0.002, R2 = 0.989, RAE = 0.057, CC = 0.998). Na+, K+, and Mg2+ significantly influenced prediction accuracy. The uncertainty study revealed a low standard error for the CatBoost model, indicating robustness and consistency. Semi-variogram models confirmed that the most influential attributes showed weak dependence, suggesting that agricultural runoff increases the heterogeneity of PO43− distribution in GW. These findings are crucial for developing conservation and strategic plans for sustainable utilization of coastal GW resources.