A Head Formulation for the Steady-State Analysis of Water Distribution Systems Using an Explicit and Exact Expression of the Colebrook–White Equation

Steady-state demand-driven water distribution system (WDS) solution is the bedrock for much research conducted in the field related to WDSs. WDSs are modeled using the Darcy–Weisbach equation with the Swamee–Jain equation. However, the Swamee–Jain equation approximates the Colebrook–White equation, errors of which are within 1% for e/D ∈ [10−6, 10−2] and Re ∈ [5000, 108]. A formulation is presented for the solution of WDSs using the Colebrook–White equation. The correctness and efficacy of the head formulation have been demonstrated by applying it to six WDSs with the number of pipes ranges from 454 to 157,044 and the number of nodes ranges from 443 to 150,630. The addition of a physically and fundamentally more accurate WDS solution method can improve the quality of the results achieved in both academic research and industrial application, such as contamination source identification, water hammer analysis, WDS network calibration, sensor placement, and least-cost design and operation of WDSs.


Introduction
Water distribution systems (WDSs) are essential infrastructures of every city and town, the purpose of which is to satisfy the water requirements of the population, of agriculture and for industry with the required quality and quantity. The hydraulic steadystate solutions of WDSs (the flows and heads) is the foundation of many, if not all, WDS academic research and industry application. Therefore, the speed and accuracy of the hydraulic simulation model that is used to find the steady-state of WDSs underpin the quality of the research outputs and industry applications. The quest for a solution method for finding the steady-state solution of a looped WDS can date back to 1936. Since then, the research community diverged into three main branches: (1) loop-based methods, (2) null-space method, and (3) range space methods.
The loop-based methods use the loop energy equations and continuity equations to model the demand-driven steady-state of WDSs. The first loop-based methods and the first WDS solution method is the Hardy-Cross method [1]. The Hardy-Cross method is an iterative manual method that uses successive approximation to solve the above nonlinear system of equations in which a set of initial flows that satisfies the mass conservation equations is successively corrected until the stopping test has been met. The loop identification and the requirement that the initial guess of flows must satisfy continuity are the factors that affect the performance of the Hardy-Cross method. The performance of the loop-based methods is dependent on the sum of the length of all identified loop. The Shortest Cycle Basis [2] is the best set of loops that minimizes the time that is required to execute any loop-based method. This is because the use of the shortest cycle basis can achieve the minimum number of non-zeros in the key matrix. Ref. [3] explored the use of minimum loops in the Newton-Raphson loop flows method and showed the time used to identify the shortest cycle basis is 100 to 10,000 times the time used to solve the network. Ref. [4] proposed two algorithms to select a set of network loops to achieve a highly sparse matrix. Although a smaller number of non-zeros in the Schur complement was reported in Creaco and Franchini [3], the substantial improvement in terms of the efficiency reported by Alvarruiz et al. [4] suggests the latter algorithm is the better practical choice. The overhead to identify the shortest cycle basis is clearly the bottleneck in the loop-based method even though the loop identification algorithm is only required to execute once for a given network topology.
The null-space methods, methods that operate in the subspace defined by a null space that is orthogonal to the column of the unknown-head node-arc incidences matrix, partitions the network into a spanning tree acyclic graph and the complementary chord tree edges. The addition of any chord tree pipe to the spanning tree graph creates a loop. The null-space method, in the context of WDS solution methods, is a special case of the loop-based methods. This is because the cycle basis created by any spanning tree graph and chord edges is a subset of the set of all cycle bases of the loop-based method, whereas the cycle basis of the loop-based method (particularly the shortest cycle basis) cannot be expressed as a combination of a spanning tree graph and the complementary chord tree edges. By sacrificing the generality, the null-space method requires fewer computation resources to find a combination of spanning tree graph and chord tree edges. Rahal [5] proposed a co-tree flows formulation in which it is necessary to (1) identify the associated circulating graph; (2) determine the demands that are to be carried by the spanning tree branches; (3) find the associated chain of branches closing a circuit for each co-tree chord; and (4) compute pseudo link head losses. Later, Elhay et al. [6] exploits the relationship between the co-tree flows and spanning tree flow by applying the Schilders' factorization [7] to permute the A1 matrix into a lower triangular square block at the top, representing a spanning tree, and a rectangular block below, representing the corresponding co-tree.
The range space methods include a collection of methods that operate in the subspace defined by the rows of the unknown-head node-arc incidences matrix. The global gradient algorithm [8] is the most widely used WDS solution method, which solves for the flows and the heads in WDSs simultaneously by exploiting the block structure of the Jacobian matrix to reduce the size of the key matrix in the linearization of the Newton method. The graph matrix partitioning algorithm (GMPA) [9] exploits the linear relationship between the flows in the internal trees and the flows in the superlinks to speed up the solution process of the GGA.
In addition to the different solution methods developed for the simulation of the WDS steady-state, network partitioning using graph theory has been another active research avenue. Network partitioning identifies sections of a WDS network that are hydraulically independent. The first graph structure that is exploited is the forest component of a WDS. A tree in a graph is any two vertices are connected by exactly one edge. Most WDSs have trees, the collections of which are called forests. Simpson et al. [10] proposed a forest-core partitioning algorithm to partition a WDS graph into its linear forest component and nonlinear core component. The flows in the forest pipes can be computed a priori and the heads in the forest nodes can be computed a posteriori by a linear process. Qiu et al. [11] proposed a bridge-block partitioning algorithm to partition a WDS graph into several linear independent blocks and bridges. The steady-state solution of one block of a WDS can be computed independently from other WDS blocks.
EPANET2 [12] is one of the most widely used WDS simulation packages, in which the GGA is used to provide steady-state demand-driven solution of WDSs. The code for EPANET 2 is in the public domain, allowing many studies to be conducted. These include the least-cost design and operation of WDSs, sensor placement, chlorine decay models calibration, contamination event detection, network vulnerability analysis, cyber-attack detection, network decontamination, and many others. Any inaccuracies that existed in the WDS solution method used will be inherited and sometimes exacerbated. One of the accuracy-related problems is the use of the Swamee-Jain equation [13], an approximation of the implicit Colebrook-White equation [14], a general equation that can be used when the Reynolds number is greater than or equal to 4000, to calculate the friction factor. It is reported in Swamee and Jain [13] that the errors involved in friction factor are within 1% for /D ∈ [10 −6 , 10 −2 ] and when Re ∈ [5000, 10 8 ]. All WDS solution methods that are currently available (1) ignored the above 1% error and (2) extended the applicability of the Swamee-Jain equation to calculate the friction factor when Re ≥ 4000 and all values of /D. This is mainly due to the solutions of the inexplicit Colebrook-White equation are required many times for multiple iterations, which is time-consuming.
In this paper, an iterative solution method that finds the steady-state demand-driven WDS solution using the Colebrook-White equation for the turbulent flow regime is proposed. This is achieved by (1) expressing the Colebrook-White equation, an implicit function, as an explicit, exact, and differentiable function to describe the friction factor in the turbulent flow regime, (2) using the Hagen-Poiseuille equation to describe the friction factor in the laminar flow regime, and (3) using cubic interpolation to fit a curve between the turbulent and the laminar flow regimes to describe the transitional flow regime. It is important to note that unlike the Swamee-Jain equation, the explicit expression of the Colebrook-White equation agree at all points, the use of which eliminates the inaccuracy that is associated with the Swamee-Jain equation. It is shown in this paper that the different steady-state solutions of WDSs are observed which can be a critical problem for some research areas (such as the least-cost design and operation, contamination event detection, network vulnerability analysis, sensor placement, district meter area, etc.). This paper is organized as follows. A description of the general WDS demanddriven steady-state problem is given in the next section. This is followed by Section 3 the derivation of the proposed head formulation. Section 4 gives an algorithmic description of the proposed head formulation, followed by the validation of the proposed friction factor equation against the Colebrook-White equation. Six case study networks are then described in Section 6, the results of which are discussed in the next section. The last section offers some conclusions.

Definitions and Notation
Consider a water distribution system that contains n p pipes, n j junctions, and n r fixed-head nodes.
The i-th node of the network has two properties: its nodal demand d i and its The matrix A 1 is the full-rank, unknown-head, node-arc incidence matrix. The matrix A 2 is the fixed-head node-arc incidence matrix. The head loss exponent n is assumed to be dependent only on the head loss model: n = 2 for the Darcy-Weisbach head loss model and n = 1.852 for the Hazen-Williams head loss model. The head loss within the pipe p is modeled by h f p = r p q p |q p | n−1 . Denote by G(q) ∈ R n p ×n p , a diagonal square matrix with elements [G] pp = r p |q p | n−1 for j = 1, 2, ....n p . Denote by F(q) ∈ R n p ×n p , a diagonal square matrix where the p-th element on its diagonal [F] pp = d dq p [G] pp q p .

System of Equations
The steady-state flows and heads in a WDS system are modeled by the demand-driven model (DDM) continuity Equation (1) and the energy conservation Equation (2): which can be expressed as where its Jacobian matrix used in the solution process is and it is sometimes referred to as a nonlinear saddle point problem [15]. This nonlinear system is often solved by the Newton method, in which q (m+1) and h (m+1) are repeatedly computed from q (m) and h (m) by until the relative differences are sufficiently small.

Derivation of the Head Formulation
Consider the continuity equations in Equation (1) and a vector of unknown heads, h, for the n j nodes in a network.
Assume the flows in pipe p can be expressed as a function of the head loss in pipe p: Let a * ,j ∈ R n p ×1 denote the j-th column of the A 1 matrix, the continuity equation can be rewritten as: is the mass balance equation expressed for every node in a WDS.
The partial derivative of f c (q 1 , q 2 , . . . , q n p ) with respect to h f can be expressed as: in which the partial derivative of f c j (q 1 , q 2 , . . . , q n p ) with respect to h f p can be expressed as: As a result, Equation (7) can be written as: where the first matrix of the matrix multiplication above can be written as ∂A T 1 q/∂q, which is A T 1 , and the partial derivative of q p (h f p ) (Equation (6)) with respect to h j can be expressed as: Substitute Equation (10) into Equation (9), we get: where the third matrix in the matrix multiplication above can be expressed as (∂A 1 h + A 2 e l )/∂h, which is A 1 . Denote by W ∈ R n p ×n p , a diagonal square matrix where the p-th element on its diagonal [W ] pp = f (h f p ). The Jacobian matrix in Equation (11) can be expressed as: Omitting term ∂d ∂h for the demand-driven analysis, this new WDS solution method can be solved by the Newton method in which dh (m+1) are repeatedly computed from h (m) using: Moving h (m) to the right-hand-side of Equation (13), we get: Equation (6) is now derived for Hazen-Williams and Darcy-Weisbach head loss models in the next subsections. (15) and can be rewritten as:

The Hazen-Williams Equation head loss equation is
and the partial derivative of q p (h f p ) (Equation (16)) with respect to h f p can be expressed as:

Darcy-Weisbach Equation
The Darcy-Weisbach head loss equation is: Equation (18) can be rearranged into The friction factor for turbulent flow regime (Re ≥ 4000) can be calculated by using the Colebrook-White equation: (20) in which the Reynolds number is Substitute Equations (19) and (21) into the RHS of Equation (20), we get: Denote by S p = h f p /L p , the hydraulic gradient in pipe p. The friction factor in pipe p can be expressed as: The friction factor for laminar flow regime (Re ≤ 2000) can be calculated by using the Hagen-Poiseuille equation: Substitute Equation (24) into Equation (18), we get: The friction factor for transitional flow regime (2000 ≤ Re ≤ 4000) can be calculated by using a cubic interpolation between the laminar flow regime and the turbulent flow regime using: is the friction factor for pipe p when the Reynolds number is 4000, and f p 4000 is the derivative of friction factor when the Reynolds number is 4000.
we can express: 0.25 ln 2 (10) , the flows in pipe p as a function of the head loss in pipe p (Equation (6)) for the Darcy-Weisbach head loss model can be expressed as where its derivative can be expressed as: and dS p dh f p = L −1 . It is important to note that special care must be taken to calculate the q (h f ) value for laminar flows as the pipe head loss will be cancelled when calculating the pipe flow: as a result, the derivative of the pipe flow with respect to the pipe head loss in the laminar flow regime can be expressed: which is a constant value. Therefore, a pipe with zero head loss will not cause any numerical problems as a result of the zero head loss. Finally, S p 2000 and S p 4000 need to be identified for pipe p for all pipes. S p 2000 can be easily identified as f 2000 = 64/2000 = 0.032 for any pipes and q p 2000 = 500 πD p ν. Therefore, The determination of S p 4000 involves the solution of f p 4000 using the implicit Colebrook-White equation for pipe p when Re = 4000 and substitute the value of f p 4000 in Equation (2). This is only time when the direct solution of the implicit Colebrook-White Equation is required.

Head Formulation Algorithm
The steps of the proposed head formulation are described in Algorithm 1. The proposed head formulation, a single-phase formulation, within the iterative phase (between lines 4 and 18 in Algorithm 1 is similar in terms of the computational intensive when compared to the global gradient algorithm, a two-phase formulation.
Meanwhile, the overhead of the proposed head formulation, particularly the computation of the S p 4000 that requires the solution of the Colebrook-White equation for the Darcy-Weisbach head loss model, can increase the computation burden of the algorithm. However, as the manufacturing limitation, there is only a limit number of pipes that is available for the construction of WDSs. Therefore, the number of combinations of the pipe roughness heights and pipe diameters is limited, so does the number of distinct values of S p 4000 . This is also true for the least-cost design problem of a WDS. S p 4000 value is only required to be computed and stored for each of the commercially available pipes a priori. This stored value can then be retrieved during the optimization phase.
Moreover, Equation (16) is just an inverse function of Equation (15). As a result, the dq dh f required in the proposed head formulation is the inverse of the dh f dq required in the GGA. Therefore, the Schur Complement in the GGA, A 1 T F −1 A 1 , is the same as the Jacobian matrix shown in Equation (12) as W = F −1 . Thus, if the proposed head formulation and the GGA starts with the same initial starting points when applied to a WDS with Hazen-Williams head loss model, they will also have the same iterative solutions and the same final solution. Calculate q (i−1) using Equation (16); 12 Calculate (17); 13 endif 14 Calculate h (i) using Equation (14);

Validation of the Proposed Friction Factor Equation
The proposed friction factor equation is validated in this section. The proposed friction factor equation and the Swamee-Jain friction factor equation will be used to compute friction factors for pipes with the value of = (1 mm, 0.1 mm, 0.01 mm, 0.001 mm), D = 100 mm and Reynolds number ranging from 4000 to 10 8 , and the computed friction factor value in turbulent flow regime will be plotted against the value of the friction factor computed by using the Colebrook-White equation. As can be seen from Figure 1, the differences are observed between the friction factors computed by using the Swamee-Jain friction factor equation and that computed by using the Colebrook-White equation as shown in Swamee and Jain [13], whereas the friction factors computed by using the proposed friction factor equation are the same as that computed by using the Colebrook-White equation.

Case Studies
The proposed head formulation was implemented in WDSLib [16] and has been applied to six WDSs. The networks used here were assigned network identifiers from N 1 to N 6 (see details in Table 1): • Network identifier N 1 is assigned to Balerma Network, first introduced by Reca and Martínez [17], is comprised of 454 pipes, 443 junctions, and four reservoirs. • Network identifier N 2 is assigned to Richmond Network, first introduced by Van Zyl et al. [18], is comprised of 834 pipes, 848 junctions, and eight reservoirs. • Network identifier N 3 is assigned to exnet Network, first introduced by Farmani et al. [19], is comprised of 2465 pipes, 1890 junctions, and three reservoirs. This is important to note that valves are replaced by pipes. • Network identifier N 4 is assigned to the large Network used in Sitzenfrei et al. [20]. This network has 4021 pipes, 3557 junctions and one reservoir. • Network identifier N 5 is assigned to Network 2 of the Battle of Network Sensors competition [21]. Network N 5 is comprised of 14,830 pipes, 12,523 junctions, and seven reservoirs. This is important to note that valves and pumps are replaced by pipes and demand patterns are removed. • Network identifier N 6 is assigned to virtRome Network, first introduced by [20], is comprised of 157,044 pipes, 150,630 junctions, and four reservoirs.
All case studies were performed on an Intel Core i7-7700 running at 3.6 GHz with four cores in C++ under IEEE-standard double-precision floating arithmetic [22] with machine epsilon mach = 2.204 × 10 −16 .  Figure 2 shows the convergence of the proposed iterative head formulation. Networks N 1 and N 2 are both relatively small network, both of which have a small number of loops (11 loops for network N 1 and 14 loops for network N 2 ). The stopping test and the continuity residual have been met for both smaller size networks. However, convergence problem has been observed for the medium size networks (networks N 3 and N 4 ) and large size networks (networks N 5 and N 6 ). It is worth noting that the relatively head differences use in the stopping test has a better convergence property that the continuity residual. This is particularly pronounced in network N 2 as shown in Figure 2b. The stopping test has been met at seven iteration, whereas the continuity residual is 10 −3 . This significant difference is also observed in Figure 2c-f for networks N 3 -N 6 .

Results and Discussion
Using junction 776 and reservoir E in network N 3 as an example, the head iterates of junction 776 is observed to be oscillating around the elevation head of reservoir E as shown in Figure 3a. The relative head difference of the nodal head at junction 775 is 3.94 × 10 −10 at 6th iteration. However, this small perturbation in the nodal head caused the flow direction reversed. This flow direction reversal happened three more times while the value of the nodal heads converges to the true value. This is because the pipe (Pipe 1865) that connects junction 776 and reservoir E is 1 m in length and 0.999 m in diameter with 0.15 mm roughness height, which means the resistance factor of pipe 1865 is a very small value. This can also be seen from Figure 3a as the head loss between reservoir E and junction 776 is 10 −8 .  Once a damping factor of 0.67 has been applied, the head convergence of this particular node is more well-behaved as a faster convergence is achieved and no head oscillation has been observed as can be seen from Figure 3b. Please note that in Figure 3b and in Figure 4 the damped Newton's method was used. Within this scheme, the derivative within the Newton method is multiplied by a damping factor between zero and one for accelerating convergence. On the one hand, the application of the damping factor caused a slower head convergence of the node that is well-behaved before its application. This can be seen from Figure 4a,b as a slower rate of convergence is observed when compared to that from Figure 2a  The convergence properties, including the relative head differences, continuity residuals, energy residuals, and Colebrook-White equation residuals, at the final iteration of the proposed head formulation after applying the damping factor is shown in Table 2.  Table 3 presents the detailed timing results of pre-processing, iterative phase and post-processing operations of the GGA and proposed head formulation applied to the six case study networks. The total wall-clock time required to apply the proposed head formulation is higher than that required by the GGA for all six case study networks. This is because the proposed head formulation requires more time to perform the iterative phase while the GGA and the proposed head formulation require a similar time to execute the preprocessing and post-processing operations. In addition, a similar amount of per-iteration runtime is required by the GGA and the proposed head formulation, which means the longer wall-clock time required by the proposed head formulation is the result of higher number of iterations required due the damping factor applied. Table 3. The number of iterations and the wall-clock time (second) required to perform the preprocessing, iterative phase, and the post-processing operations for the GGA and the proposed head formulation applied to each of the six case study networks (the number in the bracket indicates the per-iteration run time required to execute the iterative phase of the Newton method).  Table 4 shows the detailed statistics of the absolute head differences between the GGA and the proposed head formulation applied to each of the six case study networks. This head differences in the results of optimization problems can cause constraints violation of the optimal solution identified. Take the Balerma network N 1 , one of the well explored network, as an example. This problem is first introduced by [17]. Two least-cost designs presented by [23,24], both of which using the GGA formulation to find the steady-state solution in which the Swamee-Jain equation is used to model the pipe friction factors, have found to have nodal pressure violations when using the proposed head formulation as shown in Table 5, in which the Colebrook-White equation is used to model the pipe friction factors. The above two least-cost designs of network N 3 are both infeasible as the Swamee-Jain equation approximates the Colebrook-White equation. This is because (1) the errors involved in Swamee-Jain friction factor are within 1% for /D ∈ [10 −6 , 10 −2 ] and Re ∈ [5000, 10 8 ] as reported in [13], whereas GGA has extended the applicability of the Swamee-Jain equation to all values of /D and Re ≥ 4000 and (2) the errors involved in Swamee-Jain friction factor have been ignored in all optimization problems. 3.00e-3 2.31e-2 1.61e-2 9.22e-2 1.62e-2 In addition to the differences in the head solutions, the use of friction factor produced by the Colebrook-White equation in the proposed head formulation also produces different flow results when compared to that produced by the GGA as shown in Table 6 and the spatial distribution of the pipes with different flows is shown in Figure 5. The flows of the pipes in the looped component found by using the proposed head formulation for each of the six case study networks are different from that found by using the GGA, whereas the flows of the pipes in the forest component found by using the head formulation for each of the six case study networks are the same as that found by using the GGA. In addition, the number of pipes with flow direction reversal is ranging from six pipes in network N 2 to 730 pipes in network N 6 .

Network Methods
In addition to the single-period steady-state WDS simulation where the boundary conditions (pumps and tanks) are fixed, errors caused by the difference between the Colebrook-White equation and the Swamee-Jain equation can accumulate as the boundary condition in the extended-period simulation. The error accumulation is manifested as the different tank level at each time step, the time when an operation activated by the trigger level, and the pump operating points. Net3 in EPANET is used as an example to demonstrate the error accumulation described above. It is important to note that the Net3, which is a simple WDS with 92 nodes, two reservoirs, three tanks, two pumps, and 117 pipes, in EPANET has been converted from the Hazen-William head loss model to the Darcy-Weisbach head loss model. The model difference between the GGA and the head formulation is relatively small between 00:00 to 05:00, which is also observed in network N 1 . As can be seen from Figure 6, however, the tank levels start to diverge as the tank level rises. This is because the insignificant model errors for small network can accumulate over a period of time in the form of boundary conditions. It is also worth noting that when the tank operation mode changes, from filling to discharging or from discharging to filling, the tank level differences between the two model start to narrowing, because the cancelling of errors, but never reaching zero. The maximum tank level difference is 0.074 m for tank 1, 0.12 m for tank 2 and 0.37 m for tank 3. Due to the error build-up in the boundary conditions, the differences between GGA and the head formulation flow results are relatively small with no flow direction reversal before the mid-day as can be seen from Figure 7a. After 12:00, however, significant differences between GGA and the head formulation flow results are observed with the reversal of flow direction as can be seen from Figure 7c-e. Finally, flow result differences start to decrease as the tank level differences narrow. The GGA and the head formulation start with the same tank levels and end up with different tank levels over the 24 h simulation period. This boundary condition differences will keep accumulating over a longer simulation period.

Possible Applications
The main improvement of using the proposed head formulation is the use of the Colebrook-White equation addresses a fundamental error that is associated with the use of Swamee-Jain equation, an approximation of the Colebrook-White equation, in all existing WDS solution methods. It is shown in the last section that both flow and head results of the head formulation are different from that of the global gradient algorithm. The differences in flow and head results can significantly affect several research areas.

Water Quality Simulation
Water quality is one of the most important water distribution system research subject. Water quality is simulated using the chemical transportation equation (advection) and the chemical decay equation, both of which are functions of pipe flows and time for both the contaminant and the disinfectant. As a result, the differences in flow results between the proposed head formulation and the GGA, particularly the flow reversal, can significantly affect the evaluation of the chemical transportation equation. This can be a critical problem for the identification of the contamination sources, network decontamination, disinfection by-product modeling, water quality sensor placement, and others.
In addition, the water quality modeling is temporal in nature, therefore, the flow results of the extended-period simulation are normally used in the chemical transportation equation. As can be seen from Figure 7 that the differences in flow results have no correlation and the pipes with flow reversal are not always the same. As a result, there is another level of error accumulation in water quality simulation on top of the extended-period simulation as opposed to the water quality simulation using the steady-state simulation where only the initial conditions are different. These two levels of error accumulation are compounded so that the errors in the water quality modeling can be amplified.

Water Hammer Analysis
Water hammer analysis is another research area that is affected by the flows, heads, and head losses in pipes and nodes of water distribution systems. The governing equations for the water hammer analysis are the unsteady momentum equation and the unsteady continuity equation. Both unsteady equations are functions of heads, flows, and the direction of flows. The reversal of the flow direction in a pipe will significantly alter the result of water hammer analysis in a looped water distribution system.

WDS Network Calibration
The proposed head formulation can also improve the quality of the WDS calibration models. Water distribution system calibration is a process of comparing model results with field data and making the appropriate adjustments so that both results agree, and it is usually applied to the estimation of pipe roughness values and nodal demands [25]. All existing calibration studies compared the 'measurement data', which is generated using EPANET instead of using the real measurement data, against the existing model results, both of which uses the Swamee-Jain equation to model the turbulent flow regime. As a result, the errors existed in the Swamee-Jain approximation is inherited by the calibrated water distribution system. A better network calibration model can be produced by using the proposed head formulation as the differences between the Colebrook-White equation and the experimental data are smaller than that between the Swamee-Jain equation and the experimental data.

WDS Optimization and Operation
The small head differences can invalidate all existing optimal solutions of water distribution system modeled by Darcy-Weisbach head loss equation. Examples of the pressure violations in the optimal solutions found in [23,24] for the Balerma network using the global gradient algorithm using Swamee-Jain equation are shown in Table 5, both of which become infeasible if the proposed Colebrook-White equation is used.

Conclusions
This paper presents an efficient iterative head formulation for the steady-state demanddriven solution of water distribution systems for both Darcy-Weisbach and Hazen-Williams head loss model. When the Hazen-Williams head loss model is used, the proposed solution method produces the same final and iterative flow and head solutions if the same initial guess is used as the global gradient algorithm. When the Darcy-Weisbach head loss model is used, an exact and explicit expression of the inexplicit Colebrook-White friction factor equation is proposed in this study. A cubic interpolation between this explicit expression of the inexplicit Colebrook-White equation for the turbulent flow regime and the Hagen-Poiseuille equation for the laminar flow regime is generated to describe the friction factor in the transitional flow regime.
The main features of the proposed head formulation of the steady-state demanddriven WDS simulation include: 1.
friction factor for the turbulent flow regime can be calculated using an explicit and exact expression of the Colebrook-White equation without the need for an iterative method; 2.
the use of the proposed head formulation can significantly improve the accuracy of the steady-state demand-driven solution of WDSs when compared to the GGA. This is because the proposed head formulation uses an explicit and exact expression of the Colebrook-White equation to calculate friction factor for the turbulent flow regime as opposed to the Swamee-Jain equation, an approximation of the Colebrook-White equation, used in previous WDS solution methods.
The efficacy of the proposed head formulation has been demonstrated by applying it to six case study networks, the results of which have been validated using the continuity residuals, energy residuals, and Colebrook-White residuals. It should also be observed that the proposed method could be selected for analytical (e.g., Di Nucci [26]) and for numerical applications (e.g., Pasculli [27]), and in particular for improving the Wall Function.Differences between the proposed head formulation and the GGA have been observed in both the flows and heads. On the one hand, the flow differences, particularly the flow direction reversal, can be a critical problem for some research areas such as water quality simulation, contamination source identification, water hammer analysis, WDS network calibration, sensor placement, and network clustering. On the other hand, the heads differences can cause pressure violations in the WDS least-cost design when the GGA is used to perform the steady-state demand-driven analysis. It is important to note that damping factor has been applied to the Newton method to achieve a more stable convergence. The choice of initial guess, another important factor in the Newton method convergence, and inclusion of valves, pumps, and other network elements will be interesting avenues for future research.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.