Stability in the operation of water distribution systems (WDSs) is paramount to maintaining efficient and reliable water delivery. Nonlinear model predictive control (NMPC) emerged as a suitable control strategy due to WDSs’ inherent nonlinearity and cross-coupling dynamics. However, classical NMPC is formulated under a finite horizon and does not guarantee closed-loop stability. It also relies heavily on intricate model-based dynamics, a cumbersome and time-consuming process for large-scale WDSs. This paper proposes a comprehensive control strategy that employs a data-enabled model identification technique, replacing physics-based models and ensuring stability and recursive feasibility via quasi-infinite horizon NMPC. The main objective of this work is to satisfy the water demand at every time step while guaranteeing a stable pressure head and energy-efficient pump operation in the WDS. A complete stability and feasibility analysis of the control strategy is also provided. Extensive simulations validate the proposed method demonstrating (1) data-driven model accuracy with an unseen and noisy dataset exhibiting 0.01% error and (2) optimal WDS operation under nominal and robust conditions, ensuring demand compliance, cost-savings by 8% ($18k annually), and pressure head stability within 5% of the steady-state value.