The hydraulic analysis of water distribution networks (WDNs) is divided into two approaches: namely, a demand-driven analysis (DDA) and a pressure-driven analysis (PDA). In the DDA, the basic assumption is that the nodal demand is fully supplied irrespective of the nodal pressure, which is mainly suitable for normal operating conditions. However, in abnormal conditions, such as pipe failures or unexpected increase in demand, the DDA approach may cause unrealistic results, such as negative pressure. To address the shortcomings of DDA, PDA has been considered in a number of studies. For PDA, however, the head-outflow relation (HOR) should be given, which is known to contain a high degree of uncertainty. Here, the DDA-based simulator, EPANET2 was modified to develop a PDA model simulating pressure deficient conditions and a Monte Carlo simulation (MCS) was performed to consider the quantitative uncertainty in HOR. The developed PDA model was applied to two networks (a well-known benchmark system and a real-life WDN) and the results showed that the proposed model is superior to other reported models when dealing with negative pressure under abnormal conditions. In addition, the MCS-based sensitivity analysis presents the ranges of pressure and available discharge, quantifying service reliability of water networks.