Pipe Break Rate Assessment While Considering Physical and Operational Factors: A Methodology based on Global Positioning System and Data-Driven Techniques

An accurate prediction of pipes failure rate plays a substantial role in the management of Water Distribution Networks (WDNs). In this study, a field study was conducted to register pipes break and relevant causes in the WDN of Yazd City, Iran. In this way, 851 water pipes were incepted and localized by the Global Positioning System (GPS) apparatus. Then, 1033 failure cases were reported in the eight zones of understudy WDN during March-December 2014. Pipes break rate (BRP) was calculated using the depth of pipe installation (hP), number of failures (NP), the pressure of water pipes in operation (P), and age of pipe (AP). After completing a pipe break database, robust Artificial Intelligence models, namely Multivariate Adaptive Regression Spline (MARS), Gene-Expression Programming (GEP), and M5 Model Tree were employed to extract precise formulation for the pipes break rate estimation. Results of the proposed relationships demonstrated that the MARS model with Coefficient of Correlation (R) of 0.981 and Root Mean Square Error (RMSE) of 0.544 provided more satisfying efficiency than the M5 model (R = 0.888 and RMSE = 1.096). Furthermore, statistical results indicated that MARS and GEP models had comparatively at the same accuracy level. Explicit equations by Artificial Intelligence (AI) models were satisfactorily comparable with those obtained by literature review in terms of various conditions: physical, operational, and environmental factors and complexity of AI models. Through a probabilistic framework for the pipes break rate, the results of first-order reliability analysis indicated that the MARS technique had a highly satisfying performance when MARS-extracted-equation was assigned as a limit state function.


Introduction
The deterioration of pipes causing pipe failures and leaks in urban Water Distribution Networks (WDN) has become the cornerstone of water utilities throughout the world. Pipe failures and leaks occasionally take place on account of a reduction in the watertransmitting capacity associated with the pipes and water pollution in the WDNs. Water utilities are generally exposed to a large amount of costs for the replacement and rehabilitation of water mains and consequently, this issue makes it critical to assess the current and forthcoming conditions of the WDN for maintenance decision-making (Berardi et al. 2008;Berardi and Guistolisti 2021;Moslehi and Jalili-Ghazizadeh 2020;Clair and Sinha 2012;Shi et al. 2013;Shende and Chau 2019;Stephens et al. 2020). In the case of economic and social views, expenses on the pipe breaks have a significant upward trend and as a result, utility managers are put under pressure to operate annually replacement plans for deteriorated pipes which balance investment with expected benefits in risk management issues. In this way, there is a ferocious demand for obtaining reliable pipe breaks techniques for evaluating the performance of WDNs (Berardi et al. 2008).
There are a variety of techniques to evaluate the deterioration and performance of water pipes which are classified into four groups: deterministic, probabilistic, and soft computing techniques. Generally, there is no denying the fact that influential factors and output results related to the pipe deterioration techniques are significantly inextricably bound up with the applicability of these methodologies (Clair and Sinha 2012;Shi et al. 2013). In the case of the deterministic technique, most of the available predictive models were implemented by using regression analysis, composed of mechanisticempirical analysis. The major drawback associated with deterministic models is that, in the case of practical use, these models are largely limited to a particular case study. The statistical model is intrinsically restricted while using assets with a piece of insufficient previously recorded information about WDN. As a major demerit, the regressionbased techniques are not appropriate for modeling the real-world deterioration process of pipe infrastructure due to different restrictions of the sampling data (e.g., pipe structure, loading conditions, and environmental variables). In probabilistic modeling, the relative frequency of an event is considered which is applicable to determine the failure of a segment of WDN. Probabilistic approaches estimate the distribution and range of dependent parameters and additionally, results of probabilistic techniques provide priority process repair, rehabilitation, and replacement related to infrastructures. According to the above-mentioned discussions, probabilistic and deterministic approaches are grouped into statistical techniques, which are based on the mechanically-based techniques, require a long-time period to obtain pipe failure, and additionally may be costly in computations. Moreover, these techniques require knowledge extraction of experts and this causes that these models are not always cost-effective (Berardi et al. 2008;Berardi and Giustolisti 2021;Moslehi and Jalili-Ghazizadeh 2020;Xu et al 2011b).
Ultimately, soft computing techniques especially a variety of Artificial Intelligence (AI) models, were generally flexible to combine with Evolutionary Algorithms (EAs) for various purposes such as optimizing time scheduling, costs on rehabilitation, replacement, repair, pressure fluctuations, and break rates (Berardi et al. 2008;Clair and Sinha 2012;Clark et al. 1982;Kettler and Goulter 1985;Jowitt and Xu 1993;Robles-Velasco et al. 2020;Shirzad et al. 2017;Silinis and Franks 2007;Tang et al. 2019). A survey on literature review proved that soft computing models had prosperous performance for evaluation of pipe break rates because such techniques could detect well-complicated patterns of pipe break processes. Compared to the statistical models, AI models have three main capabilities over traditional techniques for pipes break rate assessment: (i) resulting in non-linear explicit equations when AI techniques are fed by a high volume of data series, (ii) selecting independently effective factors (e.g., diameter of pipe, age of pipe, material type, and loading conditions) for the evaluation of pipe failure, and (iii) efficient generalization of AI performance for unseen records of the pipe break rate. In the case of pipe failure prediction, one of the most remarkable practical viewpoints of AI models' performance is to prognosticate values of pipe break rate for other regions' WDN whose failure records are not applied for developing AI models. This property of AI models is indicative of being cost-effective and a well-guarded design of WDNs.
In this study, three robust AI models, introduced as Gene-expression Programming (GEP), Model Tree (MT), Multivariate Adaptive Regression Spline (MARS), are employed to model pipe break rates associated with urban WDN of Yazd city, Iran. In this way, a map of WDN is processed in Geographical Information System (GIS) environment, and additionally, localization of pipes is surveyed using Global Positioning System (GPS) to complete break rates databases for the case study through a comprehensive field investigation. The results of AI models are statistically quantified, and additionally, reliability evaluation is studied.

Literature Review: Overview of Existing Techniques
Since 1980, a broad range of attempts have been made to study the break rate under various factors such as the age of pipe, various materials of pipe, size of pipe, number of failure, loading conditions, and underground environment (Clark et al. 1982;Jowitt and X 1993;Rezaie et al. 2015;Shi et al. 2013;Stephens et al. 2020). Jowitt and Xu (1993) proposed an applicable technique for evaluating the effect of different pipe breaks conditions on WDNs. Savic et al. (2006) assessed the time-dependent break rate of sewer systems with the help of a robust soft computing tool, introduced as Evolutionary Polynomial Regression (EPR). They presented an evolutionary algorithm-based relationship for pipe breaks which was successful in engineering interpretation. Wang et al. (2009) employed several multivariate regression techniques to prognosticate the break rate of pipes associated with large WDSs located in three municipalities (Moncton, Laval, and Quebec), Canada. From their research, they found that results of regression models helped efficiently who are experts in academic level, municipal engineers, consultants, in better understanding trends related to the break rate of water mains. Xu et al. (2011a) applied Genetic Programming (GP) to successfully predict the break rate of WDN in Beijing, China during a 19-year period (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) in comparison with the statistical model. Aydogdu and Firat (2015) employed fuzzy clustering into topology design of Least Squares Support Vector Machine (LS-SVM) technique to estimate break rate in Malatya WDN during a six-year period beginning in 2006, Turkey. Ultimately, they found that the results of LS-SVM were more accurate than those obtained by Feed Forward Neural Network (FFNN) and Generalized Regression Neural Network (GRNN) techniques. Wols et al. (2018) used a mechanical technique, introduced as Comsima, to compute the break rate with consideration of stresses and joint rotations in a part of WDN in the Netherlands based on quite a few loading conditions on the pipe (i.e., soil, traffic, water pressure, and differential settlements). Robles-Velasco et al. (2020) applied Logistic Regression (LR) and Support Vector Classification (SVC) to evaluate the pipe breaks in a part of WDN, Spine. The results indicated that the number of unexpected breaks is likely to be largely declined. Approximately 30% of pipe break rates could have been plummeted by substituting merely three percent of the incepted pipes annually, which is a realistic alternative. Ultimately, Yazdekhasti et al. (2020) employed efficiently four Machine Learning (ML) methods, as Classification Tree (CT), Random Forest (RF), LR, and, SVM to predict the water pipe breaks in the mid-Atlantic region during 1999-2018.

Overview of Case Study
The case study of the water distribution network was located in Yazd City, central Iran. This city located from 31̊ 53´50″ N to 54̊ 22´4″ E where has 529,673 population. The first region of Yazd's WDN, which is divided into eight zones (Z1-Z8) as Azadshar, Azadegan, Esteghlal, Emamshahr, Enghelab, Shahedieh, Saber Yazdi, and Jomhouri were considered as the case study. Fig. S1 (Supplementary Materials) illustrates eight sub-regions of WDN on a GIS-derived map. Retrieving sufficient pieces of information about the length of pipe was not possible on the GIS-based map since occurrences of pipes break in the WDN had been recorded on the basis where failure has been observed.
On the other hand, there was no dedicated identification code for each pipe. In this way, to complete a dataset of pipe breaks, a handy GPS whose software was installed on the cellphone was employed to find unknown coordinates in the WDN case study. The GPS generally has five merits over other traditional land surveying techniques: (i) a higher level of accuracy than traditional land surveying methods, (ii) computations are performed highly quickly and with a high degree of precision, (iii) visibility among stations does not impose a limitation on the operation of GPS technology, (iv) it can be conveniently carried to collect accurate data, and (v) quite a few GPS systems are capable of communicating wireless for delivery of real-time information. Figure S2 demonstrates the adaptation of GPS with the cellphone. Afterward, the localization of pipes that have been deteriorated, was known. The operation of detecting coordinates has been carried out within 45-day surveying with GPS. To verify GPSderived coordinates, a handy GPS manufactured by Garmin International Company was used. More specifically, the number of coordinates detected by GPS software on the cellphone was at same as coordinates retrieved by a handy GPS. This indicated that the results of GPS software were reliable for continuing the field study.
Water pressure fluctuation is regulated by gauges installed in the WDN of Yazd city, in which water pressure associated with the first region of WDN varies 1.5-2 bar. Moreover, materials of water pipes were made of Asbestos (As) (with nominal diameters: 80,100, and 150 mm), Cast Iron (CI) (with nominal diameters: 63, 90, and 110 mm), and Polyethylene (Po) (with nominal diameters: 63, 90, and 110 mm) which various nominal diameters were dedicated to each typical material of pipes. In recent years, 4000 pipes breaks of Yazd's WDN have been accumulated from Yazd Water and Wastewater Company (YWWC) which almost 2000 events were considered for analyzing the break rates of water pipe. Generally, in the case study of WDN, 675 pipes with 92.33 km-long which were damaged, were incepted to detect causes of failure and number of failures. By observing the pipes damaged in the WDN, 827 break cases have been recorded during a year period whose allowable pressure for operation of the case study WDN was 1.5, 1.7, 1.8, and 2 bar. Age of pipes for which was in the state of damage was computed between 1 and 10 years (yr). Moreover, the depth of pipes installation was 1, 1.1, and 1.2 m. According to the breakage reports of YWWC, there is a wide range of factors (i.e., leak, poor materials, excavation operation, pressure fluctuation, settlement, rusty water pipe, heavy vehicles, and inappropriate operation of water pipe) that lead to the pipe failure of Yazd's WDN. As reported from the YWWC, since the pressure fluctuations have remained relatively unchanged during 2015-2020, the equations returned by data-driven models (i.e., M5MT, MARS, and GEP) can be efficiently applied for the future as well as for loading conditions. From 2015 to 2020, an extensive repair of Yazd's WDN has been occurred, leading to a reduction of failure records for the first region of Yazd's WDN. In the present field investigation, typical breakages of water pipes were illustrated in Fig. S3.

Effective Factor on the Pipe Break Rate
According to the previous experimental and field investigations, effective factors which play a key role in the pipe failure processes, are generally categorized into three sections: physical, operational, and environmental factors. Physical factors are related to the properties of water pipe: age, material, length, thickness, and depth of pipe installation. The operational factor is associated with hydraulic conditions of water pipe: number of failures, mean pressure, water velocity, failure type, and pressure fluctuation or time series since the last breakage. The last factor is due to the environment such as traffic, soil type, soil corrosion, temperature, rainfall, soil resistivity, soil shrinkage swell, and freezing index (Berardi et al. 2008;Robels-Velasco et al. 2020;Savic et al 2006;Seifollahi et al. 2013;Xu et al., 2011a, b;Yazdekhasti et al. 2020).
In this study, the selection of contributory variables depends on the data availability of WDNs and the environmental conditions of the case study. All environmental factors, as mentioned above, were excluded from the list of variables affecting the failure process since the last records of breakage in the Yazd's WDN were not for the reason of traffic flows, chemical properties of soil, and freezing. In accordance with reliable breakage records, physical and operational factors can be considered into account for the evaluation of the failure process. Due to limited accessibility of WDN information, length of pipe (L P ), number of breaks (N B ), depth of pipe installation (h P ), allowable pressure in the operation status (P A ), and age of pipe (A P ) were selected as effective variables on the pipe break rate (BR P ). In this study, the values of BR P associated with inception operation were computed as, For the case study of WDN, BR P values vary between 0.024 and 77.70. More specifically, information about the number of failure occurrences, date of the event, type of pipe, and nominal pipe diameter was provided in Table S1 (Supplementary Materials) for the case study of WDN during March 2014-December 2014. To develop soft computing models for the prediction of BR P values, 675 records were considered in which 80 percent of observations were dedicated to the training, and additionally, the rest of (1) the observations were assigned for testing performance. Furthermore, to develop the practicability of the AI models, 176 observations related to the Jomhouri zone (Z8) are used.

Implementation of Model Tree
The Model Tree, as one of the widely practicable Data-Mining Techniques, has the high potential of approximating the performance of function (Quinlan 1992). Development of MT is completed within two steps to building a tree-like configuration. The first step is dedicated to dividing the search space of input variables into sub-divisions. In the second phase, the trees are created by means of datasets trapped in each subdivision, and a linear regression model is thereafter fitted on each subdivision datasets. Entire search spaces related to the linear models are expressed as a set of If-Then rules. The tree includes a collection of leaves and nodes. Each leaf generally describes linear regression formulation while nodes are representative of If-Then rules associated with input search spaces. The M5 model is one of the most efficient kinds of MT, which is used in this study. This linear model might have been simplified by removing many input variables that have a marginal contribution to minimizing the value of predicted error. In MT development, applying multivariate linear regression technique, as a curve-fitting technique in each subspace, is more robust than non-linear regression analysis. This assumption has been frequently proved by various applications of MT (Sattari et al. 2020;Shamshirband et al. 2020). In the M5 technique, Greedy Search (GS) is employed to eliminate excessive input vectors, in some cases, all input vectors are neglected and as a result constant values only remain for nodes (Keshtkar and Kisi 2018). In the Model Tree, Weka software has been used to model pipe break rates. In the M5 technique, the possibility of the overfitting problem has been efficiently controlled by using two operators (pruning and smoothing) tree structure that was generated in the training stage.
In this study, pipe length, pipe diameter, and the age of pipe were assigned as input variables and then pipe break rate is expressed by M5 as, where C 0 , C 1 , C 2 , C 3 , and C 4 are weighting coefficients of linear regression-based-equation by M5. The performance of M5, implemented by Weka software, provides 17 linear models along with if-then rules. Full descriptions of M5-driven equations were presented in Table S2. As seen in Table S2, L P is the first splitting input variable with a value of 0.0235 km. More specifically, if L P is greater than (or equal to) 0.0235 km, then A P is the sole variable that plays a key role in the prediction of BR P : Table S2 indicates that both L P and A P have contributions on creating 16 rules and leading to 16 linear equations, while P A with the splitting value of 1.9 bar only produced a linear equation. It is inferred from Table S2 that P A has the lowest effects on the prediction of BR P compared to other input variables. (2)

Implementation of Multivariate Adaptive Regression Spline
MARS, as a sort of nonparametric regression model, was proposed by Friedman (1991). This flexible technique can be applied to estimate continuous numeric variables. The MARS model is capable of producing flexible, accurate, and convenient regression models to predict output variables that are continuous and binary. As major merit, MARS is capable of explaining the sophisticated and non-linear relationships governed among input-output vectors of a complicated system. The MARS fundamentally splits the datasets into different regions so that a regression technique to each sub-region is fitted (Yilmaz et al. 2018). The values of breaks among sub-regions are known as "knots", while the term "Basis Function" (BF) is applied to indicate each distinct interval of the independent variables. The basic form of a BF is expressed as, max (0,x-k) or max(0,kx) in which x and k are the independent variables and a threshold value. The overall formulation of MARS, composed of a combination of BFs with linear trend, where θ is the expected output vector, β 0 is the constant coefficient, β i is the coefficient corresponding to the BFs, m is the number of BFs. The development of the MARS technique includes forward and backward phases. During the forward stepwise approach, all possible BFs are acquired; thereafter, in the backward stepwise sense, BFs which lead to overfitting are removed to improve the accuracy level of the Eq. (5).
Equations given by the MARS model were presented in Table S3. As seen in Table S3, it can be said that L P and A P are the first influential parameters that were used to split the search space. After that, the pressure in operation is the effective variable that divides the search space into smaller spaces. In this study, Eq. (4) is developed using the forward step and 30 BFs are acquired. Thereafter, four BFs were removed through the backward step due to less probability of overfitting. The overfitting process is automatically performed in the MATLAB programming codes by "Generalized Cross-validation (GCV)" formulation along with K-fold value = 10. In this study, the occurrence of the overfitting process was successfully checked out by using GCV in the MARS development. Furthermore, it was assumed that the second-order polynomial used in the MARS model can provide more accurate predictions of pipe break rate than cubic polynomial (3rd order polynomial). This assumption has been widely used in modeling various processes for different engineering problems (Yilmaz et al. 2018).

Implementation of Gene-Expression Programming
GEP, as a major development of GP, has been proposed on the basis of the population evolutionary theorem (Ferreira 2006). GEP applies both intrinsic properties of GA and GP in a way that the simple linear chromosomes with fixed length (genome) and parse trees (GP). In GEP, the genetic parameters that are needed to be defined are comparatively the same as ones in the GP. These parameters are in a variety range: terminal set (i.e., times, minus, plus, square, tanh, exp), function set (known as mathematical functions), type of fitness function, control parameters (e.g., mutation rate, number of population, and generation), and termination conditions for GEP development. GEP includes a character string of a specific length. Individuals which are encoded as linear strings of fixed size are ultimately expressed as nonlinear entities of various sizes and configurations introduced as Expression Trees (ETs) (Saridemir 2010). In the process of GEP development, chromosomes are created randomly with a certain length for entire individuals; then, the ETs are extracted from the chromosomes and the value of fitness function is assessed for each individual. The best value of fitness function associated with individuals is assigned for carrying out the reproduction stage. The GEP development will continue with new individuals for a certain number of generations until the best solution is acquired. In the GEP model, one of the most efficient methods to prevent overfitting is the holdout data in which the split ratio of pipe failure records is 80% for the training stage and 20% for the testing stage. GEP model has been trained until it performs well for both training and testing phases. This demonstrated good generalization capability since the testing stage represents unseen data that were not applied for the training stage. Hence, a sufficiently large pipe break records to train on even after splitting was used for checking the overfitting of the GEP model. Inner functions used in the GEP-based formulation can follow exponential and tangential trends for the prediction of pipe break rates. As a rigorous assumption, this selection had already been proved by previous investigations (Xu et al. 2011a, b;Xu et al. 2013).
The GEP technique acquired the most accurate formulation for the prediction of pipe break rate. To reach this goal, genetic operations, as presented in Table S4, were acquired through the optimal evolution strategy. GeneXproTools5 was used to develop the GEP technique. GEP-based-relationship, developed by 150 generations and three genes, was expressed as,

Evaluation of Accuracy Levels for Data-Driven Approaches
The performance of the AI models is statistically assessed in this part. In this way, the in which NO is the number of observations. The R, as a standard benchmark for the evaluation of error values obtained by the predictive techniques, varies between -1 and 1. Values of ± 1 indicate the most precise prediction while existing direct and invert relationships between observations and predicted values by the predictive models, respectively. Zero value indicates no correlation between observations and output of predictive techniques. Additionally, the highest rank of accuracy value of DDR is zero. When DDR is greater than zero, the AI models illustrate overestimation. Similarly, the predictive technique provides underestimation as DDR is lower than 1. The statistical results of the AI models associated with training and testing stages were provided.
in Table 1 Table 1 indicates satisfying results for both GEP (R = 0.971 and RMSE = 0.544) and MARS (R = 0.981 and RMSE = 0.544) models whereas the M5 has no sufficient capability for the estimation of pipe break rate with R = 0.888 and RMSE = 1.0960. Additionally, DDR values obtained by the GEP (DDR = -0.951) technique prove high capability in the prediction of BR P values in comparison with MARS (DRR = -1.022) and M5 (DDR = -1.005). Due to the negative values of DDR, all the AI models demonstrated comparatively under prediction of pipes break rates. The performance of the AI models for the prediction of pipes break rate in the training and testing stages in Figs. 1-3. As seen in Fig. 1a, the remarkable quality of MARS performance in the training stage is inferred whereas, from Fig. 1b, the testing performance indicates overestimation for BR P = 7.67. Figure 2a depicts highly permissible performance for the GEP model in the training stage whereas, in Fig. 2b, insignificant under predictions for approximately BR P = 3.5-7.0 is seen in the testing performance. Figure 3a demonstrates the quality of M5 performance in the training stage. From Fig. 3a, it is seen that M5 has satisfying performance for BR P ≤ 10. M5 illustrates over prediction for BR P = 20-30 whereas Fig. 3a depicts significant underestimation only for BR P = 20. What is more, Fig. 3b indicates relatively significant under prediction for BR P = 22.3 in the testing stage. Overall, in the terms of a sound comparison, M5 could not provide satisfying performance for a large amount of BR P (77.7) in the training stage whereas MARS and GEP models have significant efficiency in this case. In this way, this factor leads to the decreasing accuracy level of M5 performance in the training phase. Additionally, variations of DDR values versus data samples for the AI models were drawn in Figs. S4-S6. As shown in Fig. S4, almost DDR values acquired by the MARS technique (Eq. (5)) were concentrated between 0 and 0.5. This means that the MARS model provided relative over predictions for the pipes break rate. As depicted in Fig. S5, almost DDR values obtained by the GEP model (Eq. (6)) are around zero levels for both training and testing phases. Ultimately, Fig. S6 illustrated that almost DDR values computed by M5 varied between -0.5 and 0. This is inferred that M5 generally under-predicted the pipe failure rate. Generally, explicit equations given by MARS and GEP techniques were statistically more reliable than M5 Eq. (3). Compared to multilinear regression equations by M5, Eqs. (5 and 6) included a more complicated mathematical structure which was obtained by using the physical and operational factors affecting the pipe failure processes. Equations (5 and 6) with a high degree of non-linearity could better understand the complexity of pipes failure rate processes rather than M5 with linear mathematical structure. On the other hand, mathematical structures of Eqs. (5 and 6) were found to be consistent when compared with available AI models-based-equations from literature. In this research, an explicit equation (Eq. (5)) obtained by the MARS model included a set of second-order polynomial functions for the pipe break predictions. The non-linear mathematical structure of Eq. (5) was consistent with those existing equations obtained by Xu et al. (2011a). An explicit relationship Eq. (6) given by the GEP model included exponential inner functions as applied in the explicit equations by Berardi et al. (2008), Xu et al. (2011b), and Savic et al. (2006), and Wang et al. (2009).

Comparison of the Present Study with Previous Investigations
In this section, comparisons of the current study with previous literature were made based on three strategies: (i) selection of typical effective variables (pipe diameter, pipe age, pipe length, for instance), (ii) complexity of empirical equations, and (iii) merits of AI models. In this section, the results of AI models were compared with those presented in the literature. The first comparison is related to the Berardi et al. (2008) investigations in which an explicit equation developed by the EPR model obtained R = 0.926. Although they did not consider the depth of pipe installation (h P ), they used both physical and operational factors. In this way, it can be said that Eqs. (5 and 6) were relatively more accurate than that reported by Berardi et al. (2008). Moreover, Wang et al. (2009) used logarithmic expressions to predict break rate pipe with different materials: Gray Cast Iron (R = 0.83), Ductile Iron without lining (R = 0.806), Ductile Iron with lining (R = 0.845), Polyvinyl Chloride (R = 0.888), and Hyperscon (R = 0.901). They did not apply operational factors (e.g., pressure fluctuation, number of the failure, depth of pipe installation) to estimate pipe break rate. By the way, equations by Wang et al. (2009) could not address well the complexity of pipe failure processes. Notably, explicit equations by AI models, in this research, are used for all typical materials of water pipes in Yazd's WDN. AI models-based-relationships (Eqs. (5 and 6)) was not limited to the typical material rather than equations reported in the  Wang et al. (2009) research. Additionally, the results of the current study were comparable in terms of accuracy level with those obtained by Xu et al. (2011b). Regardless of operational and environmental factors, they presented explicit formulations by GP (R = 0.85) and EPR (R = 0.84) models to predict the pipe failure rate. The non-linearity degree of Eqs. (5 and 6) is higher than the relationships provided by Xu et al. (2011b) and additionally, this issue causes the increase of Eqs. (5 and 6) accuracy level in comparison with those obtained by Xu et al. (2011b). In a comprehensive comparison, the Automatic Learning Bayesian Network (ALBN), designed by Tang et al. (2019), obtained a 0.82 accuracy level while considering physical, operational, and environmental factors. In the present investigation, the results of AI models' performance were statistically comparable in terms of complexity of influential factor selection with those reported in Tang  Ultimately, the results of AI models-based equations are compared with soft computing models applied in the Robles-Velasco et al.'s (2020) investigation. Accuracy levels associated with LR and SVC models were 0.794 and 0.796, standing relatively lower precision level than Eqs. (5 and 6). It seems that the complexity of SVC and LR was not sufficient to provide more accurate predictions of pipe break rate because Robles-Velasco et al.'s (2020) did not consider important operational factors such as the number of pipe failures and depth of pipe installation. On the contrary, the present study included operational factors and provided more accurate predictions rather than those obtained by Robles-Velasco et al.'s (2020).

Generalization of AI Models Results
In this section, the application of the AI models Eqs. (5 and 6) were evaluated to predict the pipes break rates in the eighth zone of Yazd's WDN. Equations (5 and 6) were fed by the database associated with the Jomhouri zone (Z8). Figure  generally over-predicts (DDR = 0.0613) the pipes break rate whereas almost the BR P values predicted by the GEP model had overestimation with DDR of 1.237. According to the statistical measures, Eq. (5), obtained by the MARS technique, has better efficiency (R = 0.878 and RMSE = 0.416) in the prediction of pipe break rate in the Z8 of WDN than those yielded by the GEP model (R = 0.63 and RMSE = 2.436). The comparatively permissible capability of the MARS model encourages practicability of MARS in acquiring preliminary pipes failure rate values in field investigations while generalizing the AI models by unseen (unreported) pipe failures datasets.

Probabilistic Analysis of Break Rate Model
The probabilistic technique is based on the acquiring negative value (known as violent) for the Limit State Function (LSF) in any system. For a given problem, LSF is the difference between available capacity for output of system and values obtained by mathematical expressions govern on system variables. Therefore, the probability of the limit state of violation is described as (Zampieri et al. 2016;Homaei and Najafzadeh 2020), where G(X) is the Limit State Function (LSF). This function assumes a negative or zero value on failure and a positive value for the safety of the system. Reliability of the system is expressed as, First, LSF needs to be defined for reliability analysis. In the current study, the LSF is formulated as, in which BR C P and BR M P denote BR P values related to the capacity of the WDN case study and BR P given by a mathematical expression. Additionally, SF is the safety level of BR P values which can be assigned greater than 1.
In this study, the results of the AI models are used to express LSF. More specifically, the most accurate equation given by AI models was applied. Accordingly, statistical results of Table 1 indicated that MARS-derived-equation had the high potential of expressing BR M P . To compute BR C P , maximum values of BR P observed in the seven zones of the WDN case study are defined. In reliability analysis, there is a need to determine random variables such as length of pipe (L P ), depth of pipe installation (h P ), allowable pressure in the operation status (P A ), and age of pipe (A P ). Moreover, typical statistical distributions govern by random variables were determined using Kolmogorov-Smirnov (K-S) test. Results of the K-S test have been given in Table S5.
According to Eq. (5), integration of Joint Probability Density Function (JPDF), introduced as LSF, is considered as the probability of failure. Values of P f are computed by using two approximation functions: First-Order Reliability Method (FORM) and Second-Order Reliability Method (SORM). The FORM and SORM employ approximations of the first and second orders Taylor expansion. In this study, FORM is applied to perform the reliability analysis. The results of reliability analysis have been expressed in terms of (10) reliability index (λ) and probability of failure in Table 2. It is inferred from

Statistical Analysis of AI Models
In this section, sustainability related to the performance of AI techniques was studied by using four probabilistic terms: uncertainty, reliability, and resilience, and F-test. These statistical indices are useful to conceptually detect variations of break rate of pipes over time.
The first criterion, uncertainty analysis with a 95% confidence level (U95) is capable of validating the results of the AI techniques for both training and testing stages. U95 restrains the uncertainty of BR P values at a 95% confidence level. Thus, the lower the values of U95, the more precise the BR P amount. Another robustness index, introduced as the reliability criterion, takes into account the overall consistency of the AI techniques. The reliability criterion is computed based on the random error values from the measurement process. The higher is the number of occurrences for which the error amount is lower than a specific value of the threshold, the higher reliable the overall consistency of the AI models would be. The final criterion is the resilience analysis considering the potential of BR P to resists stressors, adapt, and quickly recover from disruptions. Notably, a higher value of resilience  Table 3 shows that the MARS technique represented the lowest value for U95 (0.201) in the training stage in comparison with other AI models. While GEP (U95 = 0.346) and M5 (U95 = 0.346) were at the same level of uncertainty. Additionally, estimations of BR P given by the MARS model were more reliable (0.70) and resilience (0.685), compared to estimations yielded by GEP (0.649 and 0.650) and M5 (0.741 and 0.722). Table 3 indicated that, on the other hand, the M5 model was more reliable and resilient in the training stage than the performance of the GEP model. In the testing stage, the values of uncertainty with a 95% confidence level yielded by the GEP model made comparatively more resilience predictions of BR P (0.733) in comparison to the values provided by the MARS technique (0.667) and M5 (0.591). Furthermore, this trend was seen for values of reliability although break rate values by M5 were lower uncertain (U95 = 0.399) in the testing stage when compared with GEP (U95 = 0.412) and MARS (U95 = 0.465).
Additionally, F-test was applied to evaluate performance both of AI models. In the F-test, it is claimed that null hypothesis is accepted if F 0 > F α,k',n-p where α is the significant level, k' is the number of independent variables, p is the k' + 1, and n is the size of the data sample. For all the proposed AI techniques, values of α,k', and n-p are equal to 0.05, 4, and 127, respectively. From the F-test, F 0.05, 4, 127 value is 2.37 with respect to F distribution table. For the BR P prediction, with respect to F 0.05, 4, 127 , F 0 values obtained by MARS (F 0 = 1.311) and GEP (F 0 = 2.013) models have accepted the hypothesis and consequently indicated accurate estimation of BR P compared with MT approach (F 0 = 8.892). MT is indicative of being a lower capability in the BR P estimation.

Conclusion
In this research, a field investigation has been carried out to report the pipes break rate in the eight zones of Yazd's WDN. First, the key causes of pipes failure were understood, and then datasets were created to formulate the pipe failure rates by MARS, GEP, and M5 MT techniques. Results of AI models development indicated that MARS and GEP techniques stood at the same level of precision for the WDN calibration (or training stage). Equations driven by GEP and MARS techniques were used to model the pipes break rate for the WDN zones whose information had no contribution to training and testing stages. From comparison the present results with the literature, it was found that two chief factors including the complexity of AI models-based-formulations and availability of influential parameters play a key role in the improving accuracy level of AI models. Performance of the reliability analysis indicated that LSF, composed of safety level and Eq. (5), could provide robust reliability indices for various SF levels. The GEP model provided relatively more resilient and reliable estimations of pipe break rates compared to the MARS and M5 techniques.
As inferred from the performance of MT's results, the application of MT-based linear formulation is not the high capability to predict the pipe break rate. To improve the accuracy level of linear models, it is suggested that the weighting coefficients of linear regression models can be adjusted by evolutionary algorithms such as genetic algorithm (GA), particle swarm optimization (PSO), gravitational search algorithm (GSA), and Ant Colony (ACO). Furthermore, results of improved linear regression models by MT can be evaluated by various statistical tests and reliability analysis.
Ultimately, this research presented structural failure techniques for an individual pipe and the mathematical models on how is practically used in decision-making and future planning.
Authors Contributions Yaser Amiri-Ardakani; Performing the field study and collecting the data; Mohammad. Najafzadeh; Formal analysis and investigation, Writing-original draft preparation, Writing -review and editing, Resources, Supervision.
Funding No funds, grants, or other supports were received.

Data Availability
The data are not publicly available due to restrictions such their containing information that.

Declarations
Ethical Approval All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964. Helsinki declaration and its later amendments or comparable ethical standards.

Consent to Participate
Informed consent was obtained from all individual participants included in the study.

Consent to Publish
All the authors give the Publisher the permission of the authors to publish the research work.
Competing Interests There is no conflict of interest.