Production optimization of a water flooding project aims to determine the optimum production strategy for maintaining the reservoir under real-time optimal controls, and maximizing the oil production or net present value (NPV) of the project. Current production optimization method with a conventional grid-based simulator bears a high computational cost. Also, cumbersome history matching procedures are often required to build a reliable reservoir model for reducing the uncertainty of the optimized production strategies. To address the drawbacks of such conventional method, data-driven surrogate models, e.g., the correlation-based model, are developed to acquire interwell connectivity; such history-matched model can be then used for production optimization. It is noted that, however, the correlation-based models cannot handle the case where the conversion of duty occurs for some wells draining the reservoir, e.g., a production well is converted to a producer. Recently, we have developed a physics-based data-driven interwell numerical simulation model (INSIM) which can predict the production rates of individual phases and handle conversion-of-duty scenarios between any well pair. INSIM is derived by solving the material balance equation and augmenting the Buckley–Leverett theory. In this study, the INSIM is applied for the first time to optimize the production of water flooding projects. We derive and improve the INSIM by considering the active bottom or edge aquifers that are often encountered in real reservoirs. The improved INSIM can handle two types of aquifers: constant pressure aquifer and constant volume aquifer. By applying the inversion theory based on a finite-difference algorithm, we estimate the model parameters by history matching the historic production data. Subsequently, the history-matched INSIM is used for optimizing NPV of water flooding projects in both a synthetic reservoir and a real reservoir with an active aquifer. The injection rate of each injection well is treated as a tuning variable. As for the synthetic model, the simulation results show that the inversely derived model parameters are able to characterize the interwell formation properties, and the optimization result based on INSIM is in a good agreement with that obtained with the commercial simulator Eclipse 100. Similar history-matching results are also obtained when INSIM is applied to a real fracture-vuggy reservoir with an active aquifer. After production optimization of this real reservoir with INSIM, we achieve a 4% reduction in the water cut but a 20% increase in the oil production rate, which verifies the soundness of the proposed method.