Blockage detection in networks: the area reconstruction method

In this note we present a reconstructive algorithm for solving the cross-sectional pipe area from boundary measurements in a tree network with one inaccessbile end. This is equivalent to reconstructing the first order perturbation to a wave equation on a quantum graph from boundary measurements at all network ends except one. The method presented here is based on a time reversal boundary control method originally presented by Sondhi and Gopinath for one dimensional problems and later by Oksanen to higher dimensional manifolds. The algorithm is local, so is applicable to complicated networks if we are interested only in a part isomorphic to a tree. Moreover the numerical implementation requires only one matrix inversion or least squares minimization per discretization point in the physical network. We present a theoretical solution existence proof, a step-by-step algorithm, and a numerical implementation applied to two numerical experiments.


Introduction
We present a reconstruction algorithm and its numerical implementation for solving for the pipe cross-sectional area in a water supply network from boundary measurements. Mathematically this is modelled by the frictionless waterhammer equations on a quantum graph, with Kirchhoff's law and the law of continuity on junctions. The waterhammer equations on a segment P = (0, ) are given by [1,2] ∂ t H(t, x) = − a 2 (x) gA(x) ∂ x Q(t, x), t ∈ R, x ∈ P, (1.1) where • H is the hydraulic pressure (piezometric head) inside the pipe P , with dimensions of length, • Q is the pipe discharge or flow rate in the direction of increasing x.
Its dimension is length 3 /time, • a is the wave speed in length/time, • g is the graviational acceleration in length/time 2 , • A is the pipe's internal cross-sectional area in length 2 . We model the network by a set of segments whose ends have been joined together. The segments model pipes. In this model we first fix a positive direction of flow on each pipe P j , and set the coordinates (0, j ) on it. Then we impose Equations (1.1) and (1.2) on the segments. On the vertices, which model junctions, we require that H(t, x) has a unique limit no matter which direction the point x tends to the vertex. Moreover on any vertex V we require that where ν j ∈ {+1, −1} gives the direction of the internal normal vector in coordinates x to pipe P j at V , and Q j is the boundary flow of pipe P j at vertex V . In other words ν j Q j is the flow into the pipe P j at the vertex V . Hence the sum of the flows into the pipes at each vertex must be zero. Therefore there are no sinks or sources at the junctions, and also none in the network in general, except for the ones at the boundary of the network that are used to create our input flows for the measurements.
The inverse problem we are set to solve is the following one. We assume that the wave speed a(x) = a is constant. If not, see Section 6. Given a tree network with N + 1 ends x 0 , x 1 , . . . , x N , we can set the flow and measure the pressure on these points except for x = x 0 . Using this information, can we deduce A(x) for x inside the network?
The problem of finding pipe areas arise in systems such as water supply networks and pressurized sewers, due to the formation of blockages during their lifetime [3,4]. In water supply pipes, blockages increase energy consumption and increase the potential of water contamination. Blockages in sewer pipes increase the risk of overflows in waste-water collection systems and impose a risk to public health and environment. Detection of these anomalies improves the effectiveness of pipe replacement and maintenance, and hence improves environmental health.
An analogous and important problem is fault detection in electric cables and feed networks. This problem can be solved in the same way as blockage detection because of the analogy between waterhammer equations and TelegrapherâĂŹs equations [5].
Mathematically the simplest network is a segment joining x 1 to x 0 . In this setting, the problem formulated as above was solved by [6]. Various other algorithms were developed for solving the one-segment problem both before and after. See [7] for a review. After a while some mathematicians started focusing on the network situation and others into higher dimensional manifolds. Others continued on tree networks, which are also the setting of this manuscript. For tree networks Belishev's boundary control method [8,9] showed that a tree network can be completely reconstructed after having access to all boundary points. These results were improved more recently by Avdonin and Kurasov [10]. A unified approach that uses Carleman estimates and is applicable to various equations on tree networks was introduced by Baudouin and Yamamoto [11]. Also, [12] implies that it is possible to solve the problem in the presence of various different boundary conditions when doing the measurements. A network having loops presents various challenges to solving for the cross-sectional area from boundary measurements [13]. In [10] the authors furthermore show that all but one boundary vertex is enough information, and also that if the network topology is known a-priori, then the simpler backscattering data is enough to solve the inverse problem. This latter data is easier to measure: the pressure needs to be measured only at the same network end from which the flow pulse is sent.
Our paper has two goals: 1) to solve the inverse problem by a reconstruction algorithm, and 2) to have a simple and concrete algorithm that is easy to implement on a computer. In other words we focus on the reconstruction aspect of the problem, given the full necessary data, and show that it is possible to build an efficient numerical algorithm for reconstruction. From the point of view of implementing the reconstruction, the papers cited above are of various levels of difficulty, and we admit that this is partly a question of experience and opinion. We view that apart from [6], all of them focus much on the important theoretical foundations, but lack in the clarity of algorithmic presentation for non-experts. Therefore one of the major points of this text is a clearly defined algorithm whose implementation does not require understanding its theoretical foundations. See Section 4.
Our method is based on the ideas of [6,14] and is a continuation of our investigation of the single pipe case [3]. In both of them the main point is to use boundary control to produce a wave that would induce a constant pressure in a domain of influence at a certain given time. A simple integration by parts gives then the total volume covered by that domain of influence. Time reversal is the main tool which allows us to build that wave without knowing a-priori what is inside the pipe network.
A numerical implemetation and proof of a working regularization scheme for [14] in one space dimension was shown in [15]. Earlier, we studied [6] and tested it numerically in the context of blockage detection in water supply pipes [3] and also experimentally. We also showed that the method can be further extended to detect other types of faults such as leaks [16]. The one space-dimensional inversion algorithm can be implemented more simply than in [14,15], even in the case of networks. This is the main contribution of this paper. Furthermore we provide a step-by-step numerical implementation and present numerical experiments. Our reconstruction algorithm is local and time-optimal. In other words if we are interested in only part of the network, we can do measurements on only part of the boundary. Furthermore the algorithm uses time-measurements just long enough to recover the part of interest. Intuitively, to reconstruct the area A(x) at a location x, the wave must reach the location x and reflects back to the measurement location. Any measurement done on a shorter time-interval will not be enough to recover it.

Governing Equations
Networks. Denote by G a tree network. Let J be the set of internal junction points and P the set of pipes, or segments. The boundary ∂G consists of all ends of pipes that are not junctions of two or more pipes. Each of them belongs to a single pipe unlike junction points which belong to three or more. There are no junction points that are connected to exactly two pipes: these two pipes are considered as one, and the point between is just an ordinary internal point.
Each pipe P ∈ P is modelled by a segment (0, ) where is the length of the pipe. This defines a direction of positive flow on the pipe, namely from x = 0 towards x = . The pipes are connected to each other by junction points. Write (x, P ) ∼ v if v ∈ J, x is the beginning (x = 0) or endpoint (x = ) of pipe P ∈ P, and the latter is connected to v at the beginning (x = 0) or end (x = ), respectively. Equations. Inside pipe P the perturbed pressure head H and cross-sectional discharge Q satisfy where a, A, g denote the wave speed, cross-sectional area and gravitational acceleration. The sign of Q is chosen so that Q > 0 means the flow goes from 0 to . In Example 2.1 the positive flow goes from A to D, from B to D and from D to C. The wave speed is assumed constant. The pressure is a scalar: if v ∈ J connects two or more pipes (x j , P j ) ∼ v, j = 1, 2, . . ., then H has a unique value at v The flow satisfies mass conservation, i.e. a condition analogous to Kirchhoff's law. The total flow into a junction must be equal to the total flow out of the junction at any time. To state this as an equation we define the internal normal vector ν for any pipe end. If the pipe P has coordinate representation (0, ) then ν(0) = +1 and ν( ) = −1. Recall that Q > 0 mean a positive flow from the direction of 0 to . This means that ν(x)Q(t, x) is the flow into the pipe at point x ∈ {0, }. If it is positive then there is a net flow of water into P through x. If it is negative then there is a net flow out of P through x. Mass conservation is then written as Initial conditions. The model assumes unperturbed initial conditions for any pipe P ∈ P.
Direct problem. In this section we define the behaviour of the pressure H and pipe cross-sectional discharge Q in the network given a boundary flow and network structure. This is the direct problem. Recall that we have one inaccessible end of the network, x 0 , whose boundary condition must be an inactive one, i.e. must not create waves when there are no incident waves. To make the problem mathematically well defined we choose for example Equation (2.7). Other options would work too, for example involving derivatives of H or Q, but it is questionable how physically realistic such an arbitrary boundary condition is. The main point is that any inactive condition at the inaccessible end is allowed without risking the reconstruction. This is because the theoretical waves used in the calculations of our reconstruction algorithm never have to reach this final vertex, and so the boundary condition there never has a chance to modify reflected waves.
for some given functions or constants A(t), B(t) that we do not need to know.
Remark. Note that νQ = F implies that F is the flow into the network. If F > 0 fluid enters, and if F < 0 fluid is coming out.
Let us mention a few words about the unique solvability of the network wave model with a given boundary flow F . First of all we have not given any precise function spaces where the coefficients of the equation or the boundary flow would belong to. This means that it is not possible to find an exact reference for the solvability. On the other hand this is not a problem for linear hyperbolic problems in general.
The problem is a one-dimensional linear hyperbolic problem on various segments with co-joined boundary conditions. As the waves propagate locally in time and space, one can start with the solution to a wave equation on a single segment, as in e.g. Appendix 2 to Chapter V in [17]. Then when the wavefront approaches a junction, Equations (2.3) and (2.4) determine the transmitted and reflected waves to each segment joined there. Then the wave propagates again according to [17], and the boundary conditions are dealt with as in a one segment case. At no point is there any space to make any "choices", and thus there is unique solvability. We will not comment on this further, but for more technical details we refer to [17] for the wave propagation in a segment and the boundary conditions, and to Section 3 of [8] for the wave propagation through junctions. For an efficient numerical algorithm to the direct problem, see [18].
Boundary measurements. The area reconstruction method presented in this paper requires the knowledge of the impulse-response matrix (IRM) for all boundary points except one, which we denote by x 0 .

Definition 2.4.
We define the impulse-response matrix, or IRM, by K = (K ij ) N i,j=1 . For a given i and j we assume that H and Q satisfy the network wave model with boundary flow * for t ∈ R and x ∈ ∂G, x = x 0 . Here V 0 is the volume of fluid injected at t = 0. Then we set K ij (t) = H(t, x j )/V 0 (2.9) for any t ∈ R.
The index i represents the source and j the receiver. Note also that the IRM gives complete boundary measurement information: if ν(x)Q(t, x i ) = F (t, x i ) were another set of injected flows at the boundary, then the corresponding boundary pressure would be given by by the principle of superposition.

Area reconstruction algorithm
Strategy. Once the impulse-response matrix from Equation (2.9) has been measured, we have everything needed to determine the cross-sectional area in a tree network. As in the one pipe case [3], we will calculate special "virtual" boundary conditions, which if applied to the pipe system, would make the pressure constant in a given region at a given time. Exploiting this and the knowledge of the total volume of water added to the network by these virtual boundary conditions gives the total volume of the region. Slightly perturbing the given region reveals the cross-sectional area. * δ0(t) has dimensions of time −1 because δ0(t)dt = 1 and dt has dimensions of time.
Multiply Equation (2.1) by gA/a 2 and integrate over a time-interval (0, τ ), for a fixed τ > 0, and the whole network G. This gives if H(0, x) = 0 and mass conservation from Equation (2.4) applies. Let p ∈ G be a point at which we would like to recover the cross-sectional area. To it we associate a set D p ⊂ G which we shall define precisely later. Let us assume that there are boundary flows Q p (t, x j ) so that at time t = τ we have for some fixed pressure h 0 . Then from Equation (3.1) we have Denote the integral on the left by V (D p ). We can calculate its values once we know Q p , hence we know the value of the integral on the right. By varying the shape of D p we can then find the area A(p).
Admissible sets. The only requirement for D p in the previous section was that Equation (3.2) holds. Boundary control, e.g. as in [9], implies that there are suitable boundary flows Q p such that the equation holds for any reasonable set D p . However it is not easy to calculate the flows given an arbitrary D p . In this section we define a class of such sets for which it is very simple to calculate the flows. Let p ∈ G \ J be a non-junction point in the network at which we wish to solve for the cross-section area A(p). Since we have a tree network the point p splits G into two networks. Let D p be the part that is not connected to the inaccessible boundary point x 0 . The boundary of D p consists of let us say y 1 , y 2 , . . . , y k and p, where the points y j are also boundary points of the original network.
For each boundary point y j ∈ ∂D p , y j = p define the action time where TT gives the travel-time of waves from y j to p calculated along the shortest path in D p . Then set f (x j ) = 0 for x j ∈ ∂G \ ∂D p .
Definition 3.1. We say that D p is an admissible set associated with p ∈ G\J and with action time f , if D p and f are defined as above given p.
It turns out that with the choice of admissible sets made above we have This is because D p lies between p and the boundary points x j at which f (x j ) = 0. This gives a geometric interpretation to the set D p , i.e. that it is the domain of influence of the action times f (x j ), x j ∈ ∂G. If we would have zero boundary flows at first, and then active boundary flows when x j ∈ ∂G, τ − f (x j ) < t ≤ τ , then the transient wave produced would have propagated throught the whole set D p at time t = τ but not at all into G \ D p .
Reconstruction formula for the area. Now that the form of the admissible sets have been fixed, we are ready to prove in detail what was introduced at the beginning of this section, namely a formula for solving the unknown pipe cross-sectional area. For simplicity we assume that the wave speed is constant.
Theorem 3.2. Let p ∈ G\J be a non-junction point and D p , f the associated admissible set and action time.
the total volume of fluid injected into the network from the boundary in the time-interval (0, τ ) to create the waves H φ , Q φ . Then gives the cross-sectional area of the pipe at p.

Proof.
The assumption about the boundary flow implies that Multiply the equation by gA/a 2 and integrate τ 0 G . . . dxdt. This gives because H φ (0, x) = 0 by Equation (2.5). We will use the junction conditions of Equation (2.4) to deal with the left-hand side. But before that let us use a fixed coordinate system of the network G.
Let the pipes of the network be P 1 , . . . , P n and model them in coordinates by the segments (0, 1 ), . . . , (0, n ), respectively. On pipe P k , denote by H φ,k the scalar pressure head, and by Q φ,k the pipe discharge into the positive direction. Then is simply the discharge out of the pipe P k at the latter's endpoint represented by x = k at time t. Similarly −Q φ,k (t, 0) is the discharge out from the other endpoint, the one represented by x = 0. Now Equation (2.4) implies after a few considerations that Namely, previously we saw that the integral is equal to the sum of the total discharge out of every single pipe. But the discharge out of one pipe must go into another pipe at junctions (there are no internal sinks or sources). Hence the discharges at the junctions cancel out, and only the ones at the boundary ∂G remain. The boundary values are zero on ∂G \ (∂D p ∪ {x 0 }) by the definition of f and f + ∆ t . We also have zero initial values and a non-active boundary condition at Next, we will show that Recall that p is not a vertex. Hence it is on a unique pipe, let's say (0, ) and p is represented by x p on this segment. Furthermore assume that (or change coordinates so that) the point represented by 0 is in D p , and the one by is not. This implies that the difference of sets D f +∆t \ D f is just a small segment on (0, ). Let us look at the effect of ∆ t on D f +∆t which is defined in Equation (3.5).
The former implies that x is at most travel-time ∆ t from D p , and the latter says that it should not be in D p . In where we abuse notation and denote the cross-sectional area of the pipe modelled by (0, ) at the location x also by A(x) without emphasizing that it is the area on this particular model of this particula pipe.
Equation (3.11) gives the area at p, A(x p ), by differentiating. Let B(s) =

xp+s xp
A(x)dx. Then A(x p ) = ∂ s B(0) and the right-hand side of Equation (3.11) is equal to B(a∆ t ). The chain rule for differentiation gives Solving for the area. In this section we will show one way in which boundary values of Q can be determined so that the assumptions of Theorem 3.2 are satisfied. It is previously known that there are boundary flows giving Equation (3.2), for example by [9] where the authors show the exact L 2controllability of the network both locally and in a time-optimal way. However there was no simple way of calculating such boundary flows and the numerical reconstruction algorithm of that paper does not seem computationally efficient as it uses a Gram-Schmidt orthogonalization process on a number of vectors inversely proportional to the network's discretization size. We show that if the flow satisfies a certain boundary integral equation then this is the right kind of flow. Moreover, in the appendix we give a proof scheme for showing that the equation has a solution.
More recently various layer-peeling type of methods have appeard [10,12]. The method we present here is based on another idea, one whose roots are in the physics of waves, namely time reversibility. This is in essence a combination of the one pipe case originally considered in [6], and of a fundamentally similar idea for higher dimensional manifolds introduced in [14]. In the latter the author considers domains of influence and action times on the boundary, and builds a boundary integral equation whose solution then reveals the unknown inside the manifold. This will be our guide.
Let us recall the unique continuation principle used for the area reconstruction method for one pipe in [3,6]. Consider a pipe of length > 0, modelled by the interval (0, ). Let H and Q satisfy the Waterhammer Equations (2.1) and (2.2) without requiring any initial conditions. Then The lemma was then used to build a virtual solution H, Q satisfying also the initial conditions, and which would have H(τ, x) = h 0 for x < aτ and H(τ, x) = 0 for x > aτ for a given τ > 0. Without going into detailed proofs, the same unique continuation idea works for a tree network. The reason is that one can propagate H = 2h 0 , Q = 0 from one end of a pipe to the other, but by keeping in mind that the time-interval where these hold shrinks as one goes further into the pipe. Do this first on all the pipes that touch the boundary. Then use the junction conditions from Equations (2.3) and (2.4) to see that H = 2h 0 , Q = 0 on the next junctions. Then repeat inductively. We have shown Then for some x j ∈ ∂G.
We can now write an integral equation whose solution gives waves with Theorem 3.5. Let K ij be the impulse-response matrix from Equation (2.9).
If A is constant near each network boundary point we have for some function † k ij = k ji that vanishes near t = 0. Let p ∈ G \ J be and let D p ⊂ G be the admissible set associated with p, and with action time f . Take τ ≥ max f and let Q p (t, x j ) satisfy Then if H, Q satisfy the network wave model with boundary flow νQ p we have Proof. The claim for the impulse response matrix is standard. See for example Appendix 2 to Chapter V in [17] for the solution to the one segment setting, and then use mathematical induction and the junction conditions of Equations (2.3) and (2.4).
For the pressure, recall that the properties of the impulse response matrix from Equation (2.10) imply that for x j ∈ ∂G, x j = x 0 and 0 < t < 2τ . It is then easy to calculate that where we again used the time-symmetry of Q p . Summing all terms and using Remark 3.7. Assuming an impulse-response matrix that is measured to infinite precision and without modelling errors, one can show that Equations (3.14) and (3.15) have a solution. The idea of the proof is shown in the appendix. The solution is not necessarily unique, but it gives a unique reconstructed cross-sectional area.

4.
Step-by-step algorithm Solving for the cross-sectional area. In this second part of the step-bystep reconstruction algorithm we assume that the impulse-response matrix K has been calculated. It can be either measured by closing all accessible ends, or by measuring the system response in a different setting (i.e. different boundary conditions), then pre-process the measured signal to obtain the desired matrix K.
Once the impulse-response matrix has been measured as discussed in the previous paragraph, the following mathematical algorithm can be applied to recover the cross-sectional area inside a chosen pipe or pipe-segment in the network.

Algorithm 1.
This algorithm calculates the cross-sectional area using a discretization and Algorithm 2.
(1) Define k ij (t) for i, j = 0 by where δ ij is the Kronecker delta. (2) Choose a point p 1 in the network that is not a junction. The algorithm will reconstruct the cross-sectional area starting from p 1 and going towards x 0 until it hits the endpoint p 2 of the current pipe. and p(k + 1) divided by ∆ x . In other words This would be an equality if ∆ x were infinitesimal, or if A(p(k)) would signify the average area over the interval (p(k), p(k + 1)).

Algorithm 2.
This algorithm calculates the internal volume of the piece of network cut off by p: namely all points from which you have to pass through p to get to x 0 .
(1) For any boundary point where TT gives the travel-time between points. It is just the distance in the network divided by the wave speed.
and also simultanously set It is the internal volume of the pipe network that is on the other side of p than the inaccessible end x 0 .

Numerical experiments
All the programming here was done in GNU Octave [19].
Experiment 1 (Simple network with exact IRM). We will start by solving for the trivial area of Example 2.1. Consider this a test for the implementation of the area reconstruction algorithm. We start by calculating the impulse-response matrix of vertices A and B while C is inaccessible. D is the junction. Set A(x) = 1 m 2 everywhere, the gravitational acceleration g = 9.81 m s −2 and a wave speed of a = 1000 m s −1 .  4). Also, if a similar pressure pulse is incident to a boundary point with boundary condition Q = 0, then a pulse of the same magnitude (no sign change) is reflected. However at that boundary point the pressure is measured as 2M δ 0 . These considerations produce the following impulse-response matrix for time 0 ≤ t ≤ 1.6 s. We must have 2τ ≤ 1.   Let us define the action time functions for various points p in the network. If p ∈ AD then the action time is The travel-time from A to D is 0.4 s, and from B to D it is 0.3 s. Recall that the IRM has been measured only for time t ≤ 1.6 s. Hence we can solve for the area up to 400 m into pipe DC. Let the point p ∈ DC be given by the action time function where 0 ≤ t p ≤ 400 m/a = 0.4 s is the travel-time from D to p and recalling that we can solve the area only up to 400 m from D.  Let us apply Algorithm 1 next. The numerical implementation of Algorithm 2, makeHeq1, is in the appendix.  A plot of the cross-sectional areas is shown in Figure 3.  We simulate the IRM by an a finite-difference time domain (FDTD) algorithm with the following caveats: we use a Courant number smaller than one, simulate using a high resolution, and then interpolate the IRM to a lower time-resolution and use this as input for the inversion algorithm. This is to avoid the inverse crime [20] which makes inversion algorithms give unrealistically good results when applied on data simulated with the same resolution or model as the inversion algorithm uses. For numerical reasons instead of sending a unit impulse Q = νδ we send a unit step-function, and then differentiate the measurements with respect to time. We will not show the implementation of FDTD and the self-explanatory functions removeInitialPulse, medianSmooth and differentiate because they are not the focus of this already rather long article. The main point is the inversion algorithm. The simulation gives us the following response matrix, as shown in Figure 5.

add t o V t h e volume o f w a t e r from each a c c e s s i b l e end t h a t would have gone INTO t h e p i p e t o make H=1 a t t=t a u :
Now solving the inverse problem has the same logic as in Experiment 1. There are two differences, a large one and a small one. The former is that the action time functions are of course different. This is what encodes the network topology for the inversion algorithm. And secondly we must use quite a lot of regularization when solving for the area of ED. This is because Then applying Algorithm 1 gives the solution as before.
100  The solution is displayed in Figure 6. The gray uniform line represents the original cross-sectional area. The dashed line is the solution to the inverse problem.

Discussion and conclusions
We have developed and implemented an algorithm that reconstructs the internal cross-sectional area of pipes, filled with fluid, in a network arrangement. The region where the area is reconstructed must form a tree network and the input to our algorithm is the impulse-response measurements on all of the tree's ends except for possibly one. The algorithm involves solving a boundary integral equation which is mathematically solvable in the case of perfect measurements and model (see the appendix). We wrote a step-bystep reconstruction algorithm and tested it using two numerical examples. The first one has a perfectly discretized measurement, and the second one has a discretization and numerical error introduced on purpose. Even with these errors, which are very typical when doing real-world measurements, Tikhonov regularization allows us to reconstruct the internal cross-sectional area with good precision. However this is just a first study of this algorithm, and a more in-depth investigation would be required for a more complete picture of its ability in the case of noise and other errors in the data. The theory is based on our earlier work [3] on solving for the area of one pipe, and x (m) Figure 6. Solved cross-sectional areas of Experiment 2 on an iterative time-reversal boundary control algorithm [14] in the context of multidimensional manifolds. We assumed in several places that the wave speed a(x) is a known constant. What if it is not? We considered this situation for a single pipe in our previous article, [3]. In there, if the area is known and the speed is unknown, then the algorithm can be slightly modified to determine the wave speed profile along the pipe. If both the wave speed and the area are unknown then it can determine the hydraulic impedance Z = a/gA as a function of the travel-time coordinate. However this is not so straightforward for the network case because of the more complicated geometry. What one should do first is find an algorithm to solve this problem: the wave speed is constant and known, the area is unknown, the topology of the network is known, but the pipe lengths are unknown. We leave this problem for future considerations as this article is already quite long. But it is an important question, because in fact in some applications the anomaly is the loss in pipe thickness rather than the change in pipe area and this is often revealed by finding the wave speed along the pipes assuming that the area remains unchanged [21].
In Experiment 1 we did not use regularization and the solution was perfect, as shown in Figure 3. However here we had the simplest tree network, the simplest area formula, and perfectly discretized measurements.
In Experiment 2 we simulated the measurements using a finite difference time domain (FDTD) method with Courant number smaller than one. x (m) Figure 7. Solving without regularization Using Tikhonov regularization was essential, and produced the reconstruction shown in Figure 6. Without regularization one would get a large artifact, as in Figure 7. The reasons for using imperfect measurements are to demonstrate the algorithm's stability. We could have calculated the impulseresponse matrix using the more exact method of characteristics. However in reality, when dealing with data from actual measurement sensors, one would never get such perfect inputs to the inversion algorithm. First of all there are modelling errors, and secondly, measurement noise. By using FDTD with a small Courant number and also by inputting measurements with a rougher discretization than in the direct model we purposefully seek to avoid inverse crimes, as decribed in Section 1.2 of [20]. An inverse crime happens when the measurement data is simulated with the same model as the inversion algorithm assumes. Typically in these cases the numerical reconstruction looks unrealistically good, even with added measurement noise, and therefore does not reflect the method's actual performance in real life, where models are always approximations.
Several directions of investigation still remain open. A more in-depth numerical study should be concluded, with proper statistical analysis and for example finding the signal to noise ratio that still gives meaningful reconstructions. There is also the issue of measuring the impulse-response matrix. It would be very much appreciated if there was no need to actually close almost all the valves in the network to perform measurements. Hence one could investigate various other types of boundary measurements and see how to process them to reveal the impulse-response matrix used here. To make even further savings in assessing the water supply network's condition, one could also try implementing the theory from [10] as an algorithm. Their theoretical result only requires the so-called "backscattering data" from the network: instead of having to measure the full impulse-response matrix (K ij ) N i,j=1 it would be enough to measure its diagonal (K ii ) N i=1 . This means that the pressure needs to be measured only at the same pipe end as the flow is injected. So the exact same type of measurement as is done for a single pipe, as in [3], should be done at each accessible network end. For a network with N + 1 ends this translates into N measurements compared to the N 2 for the full IRM.
Lastly we comment on our algorithm. What is interesting about it, is that there is no need to solve the area in the whole network at the same time. To save on computational costs one could for example reconstruct the cross-sectional area only in a small region of interest in a single pipe even when the network is otherwise quite large. Alternatively one could parallelize the process and solve for the whole network very fast using multiple computational cores. These would likely be impossible when having only the backscattering data described above.
The steps in Algorithms 1 and 2 describe the reconstruction process in detail. In essence only a single matrix inversion (or least squares or other minimization) is required for solving the area at any given point in the network. The size of this matrix depends on how finely the impulse-response matrix has been discretized, and how far from the boundary this point is. The rougher the discretization, the faster this step. However having a poor discretization leads to a larger numerical error, and one might then need to use a regularization scheme. All the numerical examples in this article were performed with an office laptop from 2008, and they took a few minutes to calculate including simulating the measurements. The algorithm is computationally light, simple to implement, and thus suitable for practical applications.

Acknowledgements
This research was partly funded by the fullowing grants.
• T21-602/15R Smart Urban water supply systems(Smart UWSS), and • 16203417 Blockages in water pipes: theoretical and experimental study of wave-blockage interaction and detection. Furthermore all authors declare that there are no conflicts of interest in this paper.

Existence of a solution.
We give a list of the steps involved in proving that the equation in Theorem 3.5 has a solution when the impulseresponse matrix has been measured exactly and there is no modelling errors. Since the ultimate goal for our work is on assessing the quality of water supply network pipes, in the end there will be both modelling errors and measurement noise, so we will not give a complete formal proof. However we give enough details so that any mathematician specialized in the wave equation can fill the gaps after choosing suitable function space classes and assumptions for the various objects.
(1) Given p ∈ G \ J which defines the admissible set D p and action time f according to Definition 3.1, and with τ > max f and h 0 > 0, we want to solve We extend boundary data from the interval (0, τ ) to (−∞, +∞) and also absorb the interior boundary normal, by writing where is equivalent to solving the equation of Item 1 in the corresponding L 2 -based space. We equip L 2 0 with the inner product where the complex conjugation can be ignored because all of our numbers are real. (4) The equation in Item 3 can be written ash 0 = K Q p in L 2 0 . Here we define the operator K by and it is a well-defined operator if k ij is a distribution of order 0 which is the case if the cross-sectional area function A is piecewise smooth for example. It indeed maps L 2 0 → L 2 0 satisfying the support and time-symmetry conditions. (5) We will show that K is a Fredholm operator L 2 0 → L 2 0 that is positive semidefinite. Define The set Ω coincides with the set where unique continuation from the boundary holds, as in Proposition 3.4. One can show that in t ∈ R, x ∈ G \ J. (10) By Items 6, 8 and 9 we see that (11) By Items 9 and 10 and the divergence theorem we have where ν t,x is the external unit normal vector to Ω at (t, x) ∈ ∂Ω.
(12) Let us split ∂Ω next. By Item 6 we see that (13) On the first set in Item 12 the external unit normal ν t,x j is given by ν t,x j = (0, −ν(x j )) because ν was defined as the internal normal at the pipe ends. (14) Consider the second set in Item 12. On each individual pipe segment P ⊂ D p the map x → t + (x) is linear (affine). How does t + change when we go from x to x + ∆x? Recall that |∆x| = a |∆t|, and since (t + (x), x) follows the characteristics, so t + (x + ∆x) = t + (x) ± ∆x/a. Its value decreases if ∆x > 0 and the positive direction of the coordinates x point towards p. Both of these follow from Item 6 and the definition of the action time function f . By the definition of µ in Item 9 we see that ∆x Thus the normal has a slope of µ/a and so (15) By Items 7, 8 and 10 to 14 we get by the support condition of F ∈ L 2 0 given in Item 3. In detail we see that the D p -integrals vanish because of the following. Firstly we see that (a, µ) · (−µ/a, 1) = 0 making the first integral over D p vanish. Secondly, by Item 8, (QH)(t − (x), x) = 0 for x ∈ D p . The integral over the singleton (t, s) = (τ, p) is zero because QH is a function since F is too. (16) Recall Equation (2.10) and the splitting of the IRM K into the impulse and response matrices in Theorem 3.5. Thus (17) Combining Items 15 and 16 and changing the order of integration in double integrals, we arrive at for arbitrary F ∈ L 2 0 and where S, H, Q are defined by Items 7 and 9. (18) We cannot show that K F = 0 implies F = 0 unless the network is a single pipe. For a counter example choose a three-pipe network with constant unit area and wave speed, and each pipe is of unit length. Then send a pulse from one end and the same with opposite sign from another end. These pulses meet at the junction at the exact same time and cancel out the pressure there, so nothing propagates to the third pipe. Then they continue on the original two pipes until they hit the boundaries. The boundary flow F should be then chosen so that these pulses are absorbed completely instead of reflecting back. This boundary flow F gives K F, F = 0 by Item 17, but F = 0. (19) The existence of a solution follows from the following: that K is self-adjoint and Fredholm, and thath 0 ⊥ ker K , so then there is F giving K F =h 0 . We shall show these next. (20) If we assume that A(x) is smooth enough, then the components of the response matrix (k ij ) N i,j=1 have at most a finite number of delta-functions arising from the junctions of G (in the time-interval (0, 2τ )), and are otherwise a function. Hence K can be split into a multiplication by nonvanishing functions, a finite number of translations that are also multiplied by such functions, and a smoothing integral operator which is compact L 2 0 → L 2 0 . The identity and translations multiplied by non-vanishing functions are Fredholm operators, and the rest is compact. Hence the whole K is Fredholm. The space L 2 0 is Hilbert, and K is easily seen to be symmetric there because k ij = k ji . It is a bounded operator. Hence it is self-adjoint. Thus K is a self-adjoint Fredholm operator. Henceh 0 ⊥ ker K = ker K * , the latter because of self-adjointness, soh 0 ∈ (ker K * ) ⊥ = ran K = ran K because the range of a Fredholm operator is closed. In other words there is F ∈ L 2 0 such that K F =h 0 , and thus there is a solution to the equations in Theorem 3.5.
8.2. Numerical reconstruction algorithm. We describe here the numerical implementation of Algorithm 2. Its inputs include tHist, K, the discretized time-variable vector and cell containing the discretized components of the response matrix k = (k ij ) N i,j=1 and f which is a vector modeling the discretized action time function. Other inputs are a0, A0, which are vectors containing the wave speed and cross-sectional area at the boundary points of the network, and tau, g, reguparam are scalars representing τ , the constant of gravity g, and the Tikhonov regularization parameter used in the inversion of Equation (4.1). Its output, Qtau, is a cell whose components are discretizations of t → Q p (t, x j ) that solved Equation (4.1) for the various boundary points x j . We write first the N (=nBndryPts) integral equations containing the N response function components. This shall be indexed by j (pipe end where pressure measured) and i (pipe end where pulse sent from). Recall the equation, for j = 1, . . . , N and τ − f (x j ) < t ≤ τ . But recall that we must still set q(t, x j ) = 0 when t ≤ τ − f (x j ). Our numerical implementation assumes that h 0 = 1. We will write the equation above as H*q = RHS where q is a block vector where each block is indexed by i and q i = [q(dt, x i ); q(2 dt, x i ); ...; q(M dt, x i )] with M = τ /dt . Each component of RHS is either h 0 = 1 or 0. We do a piecewise constant approximation of all of these functions and equations. We discretize t=l*dt, s=k*dt with l,k=1,...,M. Set the tolerance for floating point numbers next. If two numbers are at most tolerance from each other, then they are considered the same. Without doing this, some numerical errors might happen when comparing two floating point numbers that are close to each other. This would cause errors of the order of dx in the area reconstruction. We can now start discretizing the integral equation. Recall that q(s, x i ) must be zero when s ≤ τ − f (x i ). Hence the columns of the matrix that hit the indices corresponding to s ≤ τ −f (x i ) should be zero. Similarly, when t ≤ τ −f (x j ) we want the equation to give q(t, x j ) = 0, so the part coming from the integral must be zeroed. Then add the "identity"-type part of the integral equation, and save the block into the final matrix. Finally split the vector q into a cell whose components correspond to the boundary flows at the different boundary points.