## Abstract

Modeling of pressure-dependent users’ consumption is mandatory to simulate accurately the hydraulics of water distribution networks (WDNs). Several software solutions already exist for this purpose, but none of them actually permits the easy integration and testing of new physical processes. In this paper, we propose a new Python simulator that implements a state-of-the-art pressure-dependent model (PDM) of users’ consumptions based on the Wagner’s pressure–outflow relationship (POR). We tested our simulator on eight large and complex WDNs, for different levels of users’ demands. The results show similar precision and efficiency to the ones obtained by the authors of the original model with their MATLAB implementation. Moreover, in case of fully satisfied users’ demands, our simulator provides the same results as EPANET 2.0 in comparable computational times. Finally, our simulator is integrated into the open-source, collaborative, multi-platform, and Git versioned Python framework OOPNET (Object-Oriented Python framework for water distribution NETworks analyses); thus, it can be easily reused and/or extended by a large community of WDN modelers. All this work represents a preliminary step before the incorporation of new processes such as valves, pumps, and pressure-dependent background leakage outflows.

## HIGHLIGHTS

A new simulator of pressure-dependent consumption in water distribution network.

Coded in Python, based on a state-of-the-art MATLAB model, easy to extend.

Numerical experiments on networks composed of up to 19,647 pipes and 17,986 nodes.

Instabilities handled through regularization and line search along the Newton descent.

Preliminary step for integration of new processes and first contribution to OOPNET.

## INTRODUCTION

### Users’ consumption in WDNs

Households, schools, hospitals, businesses, food-processing and pharmaceutical industries, fire departments, etc., are all significant users of water distribution networks (WDNs). Each of them has its own need of water; moreover, this need can vary over time. Over a specific period, a need corresponds to a “demand”, expressed as an outflow rate. From this demand, and depending on the capacities of the WDN, the amount of water a user actually consumes is called its “consumption”. A user’s consumption is between 0 (i.e., no water is provided to the user) and the demand (i.e., the demand of the user is fully satisfied). For example, in the last three decades, the worldwide amount of drinking water consumed daily by households and public services is estimated to be , which represents 510,995 Olympic-size swimming pools of , and 11 percent of total freshwater withdrawals (Ritchie & Roser 2017). In France, this amount is of per day, which leads to an average domestic consumption of 1481 per day and inhabitant (Lao *et al.* 2022).

The first goal of WDNs is to satisfy the demands of the users. To fulfill these needs, WDNs must be well sized at their building, and then extended according to new needs (e.g., population growing or installation of new facilities). Pumps are used to maintain enough pressure to carry water from treatment plant or natural source to users. Also, valves are used to control the pressure, and to divide large networks in several district-metered areas (DMAs) that are easier to monitor with the help of sensors. Appropriate design and equipment permit us to limit the construction and functioning costs, by reducing, for example, the initial number of pipes, maintenance operations, energy consumption of pumps and valves, number of installed sensors, and amount of water losses due to leakages (Alperovits & Shamir 1977; Gupta *et al.* 1993; Walski *et al.* 2007; Gupta *et al.* 2016; Creaco Walski 2017). To optimize the choice of decision-makers, the use of pressure-driven models of users’ consumption is mandatory (Karadirek *et al.* 2012; Berardi *et al.* 2015; Laucelli *et al.* 2017a,b; Monsef *et al.* 2018; Koşucu Demirel 2022).

### Pressure-driven modeling of users’ consumption

One approach to model users’ consumption is by assuming that the demand of the users is always satisfied. These models are called “demand-driven models” (DDMs); they neglect unsatisfied users’ demand due to insufficient pressure. Nevertheless, this solution is implemented in several softwares, such as EPANET 2.0 (Rossman 2000), its object-oriented C++ interface OOTEN (van Zyl *et al.* 2003), and Porteau (Piller *et al.* 2011).

Conversely, the “pressure-driven models” (PDMs) permit us to simulate users’ consumptions that depend on the pressure. Different pressure–outflow relationships (PORs) can be used in this kind of model (Wagner *et al.* 1988; Fujiwara & Ganesharajah 1993; Gupta & Bhave 1996; Piller *et al.* 2003; Tanyimboh Templeman 2004). Among them, the Wagner’s POR (Wagner *et al.* 1988) was proven to perform best against data (Shirzad *et al.* 2013). Simulators using a POR to model PDM consumptions include the extension of EPANET 2.0 from Siew & Tanyimboh (2012), EPANET 2.2 (Rossman *et al.* 2020), WNTR (Water Network Tool for Resilience) (Klisel *et al.* 2018), and the MATLAB simulator developed by Elhay *et al.* (2016).

Pressure-driven analyses (PDAs) can also be conducted using DDM simulators by adding artificial equipment to the network (e.g., flow and pressure-regulating valves) (Mahmoud *et al.* 2017; Vaidya *et al.* 2023). However, this approach increases the size (i.e., the number of elements) of the system, and can require specific treatments and convergence criteria to ensure the robustness and the stability of the algorithm for every test case (Gorev *et al.* 2021; Sivakumar *et al.* 2023).

The engines and the programmer’s toolkit of EPANET 2.2 are multi-platform, but they are non objected-oriented; also, EPANET 2.2’s graphical user interface, which makes its engines easier to use, works only on Microsoft Windows operating systems. The WNTR package (Klisel *et al.* 2018) was designed specifically to simulate and analyze resilience of water distribution networks under disaster scenarios; moreover, its software architecture does not permit us to extend it easily, in particular to add new physical processes like pressure-dependent background leakage outflows. Finally, the MATLAB simulator from Elhay *et al.* (2016) is efficient and works fine, but any upgrade in its code will need to be run in MATLAB, a paid-for software package not available to everyone (although interfacing with other languages is possible). In summary, there currently exists no open-source and easy/cheap to extend simulator of PDM consumption.

### The Python language

Python is a free, multi-platform, high-level and object-oriented language (Lutz 2013). Its core philosophy emphasizes readability, making the codes written in Python easy to reuse and extend, even by people who are not software developers. Many efficient libraries for scientific computing are available in Python (Hunter 2007; Chen *et al.* 2008; Hagberg *et al.* 2008; McKinney 2011; Virtanen *et al.* 2020), and Python is today one of the most popular programing languages^{1} .

Big efforts have already been made by the community to simulate WDNs with Python. For example, the software of Xing Sela (2020) allows the running of transient simulations using an adapted method of characteristics. For extended-period simulations (EPS), the EPANET-based package WNTR (Klisel *et al.* 2018), and the object-oriented framework OOPNET (Steffelbauer & Fuchs-Hanusch 2015) are stable and active solutions.

### The framework OOPNET

OOPNET (Object-Oriented Python framework for water distribution NETworks analyzes) is a convenient Python API (Application Programing Interface) of EPANET. It is aimed at students, engineers and researchers who want to run hydraulic simulations in a high-level way, apply pre and post-processing on input or output data of EPANET, and make the simulation results from EPANET more readable (Steffelbauer & Fuchs-Hanusch 2015).

OOPNET is based on several state-of-the-art Python libraries for scientific computing (Hunter 2007; Hagberg *et al.* 2008; McKinney 2011; Hold-Geoffroy *et al.* 2014; Hoyer & Hamman 2017). It permits us in particular to parse and convert EPANET’s input file to Python objects, run an EPANET simulation through a command line interface, and parse EPANET’s output files and convert them to Python objects. OOPNET is free, open-source, object-oriented, and designed to enhance collaborative effort from the community of programmers and users of WDN software. Its source code is publicly accessible on GitHub (https://github.com/oopnet/oopnet). Moreover, its authors and contributors currently work on making OOPNET more flexible, user-focused and efficient through a domain-driven design (DDD) approach (Steffelbauer *et al.* 2022). Thus, it is well suited for prototyping, integrating new physical processes, or testing new hypotheses.

### Hypothesis, objectives, and research strategy

The consumption model proposed by Elhay *et al.* (2016) seems the most effective one; however, this model is not easy to extend and it needs a MATLAB environment. Thus, our objective is to implement a new simulator that will later permit us to test new physical processes in a more convenient way and to share our developments with a wider community of users and programmers.

We think that using the Python programing language and the framework OOPNET is today the best solution for developing such a simulator. Doing so, we expect to reproduce the results obtained by Elhay *et al.* (2016) with their MATLAB implementation when running PDA, and the ones obtained with EPANET 2.0 for demand-driven analyses (DDAs). It is important to note here that EPANET 2.2 was not available yet at the start of our study; this is the reason why we decided to validate our Python simulator against the MATLAB implementation of Elhay *et al.* (2016) rather than EPANET 2.2 when running PDA. Finally, developing this new Python simulator would permit us to extend in a non-regressive way the current capabilities of OOPNET, which by now allows running EPANET simulations only.

In this paper, we will first describe the steady-state equations and the method to solve them, presenting several numerical enhancements to deal with potential sources of instability. Then, we will explain our Python implementation, the networks to simulate, and the tests and metrics proposed for validation. Finally, we will present our results and discuss them.

## METHODS

Hereafter, since we do not model any other type of demand and consumption, we will call “user’s demand” just “demand”, and “user’s consumption” simply “consumption”.

*et al.*2007, p. 29). and , as well as friction head-losses, will be expressed in mHO (metres water column), which are consistent to m. Also, as in Piller

*et al.*(2011), all flow rates, demand and consumption will be expressed in rather than in , to avoid problems of stability due to machine precision.

Finally, and in the absence of other indications, scalar parameters and variables will be denoted in italic (e.g., ), vectors in italic bold (e.g., ), matrices in italic bold upper-case (e.g., ), scalar functions in upright (e.g., ), vector functions in upright bold (e.g., ), matrix functions in upright bold upper-case (e.g., ), and sets in blackboard style (e.g., ).

### Friction head-loss

*et al.*1933) to compute the friction head-loss. To do so, denoting the inside diameter (in mm) of the pipe and its roughness coefficient (unit-less), and supposing that the pipe is cylindrical, we first compute the friction coefficient as:Then, we compute the Hazen–Williams friction head-loss along the full length (in m) of the pipe as:where is the algebraic flow rate traversing the pipe, and is the Hazen-Williams exponent (unit-less).

### Consumption

*et al.*1988) to compute the consumption. To do so, denoting the demand, the pressure-head, and and the fixed minimum and service pressure-heads (in mHO) such that , we compute the fraction of pressure-head (unit-less) as:and the consumption as:Usually, minimum and service pressure-heads are chosen respectively equal to and (Piller

*et al.*2003; Giustolisi

*et al.*2008, 2014; Elhay

*et al.*2016; Deuerlein

*et al.*2019).

### Equations of equilibrium in a WDN

**ℓ**, the number of junction nodes at which the heads are unknown, and respectively the number of tank and reservoir nodes at which the heads are known and fixed, the total number of source nodes, and the total number of nodes. Also, we denote and the unknown flow rates in pipes and heads at junctions, and the known and fixed heads at tanks and at reservoirs, and the heads at all source nodes and at all nodes, and the pressure-heads at junctions and the water levels at tanks, and , and the elevations of junctions, tanks’ bottom, and reservoirs’ surface, respectively. Then, we have the relations:at junctions,at tanks, andat reservoirs. Next, we denote the junction-pipe, the tank-pipe, and the reservoir-pipe incidence matrices as respectively , and , the source-pipe incidence matrix as , the node-pipe incidence matrix as , and the vector function of to to compute the friction head-losses in pipes. , we compute the friction head-loss along the full pipe from Equation (3), as:Finally, we denote the fixed demand at junctions, and the vector function of to to compute consumptions. , is computed from Equation (5) as:

In this study, for the sake of simplicity, we choose the same for all junctions. However, it presents no difficulty to set specific values for some junctions, to describe for example pressure-independent fire demands as in Sivakumar *et al.* (2023).

*et al.*2007, p. 29). Also, the flow rates satisfy the conservation of the mass at nodes, derived from the Kirchhoff’s first law: “the algebraic sum of the branch currents towards any node of an electric network is zero” (IEC 2023).

### Solution algorithm

We solve the system (11) with a Newton’s method that includes a system reduction.

#### Newton’s method

the initial guesses of flow rates in pipes, and the ones of heads at junctions,

and the Jacobian matrix of at iteration ,

#### System reduction

#### Initial guesses

*et al.*(2016) and , the initial guess of the flow rate in the pipe is chosen consistent with a velocity :Also, , we set the initial guess of the head at as:using the same notation as in Sections 2.2 and 2.3.

#### Convergence criterion

*et al.*(2016). Indeed, we consider that the convergence is reached as soon as, ,The only difference with Elhay

*et al.*(2016) is that our version of the criterion also handles the case where . We discuss briefly the choice of the criterion in the Supplementary Material.

Other convergence criteria are used in WDN simulators. For example, the one in EPANET 2.0 checks that the 1-norm of flow changes divided by the 1-norm of flows is less than a given tolerance (Rossman 2000, p. 94). EPANET 2.2 improves this criterion by considering also the -norm of flow changes and the -norm of energy residuals (Rossman *et al.* 2020, p. 116). Then, Gorev *et al.* (2022) enhances EPANET 2.2’s criterion by checking for flow residual in each pipe too. Finally, Sivakumar *et al.* (2023) add one more criterion that verifies the residuals on the inverse of the POR (e.g., the Wagner’s POR) to improve the convergence of EPANET 2.2’s solution algorithm. In our study, we don’t need the criterion proposed by Sivakumar *et al.* (2023) because our formulation does not consider the POR inverse function. Therefore, we choose to simply reuse the criterion from Elhay *et al.* (2016), since preliminary tests showed good performances, and to be consistent with the other technical choices made by Elhay *et al.* (2016).

### Sources of instabilities

Different sources of instabilities can hinder the convergence of the solution algorithm described in Section 2.4. We present them below, along with numerical enhancements to deal with them. These numerical enhancements are of primary importance to ensure the convergence of the solution algorithm in many common situations. However, to preserve readability we choose to describe them in the Supplementary Material.

The numerical enhancements presented in Sections 2.5.1 – 2.5.3 are the same as the ones implemented in the Porteau software (Piller *et al.* 2011).

#### Pipes with zero flow rate

#### Junction nodes with pressure-head close to the minimum or the service pressure-head

*et al.*(2003), to regularize the consumption and its derivative by respectively a cubic and a quadratic polynomial, when is close to and .

#### Initial guesses far from the solution and/or objective function with (sub-)linear gradient growth

The bottleneck of the solution algorithm is the solving of the system (22). But the algorithm can need numerous iterations (and therefore numerous solving of system (22)) if the initial guesses of the flow rates and heads at junctions are far from the solutions at equilibrium. In this case, the simulation of large WDNs can become unfeasible in reasonable computational time. Moreover, the convergence of the algorithm is no longer guaranteed if the gradient of the objective function (11) to minimize does not have a super-linear growth (Piller *et al.* 2017). Thus, to reduce the number of iterations needed by the solution algorithm and to guarantee its convergence for any network configuration, we choose to reuse the damping method proposed by Elhay *et al.* (2016) (see description of the method in the Supplementary Material).

Some formulations of Equation (11) do not need damping corrections because they consider the inverse of the POR function rather than the POR function (Todini *et al.* 2021). However, to achieve full convergence, these formulations require a supplementary criterion on the inverse POR equation (Sivakumar *et al.* 2023).

#### WDN with highly contrasted values of flow rates and/or heads

In general, at each iteration of a solution algorithm based on Newton’s method, the magnitude of the elements in the Jacobian matrix can vary strongly from one part of a network to another, according to demands, pipe resistances, tank levels, use of pumps and/or valves, etc. In the present study, when demands and/or pipe resistances differ by several orders of magnitude across the network, using the Jacobian matrix as it can cause instabilities and increase significantly the number of iterations needed to converge. To limit the difference of magnitude order between all elements of , we choose to extend the preconditioning method initially proposed by Elhay & Simpson (2011). Indeed, in Elhay & Simpson (2011), the method dealt with consumptions that do not depend on the pressures; thus, we extend it to now deal with pressure-dependent consumptions too. Using this preconditioning method, the solution algorithm becomes a quasi-Newton method.

Gorev *et al.* (2022) propose an alternative solution to deal with zero flows and low-resistance pipes in EPANET 2.2. To do so, they replace the head-loss by a linear function in the pipes where the flow rate becomes very small, avoiding the appearance of zero head-loss derivatives in the Jacobian matrix. This regularization approach preserves the convergence rate of the Newton method, which is generally better than that of quasi-Newton methods. Nevertheless, it also changes slightly the head-loss behavior (for low flow rates), it uses an absolute threshold flow rate for the cutoff of the derivatives (while the preconditioning method from Elhay & Simpson (2011) is based on the difference of magnitude order between the elements of the Jacobian), and it needs some additional checking computations. Thus, both the Elhay & Simpson (2011) and Gorev *et al.* (2022) approaches have their own advantages and limitations. We choose to use the one from Elhay & Simpson (2011) because we master it better.

### Implementation and framework

To implement our Python simulator, we use the Python programing language^{2} through its Anaconda distribution^{3}. In addition to the Python’s standard library^{4}, we use several third-party libraries:

NetworkX (Hagberg

*et al.*2008), for the creation and manipulation of the graphed data structures of the WDNs,NumPy (Harris

*et al.*2020), for the creation and manipulation of the vectors and basic linear algebra operations (e.g., computation of vector norms),SciPy (Virtanen

*et al.*2020), to find the coefficients of the polynomials defined in the Supplementary Material^{5}, and to handle sparse matrices^{6},Scikit-Sparse

^{7}, which itself wraps the CHOLMOD routine (Chen*et al.*2008), to compute the Cholesky decomposition of the Schur complement matrix,XArray (Hoyer & Hamman 2017), to store the outputs in multi-dimensional labeled arrays,

Pandas (McKinney 2011), to convert the outputs to tabular data structures, to post-process them, and to write them to comma-separated values (CSV) files,

Matplotlib (Hunter 2007), to plot the outputs,

and Sphinx

^{8}, to extract the documentation from the source code, and generate an HTML file from it.

*simulators*(https://github.com/oopnet/oopnet/tree/simulators); its integration into the branch

*main*is in progress. We use OOPNET in particular to parse and convert EPANET’s input and output files to Python objects.

### Simulated networks

^{9}; they are plotted on Figure 1. Other networks can not be plotted here, because of ownership and/or security concerns.

Network . | . | . | . |
---|---|---|---|

934 | 848 | 8 | |

1,118 | 1,039 | 2 | |

1,976 | 1,770 | 4 | |

2,465 | 1,890 | 3 | |

2,508 | 2,443 | 2 | |

8,584 | 8,392 | 2 | |

14,830 | 12,523 | 7 | |

19,647 | 17,971 | 15 |

Network . | . | . | . |
---|---|---|---|

934 | 848 | 8 | |

1,118 | 1,039 | 2 | |

1,976 | 1,770 | 4 | |

2,465 | 1,890 | 3 | |

2,508 | 2,443 | 2 | |

8,584 | 8,392 | 2 | |

14,830 | 12,523 | 7 | |

19,647 | 17,971 | 15 |

### Tests and metrics

#### Functioning and performances of the Python simulator

From the outputs of the simulations, we first check that the -norm of the residuals is less than (in if the greatest absolute value is a mass residual, or in mHO if it is an energy one). If so, we then assume that the solver reaches the equilibrium state with enough precision.

*et al.*2008), we expect that the times needed by the solver increase slower than those needed for the whole simulations when larger networks are simulated. To check this point, we compute, for each network , the ratio of time needed by the solver when compared to the times needed by the whole simulations, defined as:We implement several numerical enhancements to deal with potential sources of instabilities (see Section 2.5). To quantify the gain of these enhancements, we simulate each network for each demand multiplier , first with no numerical enhancement, and next adding successively each numerical enhancement. Hereafter, we denote this set of enabled enhancements as (“” means “in addition to all previously enabled enhancement(s)”). We then compute, for each run, the convergence order of the solution algorithm when it converges. To do so, for each enabled enhancement(s) , each network and each demand multiplier , denoting the residuals at the iteration of the solution algorithm, we compute the difference between the -norm of and , , where is the number of iterations needed to reach convergence, as:The convergence order then satisfies the relation:with a constant. After rewriting Equation (39) in scale as:we then compute by linear regression of Equation (40).

#### Comparison of the Python simulator with the MATLAB one

To validate the ability of our Python simulator to run pressure-driven analyses (PDA), we compare the consumptions computed by our simulator with the ones obtained by Elhay *et al.* (2016) from their MATLAB implementation, for the same sets of networks and demand multipliers , and the same minimum and service pressure-heads at all junctions: and . Then, we assume that the consumptions from both Python and MATLAB simulators are the same as soon as they all differ by less than . Supposing that the average flow rate in the tested networks is equal to , then represents a precision of , which corresponds to 1/10th of the best precision expected when measuring flow rates physically in a WDN (Itron 2021, p. 3).

The MATLAB simulator of Elhay *et al.* (2016) implements the numerical enhancements described in Sections 2.5.1 to 2.5.3.

Finally, it is not possible to compare the process times elapsed using each simulator, because the MATLAB simulator was run by Elhay *et al.* (2016) on a different machine than the one we use to run our Python simulator.

#### Comparison of the Python simulator with EPANET 2.0

To check that our Python simulator is also able to simulate pressure-independent consumptions, we simulate the same set of networks with EPANET 2.0^{10} , increasing sufficiently the heads at the source nodes , and reducing enough the service pressure-head , so that the demand becomes satisfied at every junction. Then, we verify, at each junction of each network, that the consumption computed by the Python simulator and the one computed by EPANET 2.0 differ by less than .

Finally, to compare the performances of the Python simulator with those of EPANET 2.0, we compute, for both simulators and for each network, the -norms of the residuals at convergence, along with the CPU times elapsed running each simulator on the same machine: an Intel Core i9 with 32 GB of memory. For these comparisons to be relevant, we set the “ACCURACY” parameter of EPANET 2.0 to the same value as the one used in the convergence criterion of the Python simulator Equation (27), that is (unit-less). The “ACCURACY” parameter is used by EPANET 2.0’s solver to determine when the convergence is reached. Indeed, according to EPANET 2.0’s documentation (Rossman 2000, p. 153): “The trials end when the sum of all flow changes from the previous solution divided by the total flow in all links is less than this number’.

## RESULTS AND DISCUSSION

This section presents the results obtained when running the tests explained in Section 2.8.

### Functioning and performances of the Python simulator

When simulating the set of networks for each demand multiplier, the maximal -norm of all residuals is equal to (in or mHO), which is much less than . Thus, the solver reaches the equilibrium with enough precision.

Globally, these metrics show that for most of the configurations (i.e., one network combined with one demand multiplier), we need to regularize the friction head-loss to reach convergence. Also, friction head-loss regularization, consumption regularization and damping correction must be enabled for the simulation of all configurations to reach convergence. If we consider only the networks and demand multipliers for which all simulations converge, then the mean convergence order is better when none of the numerical enhancements are enabled (i.e., ). Even so, the mean convergence order remains very good (i.e., ) when a part or all of the numerical enhancements are enabled. Adding preconditioning when all other numerical enhancements are already enabled does not show any improvement in these test cases. However, this does not mean that there would be no difference for other networks and configurations. Thus, since we do not observe computational overhead, we decide, for future simulations, to systematically enable all numerical enhancements, including the preconditioning.

### Comparison of the Python simulator with the MATLAB one

When simulating each network for each demand multiplier, the maximal difference between the consumptions computed with the MATLAB and the Python simulators is equal to . This difference is not significant. Thus, both simulators give the same result, which means that our simulator is able to run pressure-driven analyses in a correct way.

For any network, the number of damping corrections is related to the complexity of the network: the higher it is, the more complex is the network. In our study, we could expect that this number increases with the size of the networks, i.e., that needs more damping corrections than , more than , etc. However, Figure 5(b) does not show any obvious proportionality relation between the size of the networks and the number of corrections. For example, even if they are smaller, networks and need many more corrections than or .

Independently to their size, another source of complexity in WDNs is their degree of meshing. Indeed, Piller (1995) showed, for pressure-independent consumptions, that heavily meshed networks need more operations to be solved. Adapting the work of Piller (1995) to the case of pressure-dependent consumptions, it is possible to demonstrate numerically that are more meshed than , and therefore more complex to solve.

Finally, strong variability of hydraulic resistances and/or demands across could also explain why these networks are more complex to solve than even if they are smaller. Indeed, the presence in a network of pipes with very low resistances together with pipes of very high resistances, and/or of junctions with very low demands together with junctions with very high demands, can increase strongly the stiffness of the Schur complement of the Jacobian matrix (defined by Equation (23)).

### Comparison of the Python simulator with EPANET 2.0

We find that increasing the heads at sources by and reducing the service pressure-head to leads to fully satisfied demand at every junction of every network. With such values of and , the maximal difference between the consumptions computed by the Python simulator and the ones computed by EPANET 2.0 is equal to . Thus, our Python simulator gives the same results as EPANET 2.0 when the demands are fully satisfied, which means that our simulator is able to correctly run demand-driven analyses.

EPANET 2.0’s convergence criterion has been improved in EPANET 2.2, by considering also the -norm of flow changes and the -norm of energy residuals (Rossman *et al.* 2020, p. 116).

As expected, the elapsed CPU times of the Python simulator are slightly greater than the ones of EPANET 2.0 (Figure 6(b)). This is normal because EPANET 2.0’s solution engine is coded entirely with the C programing language, which is compiled and thus more efficient than the Python language. But these differences remain small, maybe because the CHOLMOD (Chen *et al.* 2008) library used by the Python simulator to compute the Cholesky decomposition of the Schur complement matrix is more optimized than the method used in EPANET 2.0, which is based on an older work from George & Liu (1981); however, a rigorous comparison of these two methods would require some further work.

## CONCLUSIONS

In this paper, we designed and implemented a new Python simulator to model pressure-dependent users’ consumptions in WDNs with no need of artificial elements. Our simulator produces the same results as those computed by the state-of-the-art MATLAB simulator developed by Elhay *et al.* (2016). Extensive numerical experiments on complex, large and realistic networks, composed of up to 19,647 pipes and 17,986 nodes, showed similar precision and efficiency for both simulators. Moreover, the performances of our simulator are close to those of EPANET 2.0 in case of fully satisfied demand.

The community of WDN modelers now has at its disposal a new efficient and flexible simulator of pressure-dependent users’ consumption, integrated into the generic and collaborative framework OOPNET (Steffelbauer & Fuchs-Hanusch 2015). This is the first Python simulator that includes all numerical enhancements from Piller (1995), Piller *et al.* (2003), Elhay & Simpson (2011) and Elhay *et al.* (2016). These enhancements were mandatory for the simulation of all tested configurations to converge. Among these enhancements, the damping correction will allow the study of more complex models, especially those for which the gradient of the objective function to minimize can have sub-linear growth.

A complete set of tests were driven to assess the good functioning, stability and robustness of our simulator. Also, the native properties of the Python language make code extension easier. Therefore, our Python simulator can now be reused to study trickier physical processes, such as valves, pumps and pressure-dependent background leakage outflows, while limiting the risk of software regression.

An interesting further work could also consist in comparing the results and performances of our Python simulator with those of EPANET 2.2 (which was not available yet at the beginning of this study) when running pressure-driven analyses. In particular, since the decomposition of the Schur complement matrix is a bottleneck of both our simulator and EPANET 2.2’s solution algorithms, we would like to compare the efficiency of the CHOLMOD library (Chen *et al.* 2008) used in our simulator, with that of the Cholesky decomposition from George & Liu (1981) and implemented in EPANET 2.2. The comparison of these two implementations, both written in C programing language, should not present any difficulty since our simulator has a modular architecture that permits us to change easily the solver and library to use.

## AUTHOR CONTRIBUTIONS

C. C. was involved in conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing the original draft, writing the review, and editing.

O. P. was involved in conceptualization, formal analysis, methodology, supervision, validation, writing the review and editing.

I. M. was involved in conceptualization, supervision, writing , revising and editing.

## FUNDING

This work is part of the Oriented Renewal of Pipes (ROC) French research project and was supported by the Nouvelle Aquitaine region, the Loire-Bretagne and Adour-Garonne water agencies, the Aquitaine regional public health authorities (ARS), and the European Regional Development Fund (ERDF).

## DATA AVAILABILITY STATEMENT

Please refer to the online version of this article for figures in color. Additional relevant data to the manuscript are included in the Supplementary Materials. Source code of the Python simulator is publicly accessible on https://github.com/oopnet/oopnet/tree/simulators.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

Source code from https://www.epa.gov/sites/default/files/2018-10/en2source.zip.

## REFERENCES

*Exploring Network Structure, Dynamics, and Function Using Networkx*. Technical Report LA-UR-08-05495; LA-UR-08-5495, Los Alamos National Lab. (LANL), Los Alamos, NM (United States). Available from: https://www.osti.gov/biblio/960616

*Proceedings of the 2014 Annual Conference on Extreme Science and Engineering Discovery Environment*. Association for Computing Machinery, New York, NY, USA, pp. 1–8. doi:10.1145/2616498.2616565

*International Electrotechnical Vocabulary*. Available from: https://www.electropedia.org/

*Static Water Meter INTELIS*. Available from: https://www.itron.com/-/media/feature/products/documents/brochure/intelis-water-static-meter-en-web.pdf

*1st International WDSA/CCWI2018 Joint Conference*. Kingston, Ontario, Canada, pp. 1–8

*National Report of SISPEA Data*. Technical Report Edition of Jan. 2022 from data of 2020, French national observatory of water and sanitation services. Available from: https://www.services.eaufrance.fr/docs/synthese/rapports/Rapport_Sispea_2020_VF.pdf

*2017 International Conference on ENERGY and ENVIRONMENT (CIEM)*. pp. 241–245

*Modeling the Behavior of a Network - Hydraulic Analysis and a Sampling Procedure for Estimating the Parameters*. PhD Thesis in Applied Mathematics, University of Bordeaux, France. Available from: http://www.theses.fr/1995BOR10511

*World Water & Environmental Resources Congress 2003*. American Society of Civil Engineers, Philadelphia, Pennsylvania, United States, pp. 1–15. doi:10.1061/40685%282003%29113

*Eleventh International Conference on Computing and Control for the Water Industry (CCWI 2011)*. Centre for Water Systems, University of Exeter, Exeter, United Kingdom, pp. 27–32. Available from: https://hal.archives-ouvertes.fr/hal-00653763

*Computing and Control for the Water Industry (CCWI) 2017*. The University of Sheffield, p. 327189 Bytes. Available from: https://figshare.shef.ac.uk/articles/journal_contribution/CCWI2017_F41_Why_are_Line_Search_Methods_Needed_for_Hydraulic_DDM_and_PDM_Solvers_/5364229/1

*Water Use and Stress. Our World in Data*. Available from: https://ourworldindata.org/water-use-stress

*14th International Conference on Hydroinformatics, Water INFLUENCE Water INFormatic soLUtions and opEN problems in the cycle from Clouds to ocEan*. Bucharest, Romania, p. 4. Available from: https://hal.inrae.fr/hal-03750370

*The Fourth International Conference on Engineering Computational Technology*. Lisbon, Portugal, p. 64. Available from: http://www.ctresources.info/ccp/paper.html?id=3699

*Conference: CCWI 2003 Advances in water supply management*. Imperial College, London, United Kingdom, pp. 1–8