It is natural to represent a Water Distribution System (WDS) in a graph-theoretic setting, where features such as, e.g., pipes, intersections, house services, tanks, pumps, and valves are effectively represented by means of nodes, edges and their related properties. Unfortunately, precisely modeling an operational WDS is a non-trivial and intricate task, mostly attributable to the non-stationarity of the system inputs. Indeed, a considerable number of variables and parameters involved in the model may have changed during the WDS lifespan because of, e.g., pipe degradation, valve malfunctioning, varied demand, added connections, etc., and the model developed in the design phase may not be accurate forcing a re-development of the simulation model. Updating the parameters of the model is therefore problematic, especially because of the frequent lack of available data, and aiming at recovering a detailed physics-based operational model of the WDS is very often an unrealistic task. In this talk, we discuss a surrogate non-linear model that is capable of accurately reproducing the whole water dynamics in the WDS despite being ruled by just two vectors of free parameters, which characterize the underlying weighted graph p-Laplacian operator. To construct this surrogate model, we formulate a TV-regularized p-Calderon inverse problem to identify the weight and exponent distributions, and formulate a gradient flow scheme for its numerical solution. A complex real aqueduct characterized by 138 km of pipes is used to verify the applicability of the proposed surrogate model and assess its accuracy in reproducing the real behavior.