Fractal Analysis of the Optimal Hydraulic Gradient Surface in Water Distribution Networks

The design of water distribution networks (WDNs) may be defined as the selection of a certain configuration of commercially available diameter sizes, in such a way that a group of hydraulic constraints are satisfied and a set of operational conditions are achieved (Bragalli et al. 2012; Qiu et al. 2021). In its traditional formulation, the optimized capital cost design of WDNs seeks to establish the least-cost set of diameter sizes that ensures the supply of nodal flow demands and a suitable water pressure without violating the law of conservation of mass and energy (Geem 2008; Ezzeldin et al. 2013; Saldarriaga 2016). However, additional criteria can be included, for example, water quality, reliability, and resilience (Prasad and Park 2004; Bi and Dandy 2014; Moosavian and Lence 2020; Sirsant and Reddy 2020). The complexity of the problem could be attributed to the nonlinear relationship between head loss and discharge (Saldarriaga 2016), the discrete nature of the available pipe diameters, and the extensive quantity of feasible solutions (Kadu et al. 2008). Consequently, the optimization of WDNs is a highly indeterminate problem (Reca and Martínez 2006; Saldarriaga et al. 2013).

Numerous stochastic and deterministic strategies have been applied for WDN minimum-cost optimization. The first approaches involved the use of traditional deterministic gradient-based techniques, such as linear and nonlinear programming (Alperovits and Shamir 1977; Lansey and Mays 1989), that are prone to local-optima entrapment and culminate in split-pipe sizes (Lin et al. 2007; Zheng et al. 2011; Qiu et al. 2021). Evolutionary algorithms (EAs) made it possible to overcome these disadvantages (Zheng et al. 2011). The use of metaheuristic alternatives offered a broader exploration of the search space by developing guided stochastic procedures with the use of learning functions based upon fitness measures (Saldarriaga et al. 2013; Marchi et al. 2014) yet tackling the problem from a mathematical perspective without considering the underlying fluid dynamics inside a WDN (Saldarriaga et al. 2010, 2015). One of the most prominent EAs are genetic algorithms (GAs), which are bioinspired procedures that recreate natural selection (Sivanandam and Deepa 2008). GAs generate an initial random population of possible solutions (individuals) that undergo an iterative process of recombination and mutation until a satisfactory configuration of diameter sizes is found (Dandy et al. 1996; Savic and Walters 1997; Montesinos et al. 1999; Mitchell 1999). Based on their fitness, the best individuals are selected to yield an improved offspring in terms of the objective function (Sivanandam and Deepa 2008). Note that the technique is applicable for demand-driven models, i.e., networks in which the nodal demands are known and independent of the hydraulic behavior of the system, or pressure-driven models, i.e., networks in which the nodal flow rates show a dependency on the available pressure (Farmani et al. 2007; Páez et al. 2014). Simulated annealing (Cunha and Sousa 1999; Reca et al. 2008; Cunha et al. 2019), harmony search (Geem 2008), and Tabu search (Sung et al. 2007; Cunha and Ribeiro 2004) are additional examples of metaheuristic approaches.

Within this framework, hydraulic-based power use methodologies introduced the physical phenomena inside the optimization process and lessened the computational effort by reducing the number of iterations to reach convergence (Saldarriaga et al. 2014, 2015). One of the first approaches in this matter was documented by Wu (1975) for main drip irrigation lines. Wu (1975) observed that the minimum-cost diameter size configuration in a system of pipes in series with constant and known flow demands could be obtained by predefining a parabolic hydraulic gradient line (HGL) with a sag of 15% of the total head loss. As a result, objective head loss values could be calculated for each pipe, providing enough information to determine the desired diameters (Wu 1975; Saldarriaga et al. 2010). Subsequently, Featherstone and El-Jumaily (1983) extended Wu’s criterion to looped networks, and the idea was further developed when the optimal power use surface (OPUS) methodology was introduced (Saldarriaga et al. 2010). OPUS requires the calculation of the optimal hydraulic gradient surface (OHGS), a discontinuous three-dimensional (3D) geometrical entity formed by the nodal piezometric head values that guarantee the calculation of a least-cost set of diameter sizes (Saldarriaga et al. 2005, 2010, 2013, 2020). In a similar manner to Wu’s (1975) observations, the OHGS establishes the way in which the available energy should be spent within the system to ensure an optimal design while satisfying the minimum water pressure and nodal flow demands (Saldarriaga et al. 2005, 2010, 2013). Although the resulting initial design is given in continuous pipe sizes, OPUS executes rounding procedures so that the final design is consistent with the commercial diameters and meets the hydraulic constraints (Saldarriaga et al. 2010, 2020). Several tests with benchmark networks have shown that OPUS can obtain high-quality solutions with very slight differences when compared to the minimum cost solutions reported in the literature (Saldarriaga et al. 2010). Nevertheless, it should be noted that optimization methodologies that use hydraulic-based power rely on the assumption that Wu’s criterion is fulfilled (Saldarriaga 2016). In pressure-driven models, the flow demands are unknown prior to the optimization of the WDN, thereby impeding the applicability of OPUS without supplemental adjustments (Saldarriaga 2016). In this regard, integer linear programming has been combined with some main features of OPUS, showing successful results (Páez et al. 2013, 2014). Furthermore, hydraulic-based methodologies may also be used for preconditioning purposes or as feedback to improve the performance of metaheuristic algorithms as shown by Páez et al. (2020) and Liu et al. (2020).

From a mathematical perspective, the OHGS may encompass interesting properties that could provide additional details of the underlying physics inside the optimization process (Jaramillo 2020). In terms of its structure, the OHGS exhibits a highly irregular appearance that suggests that it may not be described by Euclidean geometry, which is suitable for the characterization of regular objects. Therefore, given the mathematical complexity of the OHGS, a non-Euclidean approach, such as fractal analysis, may be used to provide a characterization of its main geometrical properties. Several works have stated that many of the amorphous surfaces found in nature or in biological and engineering applications can be analyzed under the postulates of fractal theory (Mandelbrot 1983; Zhou and Lam 2005; Ju and Lam 2009; Persson 2014; Spasic 2014). Treating the OHGS as an approximate fractal body widens the spectrum of possible examination viewpoints. For this purpose, it is necessary to understand the main concepts that are generally used when involving fractals. In a few words, a fractal can be defined as a complex and highly irregular mathematical entity whose geometrical characteristics are self-similar at any arbitrary scale (Falconer 1990). Although perfect fractals must show this property, the fractal behavior can be present in a statistical or approximate manner in real-life phenomena (Mandelbrot 1977; Falconer 1990). Measuring these entities is usually done by estimating the fractal dimension, a noninteger parameter that describes how a fractal fills space and its degree of irregularity (Falconer 1990; Jian-Hua et al. 2009). Different definitions of fractal dimension have been used in the literature, such as the Hausdorff, packing, or box-counting dimension (Edgar 2008). The latter is appealing because of its simplicity and applicability for computational purposes (Falconer 1990, 1997). In brief, the box-counting fractal dimension is based on the smallest number of sets of size $\delta$ that are required to cover the fractal (Falconer 1997). It has been stated that a scaling or power law relationship may exist between the number of sets $M_{\delta }(F)$ and the scale of analysis $\delta$, involving the fractal dimension as the exponent (Mandelbrot 1977; Falconer 1990, 1997):where $D$ = fractal dimension; $M_{\delta }(F)$ = measurement of any fractal $F$ at arbitrary scale $\delta$; and $\delta$ = scale of analysis. Plentiful numerical strategies have been adopted in the field of landscape, remote sensing, and physical geography to estimate the box-counting fractal dimension of rugged surfaces using this concept (Clarke 1986; Zhou and Lam 2005; Sun et al. 2006). Some examples of these techniques include the triangular prism, differential box counting, variogram, robust fractal estimator, probability, and variation estimator methods (Clarke 1986; Sarkar and Chaudhuri 1992; Lam et al. 2002; Zhou and Lam 2005; Sun et al. 2006; Ju and Lam 2009). Furthermore, fractal behavior is not restricted to rough surfaces. It has been shown that a WDN may show fractal properties in its topological structure that may be applied in the analysis of the resilience, reliability, water quality, and operational conditions, as well as in classification, risk assessment, and sectorization procedures (Diao et al. 2017; Di Nardo et al. 2017; Vargas et al. 2019; Kowalski et al. 2019; Iwanek et al. 2020; Kowalska et al. 2020; Giudicianni et al. 2021). In this context, WDNs have been examined employing graph theory concepts by modeling their underlying topological structure as directed or undirected graphs whose vertices must be covered with the smallest number of boxes or clusters of a certain size (Diao et al. 2017; Vargas and Saldarriaga 2019). In this case, $M_{\delta }(F)$ is understood as a set of junctions or nodes that meet certain conditions, such as proximity (Song et al. 2007; Diao et al. 2017). It should be noted that the calculation of the minimum required covering units $M_{\delta }(F)$ is usually achieved with box-covering algorithms (Diao et al. 2017). Some of the most prominent box-covering methodologies include greedy and burning algorithms that differ in the way and order in which the nodes of a network are packed (Song et al. 2007; Wen and Cheong 2021). Moreover, the calculation of the minimum number of covering units can be restricted to involve the topological properties of the network only, such as the way vertices are linked, or include hydraulic criteria, such as the flow distribution or nodal piezometric heads (Kim et al. 2007; Song et al. 2007; Vargas and Saldarriaga 2019; Vargas et al. 2019).

In this study, a numerical modeling methodology was used to conduct the fractal analysis of several benchmark and Colombian WDNs and their corresponding OHGS. For this purpose, a least-cost design and several nonoptimal designs were generated for each system. Subsequently, the designs were analyzed using the variation estimation method and a greedy box-covering algorithm to compute the fractal dimension of the OHGS and of the underlying structure of each WDN design, respectively. In terms of the OHGS, the fractal dimension is employed as a spatial metric that translates the complex morphology of the surface into numerical data that can be used to understand the relationship between the geometrical properties of the OHGS and the inherent energy dissipation phenomenon of a least-cost design. In relation to the underlying structure of the WDNs, the fractal dimension is used as a parameter that reveals the degree of irregularity in the topology, flow, and energy distribution inside a given network to characterize the intrinsic mathematical properties of the hydraulic behavior of each design. The use of several nonoptimal designs was of interest to analyze the mathematical and physical differences between the hydraulic behavior of a least-cost design and a more expensive solution. Finally, this study also intended to consider the applicability of fractal analysis as a novel approach to solving WDN capital cost optimization problems and to guide operational improvement procedures. In the context of optimization, the fractal analysis of a hydraulic gradient surface (HGS) may be used to guide the calculation of pipe diameters by considering the inherent fractal properties in the hydraulic behavior of a least-cost design. In the case of WDN operational improvement, the fractal dimension may be used as a parameter that suggests the proper fractal properties of the hydraulic behavior of a network for achieving better performance in terms of operational cost, reliability, or water quality.

The aforementioned methodology was applied to 28 WDNs. Thirteen of the chosen systems have been widely studied in the literature and are usually conceived as benchmark networks with a well-known hydraulic behavior and a reported least-cost design. The remaining 15 networks are real water supply infrastructure located in various Colombian towns and cities that were analyzed in previously published dissertations and research projects. All case studies are treated as new water supply systems in which pipe diameters need to be calculated. The networks were chosen in such a way that different topographical, topological, and energy availability conditions were tested so as to identify differences in the fractal behavior. Therefore, the 28 selected systems include demand- and pressure-driven models with flat and steep topographical conditions, one or more reservoirs, and branched and looped structures. Nevertheless, it should be noted that the analysis was restricted to gravity-driven networks, i.e., without pumping devices. As stated by Saldarriaga et al. (2020), the OPUS methodology was originally made to be applied to gravity-driven WDNs. Coding adjustments of the original algorithm are required if networks with pumps will be included in the analysis (Saldarriaga et al. 2020). Consequently, so as to maintain the current OPUS version unaltered, the analysis was only conducted in gravity-driven networks. The main properties of each network are presented in Table 1.

Each of the analyzed networks was modeled using an EPANET data file that contained all the necessary information about the spatial locations of nodes, pipes, valves, and reservoirs, as well as minor loss coefficients, emitter coefficients, maximum nodal flow rate values, and the available head on each reservoir. Most of the EPANET data files corresponding to the benchmark networks were retrieved from Wang et al. (2015), available online in the University of Exeter Centre for Water Systems. EPANET Colombian network data files were retrieved from CIACUA. Although many of the benchmark networks were originally conceived using a Hazen–Williams empirical equation to describe frictional losses, Darcy–Weisbach and Colebrook–White physical equations were preferred, as described in the proposed hydraulic calculation methodology. In this regard, an absolute roughness of 0.0015 mm was assumed for all networks, excluding Balerma in which a roughness of 0.0025 mm was used according to the literature. Furthermore, even though some of the selected benchmark networks were originally conceived to include pipes with fixed diameters and a maximum and minimum pressure requirement, every model was optimized assuming that every diameter had to be found in a completely new system and that the minimum pressure was the only piezometric head constraint. Lastly, it is important to mention that the unit cost information used to specify the coefficient ($K$) and exponent ($n$) in the capital cost equation [Eq. (2)] was obtained from different sources. For benchmark networks, the information was retrieved from the literature, whereas for Colombian networks the unit costs were retrieved from PAVCO Wavin, a pipe manufacturing company (Bogotá, Colombia). All the unit cost information is summarized in Tables S1–S14, the corresponding cost curves are showcased from Figs. S1–S13, and the resulting $K$ and $n$ values for each network are presented in Table 1.

The optimization process culminated in 28 least-cost designs, one for each WDN, and 280 nonoptimal designs, 10 for each WDN. An EPANET data file was generated for each design and is available in Jaramillo (2022). All solutions were established in agreement with the minimum pressure constraint and the fulfillment of the nodal flow demands. In addition, nonoptimal solutions were chosen to ensure a progressive increase in the capital pipe cost; therefore, nonoptimal solutions were found near the achieved least-cost design as well as much more expensive designs. As an example, the evolution of Toro network capital cost is presented in Fig. 6. The remaining capital cost information can be consulted in Tables S15–S17. Given that some modifications were carried out in the problem formulation of many of the analyzed benchmark networks in terms of the friction equation and the involved constraints, the achieved least-cost design may differ if compared with the reported record optimal design. The mean and standard deviation were also calculated for the set of diameter values of each design. The general tendency of the quantified statistical parameters is exemplified in Fig. 7 for some of the studied networks. The figure shows that the mean and standard deviation of the diameter values were generally higher in the nonoptimal designs. All mean and standard deviation values are included in Tables S18 and S19. Furthermore, Figs. 8–10 exemplify the application of the interpolation algorithm that computes the HGSs by illustrating the cases of R28, Cazucá, and New York Tunnels networks.

The calculated fractal dimensions of the OHGSs are summarized in Table 2 with their corresponding coefficients of determination from the least-squares linear regressions over the analyzed log $V(\epsilon )$ versus log $\epsilon$ data sets. Fig. 11 exemplifies the application of the fractal analysis over the OHGS and nonoptimal HGSs of El Overo network by illustrating the approximate linear behavior between log $V(\epsilon )$ and log $\epsilon$. The remaining log $V(\epsilon )$ and log $\epsilon$ graphs are presented from Figs. S14–S40. A box and whisker plot (Fig. 12) is presented to display the fractal dimension values of the OHGSs and of the nonoptimal HGSs for each of the analyzed networks and to depict the range and variability in the corresponding data sets. Eleven fractal dimension values were calculated for each WDN: one for the optimal design and one for each of the nonoptimal solutions. Fig. 12 also illustrates the interval where the 11 coefficients of determination were found for each WDN, denoting the maximum and minimum $R^2$ coefficients from the 11 analyzed HGSs in each case. If the $R^2$ coefficients were all the same in a specific network, a single $R^2$ value is specified instead of an interval. Both Table 2 and Fig. 12 show that the OHGSs and nonoptimal HGSs are fractal objects given that the $R^2$ coefficients are, in general, very close to 1.00; thus, the results reveal that the nature of the energy dissipation phenomenon in a WDN has an approximate self-similar behavior that can be characterized by the fractal dimension. In terms of the variability, Fig. 12 reveals little statistical dispersion in the data sets considering that outliers were not common and that the interquartile range (IQR) was narrow. If the computed fractal dimensions are compared with the theoretical feasible range between 2.00 and 3.00 (Bisoi and Mishra 2001), it is possible to identify that the estimated values are slightly below the lower limit in a few of the studied networks. This issue can be attributed to the numerical performance of the chosen algorithm and should be addressed in future work to avoid underestimating the fractal dimension of the surfaces. All fractal dimension values and coefficients of determination for each WDN can be individually consulted in Tables S20 and S21.

Fig. 12 also suggests that each WDN can be classified into two groups depending on the fractal dimension of their corresponding OHGS. On the one hand, it was found that in 64.29% of the case studies, the fractal dimension of the OHGS was greater than or approximately equal to the third quartile of the data set. This means that the fractal dimension of the OHGS was greater than at least 75% of the fractal dimension values of the nonoptimal HGSs in these cases. Two Loops, Two Reservoirs, Taichung, Jilin, Hanoi, Blacksburg, Fossolo, R28, Pescara, Modena, El Overo, San Vicente, Cazucá, Toro, Bugalagrande, Carmen del Viboral, Morrorico Bajo, and La Enea are among this first group. On the other hand, it was found that in 35.71% of the case studies the fractal dimension of the OHGS was between the first quartile and the median of the data set. This means that the fractal dimension of the OHGS was only greater than 25%–50% of the fractal dimension values of the nonoptimal HGSs in these cases. New York Tunnels, BakRyan, Balerma, La Uribe, Elevada, Andalucía Alta, La Cumbre, Andalucía Baja, Candelaria, and Chinú are among this second group.

Moreover, the fractal dimension of the underlying structure of the optimized networks, for each calculation criterion, are also presented in Table 2 with their corresponding coefficients of determination from the least-squares linear regressions over the analyzed log $N_B$ versus log $l_B$ data sets. Table 2 shows that the underlying structure of the optimized WDNs exhibits a fractal behavior since a strong linear relationship was found between log $N_B$ and log $l_B$ for the four calculation criteria. Although the $R^2$ coefficients were not always over 0.95, this was not interpreted as a complete absence of the fractal behavior but as the presence of self-similarity to a lesser extent. Fig. 13 shows the relationship between the fractal dimension of the OHGSs and the fractal dimension of the corresponding optimized networks, including Pearson correlation coefficients and distance correlation coefficients to characterize the statistical dependency between both quantities. The distance correlation coefficients were calculated as explained by Székely et al. (2007). Two important findings can be identified in Fig. 13. First, although the calculated Pearson correlation coefficients showed a weak linear relationship between the two parameters, a positive correlation was also identified. Second, since the distance correlation coefficients were greater than zero, neither quantity was independent, and a nonlinear relationship existed between the two.

Lastly, the fractal dimension of the underlying structure of the nonoptimal designs was calculated involving the flow, energy, and infrastructure criterion. Since the computed fractal dimension by the topological criterion does not change between the optimal and nonoptimal designs for a given network, this criterion was not used to analyze the nonoptimal solutions. Table 3 presents a summary of the main features of this data set by including the cumulative frequency at the fractal dimension of the optimized network, i.e., the number of nonoptimal designs in which the fractal dimension of their underlying structure is less than or equal to the fractal dimension of the optimized network, the standard deviation of the data set, and the minimum and maximum computed fractal dimensions for each WDN to illustrate the statistical range of the data. The cumulative frequencies show that in 64.29%, 75%, and 57.14% of the case studies, the fractal dimension of the optimized network was greater than or equal to at least six of the fractal dimension values of the nonoptimal networks when using the flow, energy, and infrastructure calculation criterion, respectively. The individual fractal dimension values and coefficients of determination for each WDN design are presented in Tables S22–S27.

Some important points can be made based on the trends depicted in Fig. 7. As for the mean, the observed behavior was expected considering that a more expensive design implies larger diameters. In relation to the standard deviation, the observed behavior can be explained considering that the calculated least-cost designs consist of a homogeneous set of smaller available diameters, while the nonoptimal solutions exhibit a greater dispersion considering that the corresponding sets of diameters involve a mix of small and large diameter values as randomly selected by the GA. The only two exceptions were Hanoi and Jilin. In these cases, the standard deviation was smaller in the nonoptimal designs considering that these solutions consisted of homogeneous sets of the larger available diameters.

Relevant findings were identified in the fractal analysis of the HGSs (Table 2, Figs. 11 and 12). In this context, the presence of self-similarity in the HGS means that the energy dissipation pattern that characterizes a specific WDN design tends to repeat itself throughout the network regardless of the scale on which the networks are analyzed. A greater fractal dimension in this context means a rougher and more tortuous energy dissipation pattern. In terms of the WDN classification suggested by Fig. 12, the observed behaviors can be explained as follows. On the one hand, for the cases in which the fractal dimension of the OHGS was generally greater than the fractal dimension of the nonoptimal HGSs, the OHGSs showed a higher degree of irregularity in accordance with the nature of the energy use in a least-cost design. Given that an optimal design is characterized by the predominance of the smaller available diameters and considering that the friction slope in a pipe is inversely related to its diameter (Saldarriaga 2016), it was expected that steeper hydraulic gradients would be found in the optimal solutions considering that smaller diameters induce higher frictional losses. As a result, the OHGSs exhibited a rougher and more tortuous morphology with sharper peaks, concavities, and abrupt head losses that translate into a higher fractal dimension. Figs. 8 and 9 display these observations for R28 and Cazucá networks by depicting the attenuation of the irregular morphology between the OHGS and the HGS of the most expensive design due to the increase of the pipe diameters and the consequent decrease of the energy losses throughout the network.

On the other hand, for the cases in which the fractal dimension of the OHGS was generally smaller than the fractal dimension of the nonoptimal HGSs, four features of these systems may be involved in the observed results. First, an elongated and stretched appearance was typical in the topological structure of these networks. Second, the nodal flow demands and pipe diameters in the corresponding designs induced the generation of HGSs with a flat structure, like a geometrical plane, as exemplified in Fig. 10. It is plausible that the performance of the variation estimation method may be compromised when no significant peaks or concavities are present in the HGSs. Third, some of these networks had an atypical connectivity, for instance, New York Tunnels and La Uribe. Finally, some of these networks had an atypical flow pattern, for instance, Balerma and Elevada, with the same flow demands in almost all nodes.

The fractal analysis of the underlying structure of the WDNs (Table 2 and Fig. 13) provided additional details to discuss. The identified self-similar behavior, as shown by the fractal dimensions of the network and the $R^2$ coefficients, can be understood as the presence of scale invariance in the spatial arrangement of nodes and pipes and how they are connected (topological criterion) (Di Nardo et al. 2017), in the way in which the available discharge (flow criterion) and piezometric head (energy criterion) were distributed throughout the networks, and in the spatial distribution of the calculated diameters (infrastructural criterion). Another aspect to consider is that a greater fractal dimension in this context should be understood as a more dispersed spatial distribution of the nodes and pipes, a more uneven distribution of the available flow and energy or a more heterogeneous set of pipe sizes, depending on the chosen calculation criterion. In addition, the positive correlation between the fractal dimension of the OHGSs and the fractal dimension of the corresponding optimized networks, as shown by the Pearson coefficients, imply that an increase in the irregularity of the cost-optimized network, in terms of its hydraulic behavior, spatial arrangement, and its diameter values, translates into a rougher energy dissipation pattern, as depicted in the corresponding OHGS. In terms of the distance correlation coefficients, the existence of a nonlinear relationship shows that the fractal behavior in the topological structure, flow, energy, and diameter size distribution in the optimized network somehow influences the energy dissipation phenomenon exhibited in the OHGS.

These findings provide a hint as to some possible applications of the fractal analysis of HGSs in the context of WDNs. First, the fractal dimension may be used as an additional fitness measure to include in new or existing WDN capital cost-optimization algorithms. Since it was found that a higher degree of irregularity in the HGS was associated with a least-cost design, the fractal dimension may be used as a parameter to maximize during the calculation of the optimal set of pipe diameters. Moreover, optimization algorithms may also take advantage of the fact that the fractal dimension of the OHGS is influenced by the fractal properties in the hydraulic behavior of the WDNs, as shown in Fig. 13, to search for suitable diameter configurations that induce a higher degree of irregularity or unevenness in the flow and energy distribution inside a given WDN. Second, the fractal analysis of the HGS may be used in the operational improvement of WDNs. Prospective work can be done to identify how the fractal behavior of the HGS of a given network should be to accomplish a better performance in terms of water quality, reliability, reduction of water leakage, and operational cost. In this case, the fractal dimension of the HGS can be used as a parameter to guide the adjustment of a given WDN via valve operation, asset management, and renewal planning to establish the ideal hydraulic behavior.

Finally, there is Table 3. The cumulative frequencies show a similar outcome to what was found in the case of the OHGS, where the least-cost design presented a higher fractal dimension in comparison with the nonoptimal HGSs. In terms of the flow and energy criterion, the results may imply that a greater degree of irregularity or unevenness in the distribution of the available discharge and energy is associated with an optimal design. In relation to the infrastructure criterion, one would suspect that the higher dispersion in diameter values for the nonoptimal designs, as established in Fig. 7, would translate into a greater fractal dimension; however, the opposite was found. This can be explained considering that the individual weight of a node, when the infrastructure criterion is used, is affected not only by the magnitude of the surrounding diameters but also by the number of pipes connected to the node, thereby inducing a spatial irregularity in the individual weights when nodes with a different number of connected pipes are present in the network. Nevertheless, it is very important to mention that, in many of the analyzed case studies, multiple designs presented the same fractal dimension. The standard deviation values, as presented in Table 3, support this statement considering that a very small dispersion was found in each data set. This situation may be attributed to two feasible causes. First, it is possible that the proposed box-covering algorithm may have difficulty detecting the difference between one design and another in terms of their fractal properties and their corresponding hydraulic behavior. Second, it is plausible that there are no noticeable differences between the fractal properties of the involved network designs whatsoever. This issue should be addressed in future work before applying these findings to the optimal design and operational improvement of WDNs.

In this paper, a fractal analysis methodology for the OHGS and the underlying structure of a WDN was presented. Multiple gravity-driven benchmark and Colombian networks with different topographical, topological, and energy availability conditions were designed so as to identify a least-capital-cost solution and several nonoptimal alternatives in accordance with the minimum pressure requirements and nodal flow demands. The OHGS and nonoptimal HGSs were successfully calculated for each network by conducting a numerical interpolation procedure, whereas the calculation of the corresponding fractal dimension was accomplished by the employment of the variation estimation method. Moreover, the fractal dimension of the underlying structure of each network design was estimated using an adapted greedy box-covering algorithm in which topological and hydraulic criteria were involved. Although the overall performance of the proposed fractal analysis algorithms was satisfactory, prospective work could be done to improve the accuracy and sensitivity of the chosen procedures. In addition, future work should also include the analysis of additional networks in which pumping devices are included.

The results show that the OHGS exhibited fractal properties that suggested a self-similar behavior in the energy dissipation pattern of least-cost design. In 64.29% of the case studies, the fractal dimension of the OHGS was greater than at least 75% of the fractal dimension values of the nonoptimal HGS; thus, the way in which the available energy is used in the least-cost design of these systems favors a steeper and sharper hydraulic gradient that translates to OHGSs with a rougher morphology. The behavior in the remaining 35.71% case studies, where the fractal dimension of the OHGS was only greater than 25%–50% of the fractal dimension values of the nonoptimal HGS, may be explained by the topological structure, nodal flow demand pattern, and the planelike morphology of the corresponding HGSs. However, prospective work should be done to identify additional details to explain the occurrence of these exceptions. In addition, the results also indicate a growing nonlinear relationship between the fractal dimension of the OHGS and the fractal dimension of the optimized networks for any of the four calculation criteria: topological, flow, energy, and infrastructure. Hence, the sparseness of the nodes and pipes, the heterogeneity in the diameter values, and the unevenness in the distribution of the available flow and energy in an optimized design affect the geometrical properties of the OHGS. Furthermore, the results show that the fractal dimension of the underlying structure of the optimized networks using the flow, energy, and infrastructure calculation criterion was greater than or equal to six or more fractal dimension values of the nonoptimal designs in 64.29%, 75%, and 57.14% of the case studies, respectively. This implies that the hydraulic behavior of a least-cost design exhibited, in most of the case studies, a higher degree of irregularity in the distribution of the available flow and energy and in the diameter values compared to a more expensive design. However, it is also necessary to consider that, in many networks, some designs exhibited an identical fractal dimension. This suggests either that the chosen algorithm may be incapable of identifying the differences between some designs or that there is no difference in the fractal behavior of the involved designs at all. Prospective work should address this issue.

The principal outcomes of this study suggest that the fractal analysis of the HGS and of the underlying structure of a given network may have useful applications in the context of WDN capital cost optimization and operational improvement. Since a least-cost WDN design exhibits a particular fractal behavior that is related to its inherent hydraulic properties, the fractal dimension may be used as an additional fitness measure to guide new or existing optimization algorithms. In relation to operational improvement applications, forthcoming research may use the fractal analysis as a tool to identify how a given WDN design should be operated, modified, or renewed in order to enhance the hydraulic performance of the network and reduce the operational cost.

Some or all data, models, and code generated or used during the study are available in a repository online in accordance with funder data retention policies. The available materials are EPANET data files, Excel data files, and Python scripts to generate figures and tables. The materials can be found in the following repository: https://doi.org/10.5281/zenodo.6634751.

Figs. 1–3, 6–12, and S1–S40 and Tables 2, 3, and S15–S19 can be reproduced by accessing the uploaded material in Jaramillo (2022). Arturo Rodriguez (Research Assistant, Master’s Student in Structural Engineering, Universidad de los Andes) ran all scripts and reproduced all tables and figures.

The authors would like to thank PAVCO Wavin for providing essential information for the development of the research work as well as for their constant support in investigation activities concerning water supply.