Shortest-Path-Based Two-Phase Design Model for Hydraulically Efficient Water Distribution Network: Preparing for Extreme Changes in Water Availability

Environmental issues can cause changes in source water availability in water distribution networks (WDNs). Thus, an efficient connection between the source and consumers is important for securing water serviceability, which can generally be achieved by minimizing energy losses. In this study, a novel two-phase design (TPD) model is proposed to design an energy-efficient WDN by maximizing a hydraulic geodesic index (HGI), which is the weighted shortest path from the source to the demand node. Before applying the TPD model for WDN design, a correlation analysis between the system HGI, hydraulic performance, and graph theory indices is conducted using 33 J-City networks to verify the proposed HGI. Next, the TPD model is used to determine the optimal layout of the grid network (Phase I). Based on this layout, the optimal diameter set is identified in Phase II. The TPD is thereafter compared with the traditional single-phase design (SPD) model, which determines the optimal layout and diameter simultaneously, and a least-cost model for each phase in the grid network layout and pipe-sizing problem. The correlation analysis clearly indicates that the system HGI with the weighted graph theory successfully determines the hydraulic performance without any hydraulic analysis. Furthermore, TPD is advantageous for designing energy-efficient, hydraulically and structurally sustainable, and resilient networks, as compared to SPD and the least-cost model. The TPD model is expected to provide a better opportunity to prepare for extreme water availability changes by enhancing the hydraulic performance and efficiency through a better connection between the source and nodes.


I. INTRODUCTION
A water distribution network (WDN) requires pipes, valves, and pumps to connect consumers with distant water sources. Fresh water is essential for humans; therefore, the connection should be maintained through proper WDN management. However, emerging environmental issues, such as climate change and extreme drought, may cause significant changes in the availability of water from water sources. This decreases the available volume and total head of water in the sources The associate editor coordinating the review of this manuscript and approving it for publication was Baoping Cai .
(generally the point with the maximum potential energy); thus, the resultant potential energy used to deliver according to customer demands may be less than that under normal conditions [1]. Under these circumstances, an efficient connection between the source and consumers is important to secure the serviceability of a WDN, which can generally be achieved by minimizing energy losses [1]- [3]. However, the determination of such a connection is challenging because of the complexity of the network connection and locations of assets (e.g., pipes, pumps, and valves) and the source (e.g., reservoir) [4]- [6]. Furthermore, as the WDN performance significantly depends on the network topology [4], [7], such topological complexity and dependence highlight the necessity of an innovative approach to system design and analysis [8].
The WDN topology can be simply regarded as a large planar graph [9], [10]; accordingly, graph theory (GT) has been introduced into WDN topology analysis [5]. GT explores the relationship between the structural properties and the eigenvalues and eigenvectors of the corresponding matrices [11]. GT has been widely applied in various areas, including data analysis [11]- [13], communication [14], [15], traffic networks [16], [17], and energy networks [18], [19]. Typically, GT represents a network in a mathematical graph as G = G(V, E), where V denotes vertices (e.g., demand nodes, reservoirs, and tanks in a WDN) with n elements, and E represents edges (e.g., pipes, pumps, and valves in a WDN) with m elements [5]. Each edge has a pair of vertices (i, j), where i = j; i, j ∈ V; and i and j are neighbors [9].
In GT, classifications can be made based on whether the edges are directed and/or the edges and vertices are weighted [3], [10], [36]. In previous studies, undirected and unweighted graphs were commonly applied. However, a few researchers have argued that a directed and weighted graph is more suitable for representing the WDN behavior according to the flow rate and direction of pipes as well as the elevation and demand of nodes [3], [10], [32], [36]. For example, if a network experiences extreme water availability changes due to a severe drought in the upstream of a service area (e.g., a 50% decrease in water volume delivered and stored in two reservoirs that supply the network), the hydraulic conditions within the network will change significantly. Changes in the distribution and direction of pipe flow rates affect the head losses and nodal pressures. Such transitions cannot be captured with undirected and unweighted graphs; thus, a directed and weighted graph should be adopted to reflect the changes and take into account accurate network connectivity. Such limitations of unweighted and undirected GTIs have been recognized by Jung and Kim [24], especially for WDN design problems.
A few examples of directed and weighted graphs have been introduced. Yazdani and Jeffrey [37], for example, adapted a directed weighted graph to be used for a WDN or to be applied to vulnerability analysis using the demand-adjusted degree of entropy. The water-flow closeness and K-shortest path introduced by Herrera et al. [32] are other examples of weighted graphs. These indices determine whether the connection of the demand node from the source is appropriate. Lee and Jung [3] proposed the concept of the source-to-node shortest path (SNSP) as an alternative measure to determine connectivity. In contrast to other measures, the SNSP approach deals with direct connections between the source and demand nodes and considers the efficiency of each connection.
Generally, WDN design studies focus on economic solutions while meeting additional criteria (e.g., minimum pressure, water quality, resilience, robustness, energy efficiency, etc.) [38]- [42]. However, most of the WDN design studies assumed that the layout of the WDN was fixed (or known) and determined the diameter of each pipe in the network. This assumption may provide a feasible solution but may fail to find better options by limiting the search space [43]. A few studies have focused on the WDN layout [44] or layout and diameter simultaneously [45]- [47]. However, these studies considered the WDN layout as an additional factor; therefore, they determined the layout and diameter simultaneously. Moreover, the effectiveness of the design has not been explored thoroughly, particularly from the perspective of energy efficiency. In addition, to the best of the authors' knowledge, the feasibility of weighted GTI to WDN layout and diameter design has not yet been discussed.
Therefore, this study introduces a two-phase design (TPD) model for WDNs using a hydraulic geodesic index (HGI), which was developed to minimize energy dissipation and to prepare for extreme water availability changes. The HGI is a new weighted GTI, which is the shortest weighted path from the source to the demand node. In other words, HGI can be used for WDS designing problem by overcoming the limitations of the unweighted GTI, as recognized by Jung and Kim [24]. The optimal layout is determined through Phase I of the TPD model using the system HGI (SHGI), while the model used in Phase II identifies the optimal diameter set for the optimized layout. The proposed model is first validated by conducting a correlation analysis between the proposed HGI, the hydraulic performance index (HPI), and the GTI using 33 real WDNs with varying characteristics in J-City, South Korea. The proposed TPD model demonstrates the optimal layout and diameter design of the grid network in B-City, South Korea. To identify the benefit of using the proposed TPD model for a WDN layout design, the optimal layout results of the model are compared with those of the layouts obtained using a least-cost model; thus, the total construction cost (TC) is minimized. In addition, the least-cost objective single-phase design (SPD) model, which simultaneously determines the optimal layout and diameter, is applied to the same grid network to investigate the performance of the proposed TPD model.

II. TERMINOLOGY
Several indices were considered in this study. To avoid confusion, the indices and their hierarchical relationships are sum-VOLUME 9, 2021 marized and described in Fig. 1. A GTI indicates a network's topological characteristics but does not consider hydraulic performance or other factors. It includes algebraic connectivity (AC), average node degree (AND), average path length (APL), link density (LD), meshedness coefficient (MC), and spectral gap (SG) [48]. In contrast, an HPI is quantified as a result of hydraulic analysis and includes indexes determining system resilience, such as the resilience index (RI) [49], network resilience index (NRI) [50], modified resilience index (MRI) [51], and energy efficiency (EE) [2]. The details of each indicator are described in the following sections. The proposed HGI, associated with the intersection between the GTI and HPI, applies the aspects of both indicators (Fig. 1). This indicates that HGI can take advantage of both indicators, specifically, the simplicity of calculation of GTI and the prediction of system hydraulics of HPI.
In GT, the term ''geodesic'' indicates the shortest path and is commonly employed to quantify GTIs, such as the network diameter or radius [33], [34]. The term ''hydraulic geodesic'' (HG) does not indicate the shortest physical Euclidian path, but uses the geodesic concept to represent the shortest hydraulic path (i.e., degree of energy dissipation between two nodes).

III. METHODOLOGY
In this study, a TPD model for a WDN layout and diameter optimization was developed by introducing a new GT-based HGI (Fig. 2). Details of the HGI and design model, as well as the background of the HPI and GTI applied in the present study are also described. For all optimizations, a revised harmony search (ReHS) was applied [52].

A. HYDRAULIC GEODESIC INDEX
When analyzing a WDN as a planar graph, it is important to focus on the connections between the sources and demand nodes rather than on the connectivity between any two arbitrary nodes [5], [53]. In this paper, the HG is proposed to consider the strength of the hydraulic connection when calculating the shortest pathway between the sources and demand nodes. Here, the hydraulic connection strength indicates a lower head loss (or energy loss) between the source and demand nodes. Although the head loss in a pipe is proportional to the physical Euclidean distance (i.e., pipe length), a longer water path length should be assigned to a node for which the upstream pipes have high head losses due to aging. To consider the hydraulic factors in determining the shortest path length, the weighted pipe length (w) for each pipe is first calculated from the Hazen-Williams head loss equation. However, because the HG does not involve hydraulic analysis (i.e., the flow rate is unknown), only the resistance coefficient of the head loss equation is considered as the weight, which is the product of the normalized roughness coefficient, normalized pipe diameter, and normalized pipe length.
where norms C i , D i , and L i are the normalized Hazen-Williams roughness coefficient (C factor), the normalized diameter, and the normalized length of the ith pipe, respectively. The maximum value of each parameter was used to normalize the weight between 0 and 1. Norm C i is normalized by the maximum C factor of a pipe, whereas norms D i and L i are normalized based on network-wide maximum values. The C factor can be considered as a static parameter [54] unless maintenance activities are employed to restore it. Estimating the C factor usually requires a large amount of data [54], [55], and an empirical equation based on only a few parameters is often useful for estimating the C factor [38], [51], [56]. In terms of practical applications, this simplification plays a major role in engineering judgment. The following empirical equation was used in this study [57]: where e 0 is the initial roughness (mm); a is the roughness growth rate (mm/year); and t is the number of years after installation; thus, C is maximum when t = 0. The pipe roughness varies over time at a rate that depends on the pipe materials, water quality, and pipe linings [57], [58]. In this study, e 0 and a are set to 0.18 mm and 0.4 mm/year, respectively [38]. The C factor calculated by the empirical equation generally provides a reliable estimation and supports the decision for engineering judgment; therefore, it can be used for SHGI calculation. Among the multiple pathways from the source to each node, HG is the pathway with the least weight summation (from (1)). It is calculated using Dijkstra's algorithm, which determines the shortest path between two nodes [59]. One of the two nodes is usually the source of the HG calculation. Although the HG can exceed 1, a higher value represents a weak connection (i.e., higher energy dissipation) with the source. Thus, the shortest pathway determined using Dijkstra's algorithm was the preferred pathway. For better intuitiveness (a higher value indicates a better result), the inverse of HG is proposed and defined as the HGI as follows: where HGI k is the HGI of the kth node; HG k is the HG of the kth node; and min HG is the minimum HG throughout the network. Here, min HG is applied to normalize all the values of HGI k between 0 and 1. It should be noted that HG is independently determined for an individual node; thus, the shortest paths are not necessarily identical if the nodes are adjacent to multiple shortest paths in the area. As the HGI is normalized using min HG, a higher HGI indicates a better connection from the source in terms of energy dissipation, and therefore a higher possibility of enduring extreme changes in the source. Finally, the average HGI k is defined as the system HGI (SHGI).
where Avg(·) is the average value. An energy-efficient network design can be obtained by maximizing the SHGI, as a higher SHGI implies that a system has a better connection (i.e., lower energy dissipation).

B. HGI-BASED TWO-PHASE DESIGN MODEL
In this study, a TPD model was developed considering the SHGI. The TPD model for Phase I (Phase I model) determines the optimal network layout of the WDN by maximizing SHGI (F 1 ). The objective function and constraints for optimization are as follows: where LO a is the network layout condition when all pipes are installed at all candidate locations; P min is the minimum pressure requirement; and P k is the pressure at the kth node. Note that candidate locations indicate all possible locations where pipes can be buried based on the locations of the demand node and their adjunct demand nodes. The decision variable for the Phase I model is the pipe status, that is, closed or open, and the hydraulic constraints (P min ) of the determined layout (based on pipe status) were checked through the hydraulic analysis program EPANET 2.0 [60]. Therefore, the Phase I model determines the status of all possible connections based on the location of each node in the WDN. From a practical standpoint, the closed pipe status indicates no pipe installation at the location.
The TPD model for Phase II (Phase II model) is a TCconstrained pipe-sizing model that seeks the optimal diameter set to maximize the SHGI (F 2 ), given the optimal layout derived from Phase 1 ( LO 1 ). Pareto optimal solutions were obtained by independently utilizing the proposed model for different TC values. Note that such a TC can be considered a limited budget for the design. The TC is calculated using a pipe cost model [61] that utilizes the parameter values determined by Jung et al. [62]. The objective function and constraint are expressed as follows: where Min(TC) is the minimum cost of the optimal layout determined from Phase I (LO 1 ); α is the cost multiplier; and TC x is the total cost of the set x of decision variables. Min(TC) is obtained through the pre-optimization of the pipe diameter of LO 1 with the objective of minimizing TC, and the model is defined as the least-cost pipe-sizing model. Thereafter, Min(TC) is utilized as the baseline for the Phase II model in maximizing the SHGI with a limited budget, which varies based on α.

C. HYDRAULIC PERFORMANCE INDICES
A GT analysis generally does not require a hydraulic analysis; this is particularly advantageous when a complex hydraulic analysis is involved. However, this also produces ambiguity in the hydraulic performance of a WDN and increases problems regarding possible unsatisfactory hydraulic requirements (e.g., minimum pressure and velocity) [24]. To explore the relationship between the proposed HGI and hydraulic performance, four indices (three hydraulic-analysis-based resilience indices, RI, and two other modified versions, i.e., the NRI, MRI, and EE) were applied. These HPIs are selected as they all consider the energy of the network while involving the concept of resilience. The resilience of WDNs has received considerable attention owing to uncertainties in disturbances and stochastic failures following such disturbances [63]. Therefore, improving resilience improves the preparedness of the network toward various disturbances. The RI is a fraction of the available energy surplus at the nodes over the maximum energy surplus in the network, which is internally dissipated to satisfy the demand and head required at the nodes. Prasad and Park [50] revised the RI to be used as a system NRI by adding a uniformity parameter to represent the loop system reliability. Jayaram and Srinivasan [51] also proposed an MRI based on Todini's 50 RI to extend its applicability to a WDN with multiple sources. Finally, the EE is calculated from the ratio of the energy delivered at the nodes to the energy supplied from the sources [2]. The EE is selected because it is based on the law of conservation of energy and is an indirect measure of the dissipated energy [2]. Therefore, when the energy dissipation between the source and nodes is low, a high EE value is obtained. The aforementioned four indices were utilized to evaluate the effectiveness of the proposed HGI in predicting changes in hydraulic performance. VOLUME 9, 2021

D. GRAPH THEORY INDICES
Numerous GTIs have been used to measure the topological characteristics of WDNs. The GTI used for a WDN can be classified as either statistical or spectral. Statistical measures quantify the organizational properties of a network based on the most frequent motifs and structural patterns and relate them to network robustness [10]. The spectral measures derived from the spectrum of the network adjacency matrix quantify the network invariants. These invariants, when considered along with the described statistical measurements, reveal useful information regarding the quality of connectedness, connectivity strength, and failure tolerance of the network (link and node connectivities) [10].
Six unweighted GTIs were selected to investigate the novelty of the proposed HGI: AC, AND, APL, LD, MC, and SG. AC represents the strength among the connections in the network, and a larger AC implies a stronger connection. AND quantifies the average number of connections per node and provides immediate information regarding network organization; a larger AND implies a more organized network. APL is similar to HG but differs by considering possible network pathways among all nodes rather than from the source to nodes; a smaller APL usually indicates a more robust network. LD provides information regarding the general connection among graph nodes in terms of ''inclusivity'' [36]; a larger LD is preferable [9]. The MC evaluates the AC in measuring the fraction between the actual and maximum possible numbers of network loops. A larger MC indicates that there are more loops in a network, which is interpreted as being more redundant [5]. SG measures the strength of network connectivity and provides valuable information regarding bottlenecks, articulation points, or bridges; a larger SG indicates a more robust network [64]. The general trends of each GTI are utilized to investigate how the HGI develops the topological characteristics of a WDN when employed in design.

E. REVISED HARMONY SEARCH
For all optimizations, the ReHS was applied [52]. ReHS is different from the original HS [65] because of changes in the harmony memory considering rate (HMCR) and pitch adjusting rate (PAR) across the iterations used to dynamically balance between exploration and exploitation, and changes to the number of solutions in the harmony memory (HM) to be replaced by the HMCR and PAR values. During the earlier iterations, high HMCR and low PAR values increase the search speed in finding the global optimum. As the iterations progress, the HMCR decreases, whereas the PAR increases, avoiding local optima and more quickly finding the global optimum. The HMCR and PAR for each iteration are updated as follows: where HMCR iter is the HMCR at iterations; HMCR max and HMCR min are the maximum and minimum HMCR; iter max is the maximum number of iterations; PAR iter is the PAR at iterations; and PAR max and PAR min are the maximum and minimum PAR, respectively. For each iteration, the amount of HM replaced by the HMCR and PAR can be calculated as follows: where N HMCR,iter is the number of new HMs in which the HMCR will apply; N PAR,iter is the number of new HMs in which the PAR will apply; and HMS is the size of the HM. In this study, an HM of 100, HMCR max and HMCR min of 0.95 and 0.3, respectively, and PAR max and PAR min of 0.03 and 0.3, respectively, are considered for all optimizations.

IV. STUDY NETWORK
The proposed TPD model was first applied to 33 WDNs in J-City (J-City networks), South Korea, to validate the HGI by conducting a correlation analysis between HGI and HPI or GTI. Then, a TPD model was demonstrated with the grid network in B-City, South Korea. Fig. 3 illustrates five representative J-City networks and a grid network with all possible connections. The TPD model was built on the Python platform, the SciPy package was used to calculate the SHGI and all other GTIs [66], and the HPI was evaluated based on the EPANET simulation. The following subsections describe the network details and case studies.

A. J-CITY NETWORKS
The TPD model was first applied to 33 J-City networks, which have a near-grid structure (Fig. 3). The number of nodes varied from 14 to 181, whereas the number of pipes was between 12 and 124. The network length is approximately 8500 km on average, and the average network demand is approximately 27.7 lps. All networks have a minimum pressure head requirement of 15 m (= 21 psi), whereas a single source supplies water to the network.
Pearson correlation analyses were conducted to check the level of linear relationship between the proposed SHGI and other indicators (Table 1), that is, between the SHGI and HPI (PCC#1), between the HGI and node hydraulic parameters (pressure or pressure head, PCC#2), between the SHGI and GTI (PCC#3), and between the SHGI and network topological characteristics (e.g., the number of links and nodes, and the total length of links; PCC#4). The impact of pipe aging on PCC#1 was investigated by increasing t from 0 to 50 years in (2) and decreasing the C factor in (1), and PCC#1, PCC#3, and PCC#4 were quantified under two different weighting conditions (i.e., weighted with w i = 1 and unweighted with w i = 1).

B. GRID NETWORK
The developed TPD model was demonstrated using a grid network (Fig. 3(f)) [67]. The network comprises a maximum of 61 pipes, (each 2000-m long) with a total of 36 nodes and is gravity-fed from a single source with a total head of 80 m. The demand nodes (all at the same elevation) are located within a 6 × 6 grid topology, and each demand node requests 94.7 lps, with a total network demand of 3409 lps. The C factor was 120 for all pipes, and the minimum pressure requirement was 28 m (= 40 psi).
First, the Phase I model was applied to identify the optimal layout of the grid network, which was compared to the network identified using the least-cost layout model. Although the least-cost layout model was applied to minimize the TC, the constraints and other optimization conditions were identical for both models. A uniform pipe diameter of 2000 mm was used to minimize the impact of hydraulic performance due to a small diameter; accordingly, only the optimal pipe locations were determined through Phase I optimizations. The resulting layouts were then compared to their topologies, that is, the SHGI, GTI, and performance HGI. A convergence test was conducted to select a reasonable number of optimal layout generations because slightly different layouts and indicator values can be obtained from each phase I optimization run. We confirmed that 1000 independent phase I optimizations were required to provide a sufficiently converged indicator value.
Next, the Phase II model determines the optimal diameter set based on the layout determined from Phase I. The proposed TPD model runs Phases I and II as a serial application. Sixteen commercial pipe diameters (i.e., 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1200, 1400, 1600, 1800, and 2000 mm) were considered as candidates. Independent Phase II optimizations were performed by varying the cost multiplier, α [see (6b)] from 1.2 to 2.0 at 0.2 increments. Again, these indicate the budget limit for WDN design and help to identify the marginal cost of the WDN design. The TPD model results were compared with those of the SPD model, which simultaneously determined the layout and pipe diameters in a single optimization run by minimizing the TC. The SPD model represents the traditional approach for determining the optimal layout and diameter of the WDN. The SPD model also uses ReHS for optimization, as described in Section III. The resulting optimal design solutions were also compared with respect to their hydraulic performance under fire flow and different demand conditions. We assumed that a fire flow of 63.09 lps (= 1000 gpm) occurs at a demand node located in the northeast corner of the network [dashed circled node in Fig. 3(f)]. Four different demand scenarios were used for the fire flow test (uniform demand multipliers were 1.0, 1.2, 1.5, and 1.8), and the minimum pressure head VOLUME 9, 2021 requirement of 14 m (= 20 psi) was applied in the result analysis according to Utah Administrative Codes R309-510-9 [68] and Utah Administrative Codes R309-105-9 [69]. The advantage of applying the TPD model for a WDN layout and diameter design was verified by examining the marginal cost of the design, pipe diameter, and flow distribution.  Table 2 summarizes several statistics (i.e., maximum, minimum, average, and standard deviations) of PCCs obtained from 33 J-City networks (PCC#1). For the PCC calculation, each pipe in each network was closed to simulate single-pipe failure conditions, and the corresponding HPI and SHGI results were collected. The HPI exhibits an average PCC value of more than 0.9, and a minimum as low as 0.61 for the MRI. In particular, it can be confirmed that the SHGI is significantly related to the EE. A comparison of the weighted SHGI to the unweighted SHGI (values in parentheses in Table 2 ) showed that the weighted approach has an improved correlation with other HPIs compared to the unweighted approach.   4 shows the box-and-whisker plots drawn using the correlation analysis results of the network with new pipes and 50-year-old pipes (identified by 0 and 50 in the figure, respectively) to determine the impact of pipe aging (decreasing C factor) on the PCC between the SHGI and HPI. As the pipes age, the PCC variance increases, while the average PCC decreases.

A. CORRELATION ANALYSIS
The correlation between the HGI and the two nodal hydraulic parameters (i.e., pressure and total head) was identified using PCC#2. The correlation between the two parameters (PCC values of 0.08 and 0.35 for nodal pressure and total head, respectively; not shown in the tables and figures) was found to be low because the HGI was not quantified based on hydraulic simulation; it was quantified according to the topological connection and weight [see (1)] and did not consider the weight based on the nodal element (e.g., elevation, demand, etc.). For a similar reason, a low PCC was found in a network with a high spatial variation in the nodal hydraulic characteristics in PCC#1.  Table 3 summarizes the correlation analysis results of PCC#3 (between the SHGI and six GTIs) and PCC#4 (between the SHGI and eight network topological characteristics). The PCC#3 results indicate that the weighted SHGI has a low correlation with the GTI, whereas the unweighted SHGI has a high correlation (either positive or negative) with most of the GTIs. This is because the GTI is based on an unweighted graph; thus, the similarity between the unweighted SHGI and GTI was higher. For a similar reason, the unweighted SHGI has a higher correlation with the network topological characteristics than the weighted SHGI in PCC#4. SHGI is a new indicator based on GT; these results highlight the novelty of SHGI by exhibiting its differences from other GTIs.
The multivariate linear regression model using topological characteristics as independent variables and SHGI as a dependent variable yielded an R 2 value of 0.735 and a p-value less than 0.05; the total length of the links, average diameter, and average elevation were considered as the parameters.
The PCC results indicate that the proposed SHGI is more similar to the HPI than to the GTI, although its calculation is considerably more similar to that of the latter than that of the former (i.e., path length calculation based on network connectivity), except when considering weight (1). An energy-efficient and resilient network can be constructed by maximizing the SHGI, which is highly correlated with EE and system resilience. It is therefore concluded that the proposed SHGI suitably reflects the hydraulic network performance even without a time-consuming hydraulic analysis and may be useful in actual practice, particularly when the network experiences extreme water availability changes.

B. PHASE 1: LAYOUT DESIGN
An optimal layout of a grid network was first identified using the Phase I model, and then compared to the layouts obtained using.

Figs. 5(a) and (b)
show the layouts obtained from the least-cost model with the minimum (TC-minSHGI) and maximum (TC-maxSHGI) SHGI values, respectively. The most inexpensive solutions tend to minimize the number of pipes and create a branched network with a single long water supply line. All 1000 least-cost layout solutions have 36 pipes (the minimum number required to connect all demand nodes). Thus, all the least-cost layout solutions have the same TC regardless of the layout, while the SHGI value is different for each solution, that is, a maximum of 0.21 [TC-maxSHGI, Fig. 5(b)], a minimum of 0.13 [TC-minSHGI, Fig. 5(a)], and an average of 0.17. Fig. 5(c) shows a representative optimal layout (SHGIlooped layout) determined using the Phase I model; the HPI, and GTI are approximated averages, as listed in Table 4. In contrast to the least-cost layout, the Phase I model produces a looped network with 46.95 pipes on an average [minimum and maximum of 40 and 55, respectively; 47 pipes are shown in Fig. 5(c)]. The construction of a looped network has an advantage over a branched network with respect to HGI and SHGI. While a node-optimal short water supply path can be developed in the former, few overlapping paths should be shared between the upstream and downstream nodes in the latter, which results in greater energy dissipation at higher flow rates. The phase I model that is used to maximize the SHGI therefore yields a dense looped layout with the highest  For a fair comparison, the nodal HGI values in the SHGIbranched with 36 pipes were compared with TC-minSHGI, and TC-maxSHGI to calculate their differences (the former minus the latter), as shown in Figs. 6(a) and (b), respectively. As expected, considerable differences between SHGI-branched and TC-minSHGI were observed. Note that discrepancies in the HGI between the SHGI-branched and TC-maxSHGI [ Fig. 6(b)] values only exist at the nodes downstream of detour pipes (flow is supplied in an energy-inefficient manner from the source-opposite to the source-toward sides), as indicated by the dashed circle in Fig. 5(b). This indicates that the Phase I model (i.e., SHGIbased) creates a layout with a more efficient flow path than the one from the least-cost objective by eliminating such detours in the network. Based on the list in Table 4, it can be confirmed that, on average, the Phase I model produces networks with higher HPIs and GTIs compared to the least-cost model. Maximizing the SHGI increases these two indicators and TC. The average TC of layouts (in the form of an SHGIlooped) determined from the Phase I model is approximately 1.3 times higher than that of layouts from the least-cost model because more pipes are installed (a uniform pipe size is used in both models; see Table 4). Interestingly, the average TC of the SHGI-branched is the same as that of the least-cost model layouts because both have the same number of pipes. VOLUME 9, 2021 However, the SHGI of the SHGI-branched is higher than that of the least-cost model layout, leading to higher HPI and GTI values in similar ranges. The results of the SHGI-branched summarized in Table 4 prove that the optimization of Phase I with a cost constraint can result in maximum SHGI and HPI values with a limited TC.

C. PHASE 2: PIPE SIZE DESIGN
The proposed TPD model was applied to determine the optimal pipe diameter of the SHGI-looped [ Fig. 5(c)] and SHGI-branched [ Fig. 5(d)] layouts using the Phase II model. Before applying the Phase II model, the least-cost pipe-sizing model was used for both layouts to identify the baseline cost [Min(TC) using (6)] of the Phase II model. The baseline costs for the SHGI-looped and SHGI-branched layouts were US$ 50.40, and the distribution of the least-cost SHGI-looped design was compared with that of the SPD model, which determines both layout and diameter simultaneously with the minimization of the TC as its single objective. The least-cost SHGI-looped design produces a consistently higher pressure, regardless of the demand scenarios considered, compared to the SPD design. Independent optimizations were conducted by varying α in (6b), from which the Pareto fronts of the SHGI-looped and SHGI-branched designs (using the TPD model) between the TC and SHGI were derived and compared to those of the SPD model (Fig. 7). The least-cost solutions of the TPD model  Fig. 7. In contrast, the Pareto front of the SHGI-branched layout exhibits a weak nonlinearity; the increase in α resulting from the SPD optimization does not increase the SHGI, leading to a flat Pareto front. The maximum SHGI (0.368), which is considerably lower than that of the SHGI-looped layout, was obtained by the SHGI-branched layout at approximately α = 1.6. The SPD model did not show any relationship between SHGI and TC and yielded the lowest SHGI, regardless of the TC.  and TC = US$ 60 million; these values return the best marginal cost. The TPD model adopts two approaches to maximize the SHGI: 1) installing a commercial pipe with a maximum diameter (2000 mm) for the source; and 2) obtaining the smallest possible HG value (shortest weighted path length) at the nodes. The SHGI is the average value of the nodal HGI, which is an inverse of the nodal HG divided by min HG in the network. The min HG value is always calculated at the right downstream node of the source [ Fig. 3(f)]. For the first approach, the length and roughness factor are fixed [see (1)], and the diameter of the source pipe should be maximized to reduce the minimum nodal HG. The maximum diameter in the TPD solutions increased from 1200 to 2000 mm as α increased from 1.0 to 2.0. For the second approach, the pipe diameters were optimized to obtain the smallest HG value possible at the downstream nodes in a tree-like layout determined from Phase I [while satisfying the pressure constraint in (6c)]. Therefore, the resulting design exhibits a smooth decrease in pipe diameter from upstream to downstream. Comparing the diameter distribution , we observe that a greater pipe flow is supplied through the large pipes that pass through the middle of the upper and lower parts of the network.
As expected, the least-cost SPD model tends to reduce the number and diameter of pipes in the network. The layouts obtained are similar to those of TC-minSHGI and TC-maxSHGI (Figs. 5(a) and (b), respectively), and no SHGI increase can be observed in the SPD model solutions (Fig. 7) when a considerable investment amount is available. Moreover, we found that the HPI of SPD solutions is similar to that of SHGI-looped solutions with a low α. Fig. 9 illustrates the changes in the HPI of an SHGI-looped layout in response to different α values. Accordingly, it is demonstrated that the TPD model is superior to the SPD model in generating resilient and energy-efficient WDN designs to prepare for extreme changes in water availability.

VI. CONCLUSION
In this study, we proposed a TPD model that employs the HGI concept. The proposed model was first applied to J-City networks to validate the HGI method, and a layout and pipe diameter design for a grid network was demonstrated. The average value of the HGI is defined as the SHGI, and the TPD model considers maximizing the SHGI as an objective function. The TPD model was equipped with ReHS as the optimization technique, and the advantages of the proposed model and SHGI were explored by comparing their responses to the HPI and GTI, which reflect the hydraulic performance, particularly network energy efficiency and resilience to extreme changes in water availability.
The analysis results show that the SHGI has a higher correlation with the HPI than the GTI, particularly with the EE, despite the SHGI being a GT-based index. Moreover, the SHGI weight can reflect the hydraulic performance without the necessity of conducting a hydraulic analysis. The correlation analysis indicates that the proposed SHGI can be a useful measure of the HPI, leading to a more energy efficient and resilient network while taking advantage of the GTI, which does not require any hydraulic analysis.
The application of the TPD model to grid network design results indicates that the TPD model can aid in enhancing the EE and resilience of a network within a reasonable TC range. In addition, unlike to other GTIs, such process does not require additional hydraulic analysis as HGI is a weighted GTI with pipe characteristics in consideration. Based on the TPD design results, two recommendations can be provided to improve the SHGI. First, to create a higher SHGI network layout, the creation of a loop and elimination of detour flow pathways must be considered. Second, when determining the diameter, locating a larger diameter near the source and a smaller diameter at a distance from the source will improve the SHGI. These recommendations can support decision-making for WDN design to prepare for extreme changes in water availability by enhancing hydraulic performance and efficiency through a better connection between the source and node.
In summary, the proposed TPD model and HGI showed a clear benefit for energy-efficient WDN design. Although only a design example is presented here, this design can be applied to various academic and industrial purposes, especially for the existing WDN, including but not limited to maintenance prioritization, identification of critical pipes, and optimal operation. Maintenance prioritization can be achieved by comparing the SHGI improvement before and after the rehabilitation or replacement of each pipe. In contrast, the critical pipe can be identified by comparing the SHGI before and after pipe breakage (i.e., closing pipe), which creates a new topology of the network. Lastly, from the perspective of optimal operation, manipulating valves and pumps would render a network topology with improved energy efficiency.
Although the novelty and usefulness of the TPD model and HGI were considered in this study, there are some limitations and gaps that need to be addressed in future studies. First, the transition of the shortest paths and SHGI values can be investigated under different conditions of water availability and source head to identify the critical water level that triggers severe water-distribution failure. This investigation can also be conducted for a set of networks with different configurations and characteristics. Second, additional weight factors based on nodal parameters (node elevation and demand and tank characteristics) can be considered to modify the proposed index. As described earlier, the correlation between the SHGI and other indices is low when the node characteristics perform a critical function. The SHGI could better capture the changing nodal and source conditions (e.g., decrease in reservoir head due to drought) by incorporating node characteristics into the index. Third, as nodal parameters and dynamic operational components (pumps and valves) in a WDN have yet to be considered using this method, the use of the HGI should be extended under unsteady conditions to account for such state changes. Furthermore, a rigorous analysis should be implemented to guide the use of the HGI for decision-making in engineering. Finally, networks with varying configurations should be examined to confirm the conclusions of this study, and the limitations described herein should be overcome to increase the usability of both the TPD model and the HGI.