Applied Pipe Roughness Identification of Water Networks: Consideration of All Flow Regimes

This article presents a mathematically rigorous unification of the pipe roughness identification problem of water distribution networks considering all Reynolds regimes, i.e., laminar, turbulent, as well as transitional flow regimes. Although the identification procedure is based on steady-state hydraulic network equations, the identified roughness parameters are also key for dynamic models that can be used for model-based controller and observer designs. While a three-cycle network simulation example serves to illustrate the presented problem formulation and solution in an extensive manner, the application on a real-world drinking water network is in focus. In addition, vital aspects, such as topology simplifications of the underlying network and the importance of the generation and selection of independent measurement sets, are addressed. We apply root-finding methods instead of methods based on optimization and thereby show that the pipe roughness identification problem may actually be applicable to identify as well as locate leakages.


I. INTRODUCTION
D RINKING water networks are of utmost importance for society all over the world. However, the share of nonrevenue water, i.e., water that never reaches a registered consumer, accounts for 25%-50% of the total amount of supplied water when put into a global measure (and up to 75% in the emerging markets) as stated by the International Water Association [1]. Consequently, the potential contributions of control system technologies to water distribution systems are manifold and cover, e.g., the coordinated control of distributed actuators (valves and pumps), and the detection of leakages in the water networks [2].
The basis for advanced techniques is usually formed by mathematical modeling of the underlying distribution network. In the literature, two classical modeling approaches are widely used.
1) A description based on the Navier-Stokes equations allows to account for the full transient case [3] with the goal to model effects as detailed as possible. The underlying partial differential equations are, e.g., solved by the method of characteristics or the method of finite differences (see [3]). 2) A second modeling approach uses the relations for the conservation of mass and energy to deduce a description of the main properties under steady-state conditions, as, for example, described in [2]. The latter approach is particularly useful for networks with a large number of pipes and nodes to keep the modeling complexity low. Consequently, the steady-state modeling is well suited for tasks such as pressure control, calibration, consumer demand prediction, and leakage detection in large hydraulic networks and is, thus, also implemented in standard simulation tools for water distribution networks such as EPANET [4], [5] provided by the United States Environmental Protection Agency. In both cases, pipe roughnesses and pressure head losses are essential for the mathematical description of the water network because the water flow distributes along the networks' branches according to the principle of least actions. The corresponding pressure drop due to friction at the pipe's inner surface strongly depends on the flow regime, i.e., on the laminar, turbulent, or transitional (intermediate) regime with the corresponding flow dynamics, see [2], [6], [7] for details as well as [3], [8] for additional aspects of steady-and unsteady-state friction losses, which are usually neglected.
Kaltenbacher et al. [9], [10] formulated dynamic models based on the rigid water column theory. This can be seen as an extension of steady-state relations mentioned in 2) to account for the dynamics in pipe flows, depending on the friction along the pipes (head losses), the pressure heads at the nodes (also at source nodes), and the actual nodal consumption. However, the modeling complexity is still low when compared to the approaches in 1) since ordinary differential equations are utilized for the mathematical description of the fluid flows and pressure heads in the network. Pressure-driven [9] and demand-driven versions [10] exist, depending on how consumer nodes are handled in the mathematical formulation (see [11] for an in-depth description). It is also shown in [10] that the dynamical model's equilibrium is equivalent to the steady-state equations used in EPANET.
For dynamical models as well as for steady-state models, the identification of the friction/roughness coefficients is crucial since it provides the foundation of model-based techniques for the control and observation of nodal pressures and flows [7], [12]- [14], [15]. For example, Kumar and Kumar [16] compared different controllers for water distribution networks and This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Terra et al. [17] presented an H ∞ controller design for such applications.
It is shown in [11] that an observer for the considered water distribution networks specified by the pressure-driven dynamical model exists if the pressure head losses along pipes (and thus also the pipe roughnesses/friction parameters) and the nodal consumption are known. The analysis is based on the concept of the strong zero detectability (see [18], [19], and references therein). In the linear case, strong zero detectability is equivalent to strong * detectability in the sense of Hautus [20]. As a result, the identification of the roughnesses is also necessary for techniques of fault (or leakage) detection and isolation, where the actual behavior of the system under investigation is compared to the nominal one. However, observer design and fault detection techniques are not the focus of this actual work.
To sum up, the identification of pipe roughnesses (and consequently also the friction losses) has a significant role in many domains of control system technologies. Nevertheless, the identification of hydraulic friction parameters is especially challenging because only a very limited number of sensors is present. A vast amount of sensor data is important, e.g., for the localization of leakages that may have a similar impact on the steady-state behavior as increased friction parameters [21]. A positive aspect of the underlying problem is that the roughnesses, i.e., the pressure head loss coefficients, change very slowly over time in real applications. As a result, real-time requirements are not the main issue and the reidentification of the parameters can be triggered if, e.g., new independent measurement sets are available. This is detailed for the real-world drinking water system described in the following.
One way to tackle the identification task is to formulate an optimization problem with a certain (e.g., quadratic) cost to be minimized subject to constraints provided by the model description as, e.g., covered in early works as in [22] and later on in [2]. The underlying steady-state equations for the turbulent region are often not used directly but collapsed into a scalar relation in such approaches. As a result, the components of the gradients of the cost function are not independent any more as it is the case for the original steady-state equations, which makes it harder to find the global minimum. Consequently, solvers, e.g., based on genetic algorithms, are widely applied [23], [24]. The main goal for such algorithms is to reproduce certain signals of interest in the network rather than the true pipe roughnesses. This is often referred to as calibration of water networks in the literature (see [23] and references therein).
In [25], different approaches to directly use the underlying steady-state equations instead of formulating an optimization problem are presented. A comparison of Newton-Raphsontype approaches and successive approximation approaches led to similar results for the considered algorithms, where only the turbulent region is considered. This work was one of the starting points for the present work.
In [26], we primarily focus on the formulation of the problem and deduce concrete conditions that allow to uniquely solve the identification task. Instead of merely calibrating the network, the emphasis is on the identification of the true pipe roughnesses and not measured pressure heads. A set of assumptions [26, Table 1] has thereby been introduced and solution algorithms have been proposed in [26] and extended in [27]. The roughness identification problem formulation in those papers assumes that each pipe flow in the network has to be in the turbulent regime.
This restriction is overcome in this present article using the results of Kaltenbacher et al. [28], where they mathematically formulate the transitional Reynolds region between laminar and turbulent flow in steady state. The existence of an explicit function of the steady-state water flow γ in the transitional region [28], satisfying the boundary and gradient conditions [28, eqs. (i)-(iv)] to the other flow regimes, is exploited in the following. In particular, the satisfaction of gradient conditions at the boundaries enables one to propose a root-finding algorithm for the roughness parameter identification task to autonomously decide whether a pipe flow (in the corresponding measurement set) is in the laminar, turbulent, or even transitional regime. In contrast, EPANET applies a simple cubic interpolation for the transitional region [29] only, resulting in nonsmooth transitions (particularly with respect to the roughnesses) between the flow regions as explained in [28] and does not provide any possibility to identify the parameters.
The authors' previous work about the smooth description of the transitional flow together with the introduced algorithms makes it possible to come up with an advanced identification approach for real water distribution networks, where the actual flow regimes in the considered pipes can be reconstructed. The present work provides the following main contributions.
1) A roughness identification procedure is introduced that allows to account for all three, i.e., laminar, transitional, and turbulent, flow regimes. This is essential for real-world applications in which flow directions are not known a priori in cycled networks. Hence, all flow regimes might be traversed in a smooth way in the course of the solution finding the pipe roughnesses and not measured pressure heads. 2) Two algorithms, based on Newton's method and the tensor approach, are introduced for the underlying identification task. A simulation example of a water network with low complexity is presented to provide insights to the algorithms. 3) The applicability to the real-world drinking water distribution network of Graz-Ragnitz, Austria, is shown. It illustrates the steps to follow for the roughness parameter identification by utilizing the proposed root-finding approach considering all flow regimes. Practical issues related to topology simplifications of the considered water distribution network as well as the importance of the generation and selection of independent measurement sets are discussed in detail. Consequently, the proposed identification procedure forms the basis for advanced model-based control techniques and fault detection and isolation approaches, which consider all three flow regimes.
Notation: Vectors and matrices are highlighted bold and italic and are symbolized using lower and upper case letters, respectively. Symbols 1 x and 0 x represent a matrix or vector filled with ones or zeros, i.e., 1 3 = [1 1 1] T and  [26], [30] It is also utilized to show elementwise exponentiations in a more compact form, e.g.,

II. PRELIMINARIES
This section summarizes preliminary definitions and relations relevant for the identification problem. Details can be found in [26] and [27].

A. Network Hydraulics
One distinguishes between I = {1, 2, . . . , n j } inner nodes/vertices and S = {n j + 1, . . . , n j + n s } source nodes of the network, which is represented by a connected graph [26,Assumption 1]. The intersection I ∩ S = {} is empty and the combination I ∪ S = N = {1, 2, . . . , n j + n s } yields the entire set N of all nodes. Source nodes characterize nodes, where the associated pressure heads [h s ] i = h (i+n j ) for i = 1, 2, . . . , n s (in m) are known, as it is the case for reservoirs and pumps. Vector [h] i = h i∈I (in m) combines all nodal pressure heads at the inner nodes. They can be split into n p (out of n j ) measured heads and not measured heads and its complementary partC h consists of unity vectors e i ∈ Z n j {0,1} with indices i ∈ P = {p 1 , p 2 , . . . , p n p } ⊆ I and i ∈P = {p 1 ,p 2 , . . . ,p n j −n p } ⊆ I, respectively, while respecting P ∩ P = {} and P ∪P = I. As a result, C T h C h +C T hC h = I n j has to hold (see [10,Lemma 1]). In addition, the geographical elevation at each inner node [z] i = z i∈I (in m) has to be considered such that the source (pressure) heads h s are increased by the nodal elevation at the corresponding source nodes S.
The nodal consumption [q] i = q i∈I actually is volumetric flow rates (in m 3 /s) and is a linear combination of P = {1, 2, . . . , n } edge/pipe flows (i.e., volumetric flow rates) symbolized by [x Q ] i = Q i∈P so that where A ∈ Z n j ×n {−1,0,1} denotes the incidence matrix related to the inner nodes I. In contrast, the complete incidence matrix also considers source nodes S with respect toC s . Influent flows to inner nodes I are counted positively and effluent flows to inner nodes I are counted negatively in (2). However, flows effluent of source nodes S are counted positively forC s because they are influent to the I inner nodes by definition. Examples how to build those matrices can be found, e.g., in [26,Secs. 2 and 6].
One also has to account for n c = n −n j linear independent cycles in the network's graph by considering the cycle matrix S ∈ Z (n −n j )×n {−1,0,1} (see [26,Sec. 2]). If [26,Assumption 1] applies, cycle matrix S and the transposed incidence matrix A T are orthogonal, i.e., S A T = 0, where rank(S) = n − n j .

B. Colebrook-White's Turbulent Flow
Turbulent pipe friction is considered via the Darcy-Weisbach formula together with Colebrook-White's implicit friction factor for a Reynolds number of Re ≥ 4000 (see [26] yields It expresses the turbulent flow in each pipe P as a function of the corresponding roughness and the head loss h. Please note that , h, pipes diameter d, cross section A, k = (l/(2dg A 2 )), and length l may differ for each individual pipe. The pipe indices are omitted in (3) to improve readability. The gravitational acceleration is denoted by g ≈ 9.81 m/s 2 . The water density ρ and the dynamic water viscosity η are assumed to be constant although they may vary with temperature.

C. Laminar Flow
A smooth distribution of the velocity profile over the pipe's cross-sectional area can be seen in the laminar flow condition. Applying the relation of Hagen-Poiseuille for the laminar friction factor, one obtains the laminar flow which is an exact mathematical result stemming from the fluid equations of motion [6].

D. Transitional Flow
The authors derived an explicit function [28,Eq. (14) and (43) depending on the pipe's roughness and the head loss h along the pipe, which not only satisfies the boundary conditions to the turbulent and laminar flow but also the corresponding gradients to a sufficient degree of accuracy in the physical relevant range where the relative roughness remains between 0 ≤ /d ≤ 5% according to the Moody chart. Detailed instructions to obtain γ can be found in [28,Sec. 6] and are not included here due to space limitations. Using abbreviation b = 4000(Aη/ρd), the boundary from laminar to transitional flow is located at h = b/w and thus constant over h and , whereas the boundary from the transitional to turbulent flow is described by the curve h = q() by means of function for all h and /d ∈ [0, 5%]. Using abbreviation a = 4000(Aη/ρd), it can be shown that for all h and /d ∈ [0, 5%] with utilizing parameters q 0 , q 1 ∈ R that are determined as indicated in [28,Sec. 6].

E. Combined Flow
Given (3)-(8) while using abbreviation b = 2000(Aη/ρd), the steady-state flow through the pipe can be explicitly expressed by Function (9) is central for the considered roughness identification problem. The mathematical descriptions in the laminar and turbulent flow regime have been validated separately in the literature. In the laminar flow regime, the results are exact representations based on dynamic fluid equations (partial differential equations) as, e.g., detailed in [2] and [3]. The equations by Colebrook-White have been validated for the turbulent region as stated, for example, in [28]. It was shown in [28] for the transitional water flow that the introduced mathematical description improves the existing formulations (e.g., by using cubic functions as in [4]) with respect to the gradient on the boundaries to the laminar and turbulent regime.

F. Network and Sensor Configuration
In order to identify a large number of unknown roughnesses ∈ R n ≥0 and the not-measured pressure heads h (i) for all i ∈ M = {1, 2, . . . , n m }, several sets of linear independent measurements [26, Assumption 5] denoted by M are necessary. The minimal number of measurements is n m,min = n /n p [26, eq. (14)] leading to as many equations as unknowns. Linear independency can be achieved by the variation of the source heads h (i) s and/or the nodal consumptionq (i) along the measurement sets i ∈ M, where the hydraulic network has to be in steady state [26,Assumption 6]. In summary, the following holds.

Assumption 1 (Properties of Network and Measurements):
It is assumed that the following conditions hold. 1(a) The network's graph is connected and has no self-loops. 1(b) The hydraulic friction functions are continuous and strictly monotonically increasing. 1(c) Minor losses are negligible. 1(d) Pipe dimensions, source pressure, and consumption are known. 1(e) Also, a sufficient number of independent measurement sets is available. 1(f) Those measurements are in (hydraulic) steady state. 1(g) Mean-free, white measurement noise with a magnitude smaller than the head loss along each pipe is considered. All assumptions related to the identification problem are more detailed in [26, Table 1], where the same numbering is used. Assumptions 1(a), 1(b), 1(d), and 1(g) are not very restrictive and Assumption 1(c) is usually presumed in calibration methods. However, providing steady-state measurements [Assumption 1(f)] might be challenging in real applications with a large number of pipes. Hence, the proposed identification procedure should be successively applied to different sections of a water distribution network. This is possible by closing valves to isolate a specific section of the network as explained for a real water distribution system in Section VI. The most challenging assumption to be fulfilled is Assumption 1(e), which is related to the number of independent measurement sets. To satisfy Assumption 1(e), additional outflows of the water distribution network can be initiated by using hydrants (which are distributed over the network) to obtain enough independent measurement sets in a steady state (see Section VI for details).

III. PROBLEM STATEMENT CONSIDERING
ALL FLOW REGIMES In this section, the overall parameter identification problem is formulated. We thereby distinguish flow regimes for individual pipes combined in vectors for all j ∈ P and i ∈ M to specify the turbulent flow (10a), the transitional flow (10b), and the proportionality factor (10c) for the laminar flow for all pipes. Entities γ j and f t, j did receive a dedicated index in (10) to highlight that all relations in (10) involve parameters that also vary along P. Also, note that γ (, h (i) ) has to be replaced with sign(h (i) ) γ (, |h (i) |) for negative h (i) (see [28,Remark 2]). One has to ensure that each pipeflow is unique among each measurement set. Hence, the corresponding sets of pipe flow indices along the flow regimes. To be able to formulate the problem in a unified and compact fashion for all flow regimes, flow separation matrices . .
and flow separation vectors are introduced. They are constructed from unity vectors e j ∈ Z n {0,1} ∀ j ∈ P using the number of pipes of measurement set i in the turbulent n (i) t , transitional n (i) γ , and laminar n (i) L region. The number of pipes in the different flow regimes may vary between all i ∈ M. Consequently, it follows that where I n is the identity matrix with dimension n (i) t +n (i) γ +n (i) L = n in which n t , n γ , n L may also vary in each measurement set. Please do not confuse the number of pipes n with the number of laminar flows n L in the i th measurement set. Although the flow separation matrices actually depend on the variables and h (i) N , they are only comprised of zeros and ones, i.e., their derivative always yields the zero matrix.
Finally, applying (10) and the flow separation matrices (12), one receives the equation set at nodes with no sensors in measurement sets i ∈ M that are the unknowns to be found. Please note that the flow (9) is a nonconvex function by nature, whenever all three flow regimes are considered. The variation of initial conditions primarily intends to relax Assumptions 1(e)-1(g), while it also helps to handle the nonconvex nature of function (9) as explained in detail next.
One does not have to take the flow regime case separation into account explicitly, and the consideration of matrices is sufficient. The consideration of all flow regimes facilitates the transition between laminar and turbulent flow regimes in the solving of problem (14). This allows individual pipe flows in P, in principle, to never become turbulent (or even transitional) in any of the i ∈ M measurement sets. Hence, the following assumption is imposed in addition to Assumption 1.
Assumption 2: All individual pipe flows [x (i) Q ] j must at least once be in the turbulent regime in any of the i ∈ M measurement sets, i.e., In the following, numerical algorithms are introduced to solve (14) while incorporating the different flow regimes.

IV. SOLUTION ALGORITHMS
This section briefly introduces the algorithms presented in [26] and [27] and their extensions to solve problem (14). However, an in-depth description of these algorithms can be found in [27,Sec. 7.1] and [26,Secs. 5 and 6]. It is important to emphasize that the authors use numerical root-finding methods on the basis of Newton-Raphson algorithms rather than optimization-based methods that are very common in this research area [13], [31], [32]. This is possible due to the actual formulation providing, at least, as many unknowns as equations.
Generally, two different types of algorithms have to be distinguished. Type (I) algorithms handle the classical Newtonian iterative solution of a nonlinear yet smooth equation set with step length variation as proposed in [26,Algorithm 1]. Since these algorithms of type (I) still require to start in the vicinity of the real root, denoted by , for convergence, a type (II) algorithm then launches the type (I) algorithm (e.g., [26,Algorithm 1]) several times with different initial values. This also accounts for the nonconvexity of the flow (9), which is not a specific feature of the proposed approach but a direct consequence of the considered physical phenomena (see [28]).
The intermediate best result in terms of the residual v(x) = || f (x)|| L 1 along the different initial values is denoted by x + . In this context, function f (x) follows from (2) and (14) such that: . . .
The combination of these two different kinds of algorithms is necessary for solving as exemplified in the authors' previous publications on the full turbulent case [26], [27] to relax the necessity to start close to real root of . Referring to [27], the extension of second-order derivatives in the determination of a search direction is studied. Using the iteration index k to denote steps in the iterative solution finding, a Taylor series is truncated after the linear term in the Newton-Raphson Algorithm 1 Type (I) Algorithm: Modified Line Search With Step Length Variation (Newton or Tensor Approach) algorithm to obtain a search direction In contrast, the tensor equation as presented in [27, eq. (12)] for the full turbulent case also considers the second derivatives in the Taylor series so that where H([ f (x k )] i ) denotes the Hessian of the i th component of function f (x) (15). Although performance improvements could be achieved with respect to the residual v(x) = || f (x)|| L 1 , this involves considerable effort in terms of the implementation and computational demand. This is particularly true since (16) is also solved iteratively with the MATLAB built-in fsolve(.) function in this work. Details of this tensor method are provided in [27], whereas also, Algorithm 2 Type (II) Algorithm: Variation of Initial Values a customized method to vary the step length μ in x k+1 = x k + μx k is applied.
Our previous results from [26] and [27] can be extended to the actual problem formulation, i.e., the solution of (14), by replacing partial turbulent flow derivatives p (i) X (according to the definition in [27] as vector) with respect to X ∈ {, h, 2 , h, h 2 } with where X ∈ {, h, h 2 , 2 }. In addition, the partial derivatives γ (, h (i) ) and γ h (, h (i) ) have to be replaced by for all i ∈ M to also consider negative head losses, i.e., h (i) < 0. The gradient of γ with respect to depends on the sign of h (i) as can be seen in (9), which is not the case for γ h . Consequently, one can consider all three flow regimes in the proposed identification approach although the actual regime (per pipe) is not known a priori. Explicit relations for the first and second derivatives needed for the Newton and tensor method are stated in the Appendix. The algorithms used in this article follow from the algorithms in [26] and [27], where (17) and (18)   Hydraulic network with one reservoir, eight pipes, three cycles, and pressure sensors at nodes k = 2, 3, 4 (red), see [26], [27]. The presented algorithms are first applied to a three-cycle network for illustration purposes and then to the real-world water distribution system of Graz-Ragnitz in the following.
V. THREE-CYCLE NETWORK EXAMPLE A hydraulic network with low complexity is used to show the main steps to apply the algorithms described in Section IV for solving equation set (14).

A. Configuration
The network shown in Fig. 1 consists of a reservoir, i.e., n s = 1, that provides the source head h s , n = 8 pipes, n j = 5 nodes (n p = 3 pressure nodes), and n c = 3 cycles. For analysis purposes, two measurement sets additional to the minimal required n m,min = n /n p = 3 linear independent measurement sets M were generated to investigate performance improvements by accounting for more than n m,min measurement sets. In addition, (mean-free white) measurement noise is added to the pressure sensors and the consumption data q 2 , q 3 , and q 4 (e.g., records of fireflows). Fig. 2 shows the consumer pattern that is chosen for the simulation example. Steady-state simulation results of the dynamic PD u model, as presented in detail in [11], are taken as artificial measurement sets. This dynamic model applies a pressure-driven consumption approach, where the consumer demand is determined by the pressure at respective nodes instead of a priori values for the consumption (cf. Fig. 3). The set of "measured" values in the steady-state time frames  shown in Figs. 2 and 3 is utilized for the subsequent pipe roughness identification by averaging with the aim to mitigate noise effects, which should also be done when using realworld data. Steady-state pipe flow 5 in measurement set 3, i.e., Q (3) 5 , as well as steady-state pipe flow 8 in measurement set 4, i.e., Q (4) 8 , are laminar as can be seen in the zoom window in Fig. 4. All other flows are in the turbulent regime. See Figs. 2 and 3 for the time windows that correspond to the different measurement sets i .
Two particular reasons, which render this current configuration interesting, shall be highlighted. First, it will be clarified if the derived transitional flow γ (5) is eligible to provide the transition from turbulent to laminar flow in solving (14). Second, insight will be given to which extent mean-free measurement noise in y h and q interferes with the roughness identification.

B. Results and Discussion
Tables I-VI represent the results of the roughness identification, i.e., the solution of equation set (14) using either the Newton approach (see [26,Algorithm 1] in combination with [26, Algorithm 2]) or the tensor approach (see [27,Algorithm 3] [26], [27] for a detailed description of the algorithms). The computational duration stated in the tables result for an implementation of the proposed algorithms in MATLAB 2 on a low-power notebook from 2015 with an Intel dual-core i5 processor and 8-GB RAM. The focus is not on real-time requirements since the roughness parameters change very slowly over time. Thus, the use of a standard notebook is sufficient in the present case to point out the main properties of the proposed approach.
The intermediate best result in terms of the residual v(x + ) is highlighted in x T h N ,0 are the estimated mean values for the pressure heads, where the knowledge of the source, the pressure sensors, as well as the network structure ( Fig. 1) are exploited.
In Tables I-VI, different numbers of measurement sets are used to show their impact on the identification of the pipe roughnesses and not measured head losses h N,1 and h N,5 at nodes 1 and 5 (see Fig. 1) in different measurement sets, respectively.
One can see that the results are noticeably consistent, cf. values of the identified parameters in the columns with the lowest values for v(x + ), which are marked in green. The residual of the results is approximately two orders of magnitude smaller than the one of the real root v(x * ) (especially in Tables I and II) as the exact solution of without measurement noise. Concerning the regime case separation, the tensor as well as Newton method managed to successfully identify Q (3) 5 (in case of n m = 3) as well as Q (3) 5 and Q (4) 8 (in case of n m = {4, 5}) as laminar and all other steady-state flows as turbulent by evaluating (9) for x + or by checking the flow separation vectors (13). This underpins the applicability of the transitional flow characteristics γ (5) for the complete roughness identification problem (14). Note that the consideration of laminar flows improves the reconstructability of original roughnesses as the known linear (laminar) terms indeed facilitate solving . This can be seen when comparing tables considering n m = 3 with ones considering n m = 4 measurement sets as the results in the former ones are much closer to the original root. In addition, the inclusion of measurement set i = 5 in the problem setup aggravates the solution finding as additionally uncertainties come into play.
Comparing tensor with Newton results for all n m , no significant improvements could be achieved when considering the second-order derivatives concerning the tensor method. Actually, a slightly smaller residual was found by the Newton method for n m = 5 as can be seen by comparing Table V with Table VI. This has three reasons. First, the second derivatives of γ (see the Appendix) are neglected in the solution finding of the search direction via the tensor equation (16). Second, when solving (16) iteratively with the MATLAB built-in fsolve(.) function, it occasionally happens that a search direction is taken, which has a comparably high residual measurement set i in solving iteration k of . In this context, one can observe that the tensor method indeed takes more iterations on average to converge or abort. Third, as long as measurement sets are sufficiently independent, the application of the Newton method might presumably suffice in combination with [27,Algorithm 4] to obtain a solution that features an even smaller residual v(x + ) than the original root v(x * ). Tables I-VI, the comparably high deviation in the solution of the not-measured pressure head at node 1, i.e., h (i) N,1 ∀i ∈ M, is also noteworthy. The deviation in h (i) N,1 is inherently connected with the deviation in the roughness of the first pipe. Although the estimated roughness of the first pipe, i.e., 1 , is close to the original one * 1 in all the above tables (particularly for the case n m = 4), the error in h (i) N,1 ∀i ∈ M is visible. The major reason for this can again be found in the nonlinear friction relation, where comparably small differences in roughness lead to a substantial deviation in the head loss and thus in the nodal head due to the high flow in pipe 1.

In view of the results in
In sum, the collection of suitable measurement sets is key for the reconstruction of the pipes' roughness. The best indicator for the quality of the measurement sets is given by the differences among y h + C h z (cf. Fig. 3) and the differences in the consumption q (see Fig. 2). However, when considering minor losses (violating the assumption about the minor losses), i.e., nonzero minor loss parameters [9], and measurement noise that features a nonzero mean, solving becomes substantially more challenging in a sense that one commonly receives estimated roughnesses outside their physical range (above 5% of the pipe's diameter). These solutions, which lie outside the considered physical range, sometimes indeed have comparably low residuals, due to the aggressive search by [27,Algorithm 4]. Nevertheless, in the opinion of the authors, solving problem (14) offers the best chance to reconstruct individual roughnesses per pipe. However, it is very common to group pipes' roughnesses for the solution finding, meaning to assign a single roughness value to multiple pipes in order to compensate for additional unknowns.
The three-cycle network example gives valuable insights about the applicability of the presented roughness identification scheme as follows.
1) Adverse effects due to measurement noise have been studied. 2) The consideration of laminar flows turns out to facilitate the solving of .
3) The functionality of γ (5) to allow a transition from turbulent to laminar flow (or vice versa) in the solution finding of has been demonstrated. 4) The analysis of considering more than n m,min measurement sets reveals that it does not always facilitate solving , but a careful selection of measurement sets may be necessary. In Section VI, the proposed approaches are applied to realworld measurement data.

VI. ROUGHNESS IDENTIFICATION ON A REAL NETWORK
The application of the proposed roughness identification approach is now shown for the real-world drinking water distribution network of Graz-Ragnitz, Austria [33]. The measurements used for this analysis were conducted prior to the involvement of the authors with the help of the local water utility in the course of a project with the Institute of Urban Water Management, Graz University of Technology.

A. Topology and Its Simplification
During the time frame where measurements have been recorded, the network according to Fig. 5 was isolated by valves from other parts connecting itself to a larger distribution network. The topology of this isolated network (part) as presented in Fig. 5 consists of the elements stated in Table VII. Hence, given the number of n = 627 pipes and n p = 13 pressure sensors, one would need to produce at least n m,min = n /n p = 48 independent measurement sets with only n q = 4 fireflows. It is reasonable to assume that there is no way to produce so many independent measurement sets, especially on a high fireflow level necessary to produce enough head loss. The only chance is to simplify the topology to the highest extent possible and thereby reduce the number of pipes and nodes in order to result in a sufficiently small n m,min = n /n p . This is explained in detail in the following.  6. Result of combining adjoining pipes with the same diameter. Nodes with more than two connections, which are used to simplify the topology, are highlighted in green.

1) Dead Ends:
The first step of the topology simplification involves the removal of dead ends, which also include nodes that represent real consumers (not hydrants) in the original topology. However, as it has to be assumed that they do not retrieve water during the identification measurements (at night), they have to be removed as no corresponding roughness values can be identified.
2) Combine Adjoining Pipes With the Same Diameter: The next step involves the detection and then combination of adjoining pipes with the same diameter by removing the node in between, i.e., the new combined pipe's length has to be adjusted. Pressure and hydrant nodes have to be excluded from this removal. In Fig. 6, one can see the result of the iterative rejoining of pipes with the same diameter. In this process, the x-and y-coordinates of the original nodes have been preserved, which causes, for instance, the small cycle at "HG3933" seen in Fig. 5 to seemingly disappear in Fig. 6. This cycle is still present in the graph of Fig. 6 at this point of simplification, and however, one would have to zoomed-in view substantially.

3) Remapping Pressure Nodes and Hydrant Nodes:
This section lists four measures to simplify the graph's topology.
1) As pressure nodes with no hydrants on the same link do not have any head loss along this link, it is feasible to remap the pressure nodes to the next respective junction, i.e., a node with more than two connections. In the course of this process, one has to carefully consider the difference in elevation between the pressure node and the final junction node as the height difference has to be considered in the corresponding pressure readings. Thus, the pressure nodes, i.e., the nodes with pressure sensors, are moved to the green junctions in Fig. 6. 2) It turned out that two deployed pressure sensors, namely, "HG3445" and "HG4215," are redundant for all measurements sets, i.e., "HG3420" and "HG3445" or "HG4162" and "HG4215" have the same pressure readings as there is no flow in the links of these two sensor pairs. Hence, the corresponding nodes are removed from the graph. This then enables to move pressure nodes "HG3420" and "HG4215" to the next higher junction with more than two connections. 3) In analogy to 2), the (small) cycle (seen in Fig. 6) enclosing pressure node "HG3933" does not have any hydrant to cause water to flow into this network part.
As a result, there should be no head loss in this entire cycle, leading to a uniform head distribution up to the next junction. As side note, this is of course only feasible if no background consumption/leakage occurs there. Effectively, the pressure readings of "HG3933" are moved to the next higher junction also accounting for differences in nodal elevation. 4) Similar to 1), as hydrants are also located on dead-endlike links, they are moved to the next higher junctions that are colored in green. This is only feasible due to the assumption to take the measurements in steady state only.

4) Final Graph:
The final graph's topology used for the roughness identification is shown in Fig. 7. Mind that the identifiers (IDs), e.g., "HG3537," are actually used to denote individual nodes of the original graph in Fig. 5, meaning that the highlighted pressure and hydrant nodes of the final graph in Fig. 7 do have different names, i.e., IDs, internally. However, in order to not confuse the reader, the original IDs have been preserved to support identifiability with the corresponding sensors. Effectively, one can see that the graph has been simplified substantially in comparison to the original one in Fig. 5, as shown in Table VIII. Given the simplified graph, at least n /n p = 7 independent measurement sets are needed to reconstruct suitable roughness values for n = 74 (combined) pipes in the drinking water distribution network.  (FIG. 7)   Fig. 8. Head readings of all (relevant) pressure sensors according to the graph in Fig. 7. The steady-state time frames used for averaging are highlighted with black dashed lines at the starting and end points. For instance, the fourth time frame can be found between t ∈ [5080, 5140] s. This is a significant improvement compared to n m,min = 48 when considering the original topology in Fig. 5.

B. Sensor Readings and First Assessments
Generally, hydrants are equipped with flow sensors in order to record the fireflow, treated as "consumption" q, while the pressure sensors record the nodal pressure heads y h during the minimum night flow where the real, unknown consumption (also including losses) is lowest. In Fig. 9, one can see the individual fireflow measurements at the n q = 4 hydrants, whereas Fig. 10 compares the sum of all fireflows of the four hydrants with the total network's inflow that was measured separately at the link where the source connects to the rest of the network (cf. topology in Fig. 7). Counting the peaks in Fig. 9, one can see that there are potentially 15 measurement sets available. However, it will turn out that approximately half of these sets are not suitable for identification. In Fig. 8, one can see all relevant head readings of the pressure sensors with respect to the IDs found in the graph of Fig. 7. Considering Figs. 8 and 9, observations and the applied approaches are summarized in reference to the assumptions made in this work: 1) Assumptions 1(a), 1(b), and 1(d) are fulfilled for the current application. 2) It is, in principle, possible to cause sufficiently large head loss by fireflows to satisfy Assumption 1(g) and thereby provide sufficiently independent sets, see Assumption 1(e). 3) Referring to head readings, Assumption 1(f) to only measure in steady state does seem legitimate yet in very small time frames only. The numerous occurrences of oscillating heads and peaks indicate frequent transient events.   Fig. 10 shows the presumably major issue of these measurements as the unknown background consumption, i.e., the difference between blue and red curve, makes up to approximately 40% (between t ≈ [2900, 2960] s) of the total consumption in the corresponding measurement set. The only possibility to circumvent this is to avoid high-erroneous measurement sets, with respect to the background consumption, from (14) and to take only those were the error is smallest, e.g., the first time frame at t ≈ [2380, 2450] s, as shown in Figs. 8 and 10. 6) Three time frames of comparably small fireflows around t ≈ 6000, 7000, 12 000 s at a level of 1 l/s in Fig. 9 actually do not cause any notable head loss as shown in Fig. 8 and are therefore completely unsuitable for the consideration in . 7) The chosen n m = 8 time frames, i.e., measurement sets, for the averaging of measurements had to be chosen considerably small to avoid any possibly transient events, see time frames between the dashed lines in Fig. 8. 8) The minimum night flow, which gives a good estimate of the total water loss, can be identified in the range of 1.2 l/s according to Fig. 10. These background losses are potentially problematic for any identification scheme that is based on the satisfaction of nodal Kirchhoff equations, see (14). However, this is a general problem with all sorts of roughness identification algorithms as these background losses at the time of the identification measurements ultimately lead to mistakenly increased friction parameters. One can hardly detect and localize leakages in an uncalibrated hydraulic network, meaning that all leakages prior to the identification (identification of friction parameters) are likely to remain invisible to the observer. This certainly depends on the size of the leakage nonetheless. If it exceeds a critical size, the detection and localization become feasible at some point.

C. Initial Values
A crucial point for the considered problem is the selection of the measurement sets and the stochastic variation of the initial values for the launch of the roughness identification concerning (14) as described in Section IV. This is due to the fact (see also [28]) that the considered flow (over all three regimes) is characterized by a nonconvex relation (9) due to the underlying physical principles.
In particular, the selection of suitable initial values for the nodal pressure head values at nodes with no sensors concerning h (i) for all i ∈ M = {1, 2, . . . , 8} is challenging and leaves a lot of potential variables to be adjusted. However, the stochastic variation of initial values as explained in Section IV helps to relax the requirement of already starting with values close to the ones to find.
Similar to the three-cycle network in Section VI-B, the initial values for h N are chosen by means of averaging heads between pressure nodes. Details about the choice of initial conditions can be found in [11]. As in the previous example, the initial roughness values for launching the identification are selected as 1% of the corresponding pipe's diameter. Mind that according to the proposed algorithm, these initial roughnesses are already varied before launching Newton or tensor algorithms.

D. Results and Discussion
The obtained roughness results are presented when applying a few iterations of [27,Algorithm 4], which itself launches [27,Algorithm 3]. This process is then stopped after a few iterations because a clear trend can be observed, which will be subject to discussion.
To begin with, the equation set (14) to be solved features 576 (nodal) equations and 562 variables from which only 74 elements are the roughness values of particular interest. The residual could be reduced from a total error of v(x 0 ) = 33.2020 m 3 /s to v(x + ) = 8.4101 m 3 /s at the 576 nodes considered in solving. Looking at the average error on each node over all eight measurement sets, a presumably more distinctive value for analysis, the error could be reduced from 57.6 to 14.6 l/s. Considering that the fireflows never exceed the 16 l/s mark (see Fig. 10), the error is in a range, which raises questions about the quality of the data. To put the residual in a different perspective, 20.7 nodes exceed the 5-l/s error and 29.1 nodes exceed the 1-l/s error from the total of n j = 74 on average over all eight measurement sets. Fig. 11 shows the results for the obtained relative roughness values directly on the corresponding pipes and highlights 32 pipes that have a relative roughness greater or equal to Fig. 11. Obtained relative roughness values in percent along the pipes of the simplified network's graph highlighting in red those which are greater or equal to /d = 6%. Nodes with an error greater or equal than 10 l/s in one of the eight measurement sets are highlighted in red. Nonpressure nodes, whose head values are not inside their considered range in the respective measurement set (cf. Fig. 8), are mentioned in the legend. /d = 6% (red lines). Fig. 11 also highlights nodes with an error greater or equal than 10 l/s in one of the eight measurement sets (red dots). One can clearly see that there is a correlation between high-erroneous nodes and unrealistically high roughness values. Please note that we aim to find the roughness parameters for all pipes in the considered water distribution system. However, the real roughness values, which would allow a comparison to the identified values, cannot be found for real applications with a reasonable effort. Under the assumption that leakages can be present at each of the n j inner nodes, one would need n j pressure sensors and n flow sensors. In the other extreme case, where no leakages are present, one would need n j pressure sensors and n c flow sensors. Neither case is realistic in real water distribution networks.
Effectively, two major issues make it difficult to solve and obtain a smaller residual. First, the occurrence of the comparably high background consumption/leakage seen in Fig. 10 causes inherently higher than actual roughnesses. Due to the nonlinear friction relations, an unrealistically high roughness might be necessary to compensate for the unaccounted flow at respective nodes. Second, it turned out that apart from the background leakage, a valve in the lower left corner in Fig. 11 was closed partially during the measurement. Interestingly, however, one can actually see that the exact area of the partially closed valve was identified by the presented roughness identification scheme. This demonstrates that the developed approach may also be suitable to identify network parts inconsistent with the measurement data, e.g., due to leakages. The area close to pressure sensor "HG3933" happens to be red-colored because of the significant background leakage/consumption in the removed cycle. Nevertheless, this cycle had to be removed as no hydrant was opened in this cycle during the measurements.
The results for h N remain in their physical range, except at nodes that are highlighted in Fig. 11 in the respective measurement set (see legend). This range violation is very minor for the yellow and orange colored nodes, for instance, the orange colored nonpressure node results in h N + z = 497.9588 m, which exceeds the identified upper limit h N +z = 497.9061 m by a mere 5.27 cm. Also, regarding the orange and yellow nonpressure nodes of Fig. 11, heads h N tend to leave their intended range at exactly those nodes adjacent to pipes with very high roughness values, which indicates a correlation between them. The red-colored node of Fig. 11 has h N + z = 476.3744 m but should actually be above h N + z = 480.7406 m. Knowing that there was a partially closed valve as shown in Fig. 11, the developed identification algorithm is able to confidently find discrepancies in the measurement data.
When taking the results as initial value for the next iteration of, e.g., [27,Algorithm 3] and relax the requirement for the not-measured pressure heads h N to remain inside their selected limits x h N ∈ [h N , h N ], one indeed manages to obtain smaller residuals. However, the roughnesses exceedingly leave their physical range and therefore do not represent results, which can be seriously applied. The effect of nonzero minor losses thereby violating the assumption about negligible minor losses may be of particular interest in this context.

1) Flow Regimes:
For the considered water distribution network, there are n n m = 592 pipe flows to be categorized among the three flow regimes. Thereby, only 238 pipe flows have been identified to be in the turbulent regime, 219 in the transitional regime, and hence 135 in the laminar regime. This result underlines the importance to distinguish among the three flow regimes in the solving process of . Violating Assumption 2 [1(a)], pipes have flows, which never reach the turbulent regime in any of the eight measurement sets. According to the results of the roughness identification, there is only a single pipe whose flows never become turbulent nor transitional in any of the measurement sets. As shown in Fig. 11, the obtained roughness of this pipe then also exceeds the physically reasonable values, as the variation of this roughness value has no effect on the residual of in the solving process. Looking at Fig. 7 for instance, it is clear why this pipe can neither have turbulent nor transitional flow as there cannot be any flow between "HG4118" (a hydrant) and "HG4150" (a pressure sensor) according to the hydrant configuration.
This is a remarkable result nonetheless and proves the applicability of γ (5) with respect to to properly distinguish among flow regimes.
In analogy, the pipe whose flow never reaches the turbulent regime in Fig. 11 has a very high relative roughness of 43.3% according to Fig. 11. Knowing that the smoothness and validity of γ are only preserved in the range of /d ∈ [0, 5]% of the pipe's diameter (see [28]), the root-finding algorithm for particularly struggles to handle roughness values greater 5% when being in the transitional regime. Once in the transitional regime and outside the physical limitation, the algorithm cannot again converge to the turbulent regime.
2) Linear Independency: Recall that n m ≥ n m,min = n /n p is only a necessary condition to satisfy the assumption about the number of linear independent measurements. Apart from that, the best indication about the linear independency is given by the rank of the Jacobian, which is a function of the solution in the iteration step k, i.e., J (x k ). This Jacobian has size J ∈ R 576×562 for the eight measurement sets. One obtains rank( J (x 0 )) = 545 in x 0 , and in the solution corresponding to the roughness values in Fig. 11, it reaches rank( J (x + )) = 549, i.e., its deficiency concerns at least 13 variables among n + n m (n j − n p ) = 562. This has two reasons. First, as discussed previously, there are pipes whose flows are never in the turbulent regime making it impossible to find the corresponding roughnesses. Second, although there is, in principle, enough head loss along some pipes according to the final graph in Fig. 7, this presumably does not apply to all n = 74 pipes when inspecting Fig. 8. Mind the assumption about the measurement noise in this context. The only possibility to improve on this is to consider further and/or different measurement sets with other fireflow configurations, whereas additional fireflows provided by other than the n q hydrants may be particularly valuable in this regard.

VII. CONCLUSION
The basic applicability of the developed roughness identification on a real drinking water distribution network was demonstrated, considering all three flow regimes. The difficulties and uncertainties when using real-world data were discussed, whereas the problems in the solution finding and results could be traced back to inconsistencies in the data and, for instance, background consumption. The legitimacy of Assumptions 1 and 2 have been studied and turn out to be attainable to a large extent. Nevertheless, the validity of the neglection of minor losses concerning Assumption 1(c) has yet to be verified in greater detail. In short, the authors believe the developed identification scheme not only to be eligible for roughness identification but for the detection and localization of leakages in water distribution systems, which is subject to further research.

A. First Derivatives
The derivatives of turbulent flow are given in [26]. In addition, the partial derivatives of the transitional water flow γ (5) with respect to the roughness , i.e., [γ ] j , and its derivative with respect to the pressure head loss, i.e., [γ h ] j , are shown in for all j ∈ G (i) ⊆ P ∧ i ∈ M are considered, where A = [a 1 a 2 · · · a n ] in this context. As first-order γ -derivatives (19) and (20) are already complex expressions, the secondorder γ -derivatives, i.e., γ (i) X , j for X ∈ { 2 , h, h 2 }, have been neglected in this work. Details can be found in [27], where p (in the sense of [27], i.e., p only applies for the turbulent regime) has to be replaced by γ.