Current analyses of groundwater flow and transport typically rely on a physically-based model (PBM), which is inherently subject to error and uncertainty from multiple sources including model structural error, parameter error and data error. The model uncertainty can be difficult to quantify, and is propagated to the prediction. In this study, complementary data-driven models (DDMs) are used to improve prediction of groundwater flow models. The DDMs, trained with the historical residual of the PBM, have the capability to compensate for the defects of PBM. Five machine learning techniques, instance-based weighting (IBW), locally weighted regression (\textit {loess}), decision trees (DT), artificial neural networks (ANN) and support vector regression (SVR) are employed to construct the DDMs, and their performance of enhancing the prediction of the PBM is compared. Before the DDMs updating, cluster analysis is implemented on the dataset to improve the robustness and efficiency of the framework. The framework is tested in two real-world case studies based on the Republic River Compact Association (RRCA) model and the Spokane Valley Rathdrum Prairie (SVRP) model. The DDMs reduce the root-mean-square errors (RMSE) of the temporal, spatial and temporal plus spatial head prediction of the RRCA model by 82\%, 60\% and 48\% respectively. In the SVRP case study, the DDMs reduces the temporal head forecast of the PBM by 79\%. Localized DDMs that are conditioned on each cluster outperform global DDMs without clustering. It is also demonstrated that clustering significantly reduces the computational cost of training and cross validation of the DDMs. After clustering, the run-time of DDMs is negligible comparing with the PBM, which makes the framework very computationally efficient.