Conceptual Water Main Failure Risk: Self-Excitation, Pipe Age, and Statistical Modeling Performance

SPPs can be characterized by their conditional intensity function, which estimates the event occurrence rate at any point in space and time. SPPs typically map to two-dimensional Euclidean space, although significant progress has been made recently in their application to linear networks (Baddeley et al. 2021). Because this model applies to a network of pipes whose composition is variable over time (pipes can be installed or abandoned/removed and occupy the same physical location at the same time), the spatial coordinates refer to each unique pipe ($i$), with spatial information associated with each pipe included as covariates. Each unique pipe ($i$) consists of a continuous segment of pipe as found in the geospatial pipe data set provided by a water utility. The model presented here is a SPP with pipes $i\in P=\mathbb{Z}[0,n]$, events on pipes $j\subset P$ with event times $t_{j,k}\in \mathbb{R}[0,T+1)$ (where $T$ is the last time step in the study and 0 is the beginning of the study observation window, not the age of the pipe) with event number $k\in \mathbb{Z}[0,\infty ]$, time steps $\tau \in \mathbb{Z}[0,T]$, and pipe history ${\mathcal{H}}_i$. The event times ($t_{j,k}$) consist of real numbers and are thus not restricted to the integer domain of the time steps for the time-varying covariates. For example, given a group of five pipes ($P=\mathbb{Z}[0,4]$), one might observe the sequence of breaks shown in Table 1.

The conditional intensity for pipe $i$ at time $t$ given its history ${\mathcal{H}}_i$ takes the formwith ${\mathcal{H}}_i$ implied hereafter. The background intensity function $\mu (i,t)$ is shown in Eq. (2), where ${\beta }^{\tau }$ and $\beta$ are vectors of $n_1$ and $n_2$ parameters, respectively, and $X^{\tau }$ and $X$ are matrices of time-variant and time-invariant covariates, respectively.

The triggering function $g(·)$ accounts for the self-excitation of breaks on the same pipe and is commonly specified in the exponential form, i.e., where $\alpha$ controls the increase in intensity due to an event and $\gamma$ controls the decay (Laub et al. 2015).

This model assumes a background intensity based on spatiotemporal covariates and a self-exciting triggering intensity that decays exponentially in time, a specification also assumed by Reinhart and Greenhouse (2018), although with the addition of a Gaussian spatial triggering term not considered here.

The parameters of SPP models are estimated by maximizing the log-likelihood, which is adapted from Reinhart and Greenhouse (2018) aswhere $\widehat{\mathrm{\Theta }}=⟨\alpha ,\gamma ,{\beta }_k^{\tau },\dots ,{\beta }_{n_1}^{\tau },{\beta }_k,\dots ,{\beta }_{n_2}⟩$; ${\tau }_{i,o}$ is the first time step of the study in which pipe $i$ is active; and ${\tau }_{i,f}$ is the last time step of the study in which pipe $i$ is active. Maximization is performed using either sequential least-squares programming (SLSQP) or the modified Powell algorithm as implemented in the Python package SciPy (Virtanen et al. 2020). As noted earlier, the expression for the log-likelihood given by Eq. (4) assumes that the complete break history is known (Reinhart 2018). To improve problem solvability and for consistency with the conceptual model shown in Fig. 1, $\alpha$ and $\gamma$ were limited to positive nonzero values.

Confidence intervals for the fitted parameters were produced with the asymptotic covariance matrix estimatorwhere the list of all $j$ is the list of all failure events (with more than one entry per $j$ when a main has broken more than once); the value of $j$ is the pipe identification number; and the matrix-valued function $\mathrm{\Delta }(j,t_j)$ has $m$ rows and $p$ columns (where $m$ and $p$ refer to a given parameter within $\widehat{\mathrm{\Theta }}$) and is given by

This method was suggested by Rathbun (1996) and used by Reinhart and Greenhouse (2018). There are two other options for obtaining confidence intervals: (1) calculating the Hessian of the log-likelihood and (2) a parametric bootstrap method. The bootstrap approach is computationally intractable for large data sets, and the Hessian method was not as robust as using Rathbun’s estimator in Reinhart’s analysis (Reinhart and Greenhouse 2018). The 95% confidence intervals for the parameter estimates are produced by multiplying the estimated standard errors by a critical value of 1.96 (Hazra 2017). We assume that regularity conditions apply: the true parameters must be interior to the parameter space, the log-likelihood function must be thrice differentiable, and the third derivatives must be bounded.

Model selection was performed by fitting the model with all available covariates and subsequently removing those covariates with the most uncertain coefficients based on a Wald test (Agresti 2007). Parameters whose estimates have $p$-values less than 0.1 are included in the model.

Several published models were evaluated and compared to the proposed SPP model using identical estimation and simulation procedures and covariates. The model given by Scheidegger et al. (2013) assumes the time to the first failure is Weibull distributed, while the subsequent failures are assumed to follow a single exponential distribution. The conditional intensity for this model, modified to include covariates following Scheidegger et al. (2015), is given bywhere $n(i,t)$ is the previous number of breaks at time $t$. Scheidegger et al. (2013) modified the maximum likelihood estimation to account for missing pipes; however, the same maximum likelihood estimation approach is used here for all models to ensure the comparisons are valid. While Scheidegger et al. (2013) used Bayesian methods, the data sets used here are sufficiently large to allay concerns over discrepancies between the two approaches. Unfortunately, significant missing failure history means that an assumption is required to implement Scheidegger et al.’s (2013) model: the first observed failure is assumed to be the first failure [i.e., $n(i,t)=0$].

Kleiner and Rajani (2010) implemented a nonhomogeneous Poisson model that is essentially the same as the SPP model but with a constant parameter (${\alpha }_o$) and without the triggering function. Their model is given by

Lin and Yuan (2019) proposed a type of two-time-scale point process model that incorporates covariates, pipe age (included here in $X_{i,k}^{\tau }$), and time since the last break in a conditional intensity function similar in form to

We assume that if there are no previous breaks, then the term accounting for break history has a value of 1.

To evaluate the performance of the models with and without self-excitation, modifications were made as follows. Scheidegger et al.’s (2013) model was compared with and without the triggering function shown in Eq. (3). Kleiner and Rajani’s (2010) model was compared with the SPP model because they are nearly identical aside from the inclusion of self-excitation in the SPP model. Lin and Yuan’s (2019) model was compared with and without its original self-excitation function as shown in the right half of Eq. (9).

Simulations of failure events were obtained by implementing the method outlined in Section 5 in Zhuang et al. (2004), which involves first generating a background catalog of events based on the specified background intensity ($\mu$) and subsequently generating the so-called children for each generation of events based on the specified triggering function’s ($g$) normalized probability density function. The background nonhomogeneous stationary Poisson process was generated using Lewis’s thinning algorithm (Lewis and Shedler 1979; Zhuang and Touati 2019). Rather than simulating the break behavior for all pipes over the entire study period in a single simulation, which could lead to error propagation, the method used here draws upon the approach given by Algorithm B in Section 3.3 of Zhuang (2011) and simulates breaks in a stepwise fashion, with the true break history used to simulate future breaks. For example, given the true break history up to, but not including, time $t$, the break behavior for $[t,t+\mathrm{\Delta }t]$ was simulated using the true break history and the method outlined earlier. The training period consisted of monthly stepwise simulations performed until January 2018. The rest of the study period (from January 2018 through December 2020) was used as a single forecasting period (i.e., the testing period). The same simulation procedure was used for all models. However, the algorithm used to generate the children was only implemented where the triggering function is additive [e.g., see Eq. (1)]. Further details regarding the simulation algorithm may be found in the Supplementary Information.

In addition to providing information on factors related to pipe breaks, the models presented here provide additional information in the form of pipe-level break probabilities. Each simulation based on a fitted model is a realization of the modeled stochastic process for each pipe segment. Whether or not a given pipe breaks in each simulation depends on its conditional intensity over the forecasting period (break probability) and the randomness of the simulation. Thus, when many simulations are run over a certain time period (i.e., a Monte Carlo simulation), the pipes that are simulated as breaking are likely to break only in some fraction of those simulations. If a pipe breaks in 60% of the simulations, it is assigned a break probability of 0.6. This value may be used to classify the pipe segment as either breaking or not breaking during the forecast interval, depending on the classification threshold. This information may be useful to water utilities seeking accurate, pipe-level predictions for pipeline replacement programs. The pipe break predictions from the training period were not used for this aspect of the analysis.

Classifiers are often evaluated using receiver operating characteristic (ROC) curves, which plot the true positive rate versus the false positive rate. However, ROC curves are flawed measures of classification performance when the data are highly imbalanced (i.e., a small number of positives versus negatives) (Lever et al. 2016). Because main breaks are rare when considering the total number of pipes that are active within a distribution system, precision-recall (PR) curves are more appropriate for measuring the performance of a main break classifier. In plots of precision versus recall, superior classifiers produce curves closer to the top right. The area under the PR curve (AUC-PR) is a way to quantify the relative performance of each classifier, with greater AUC-PR indicating better performance. A random classifier will result in an AUC-PR approximately equal to the positive rate in the sample (Sofaer et al. 2018). AUC-PR values were calculated using the trapezoidal rule as implemented in the Python package NumPy (Harris et al. 2020).

The distribution network and main break data used for this analysis were provided by East Bay Municipal Utility District (EBMUD), a water agency serving 1.4 million customers in California, USA. Main break data were provided in the form of a geospatial data set containing the location and date reported for each break in addition to the attributes of the pipe on which the break occurred, the failure mode, comments by the maintenance crew, and other identifying information. Additionally, a geospatial pipeline database was provided that gives the location and attributes of each pipe, here defined as the unit of pipe assigned a unique identification number by EBMUD in the provided data sets, including material, diameter, length, installation year, life cycle status, life cycle status change date, pressure zone, and others. However, there is no direct link between the main break data set and the provided pipeline data set, meaning we cannot assign breaks to pipes with full certainty. This means that the two data sets, one for breaks and one for mains, must be combined to obtain the best available estimate of the break history for each pipe segment and to infer missing pipes.

To obtain the best available estimate of each pipe’s break history, an automated algorithm was used to match breaks to pipes and if necessary infer pipes based on the available evidence. EBMUD has confidence in the geospatial fidelity of the mains to within approximately $\pm 1\mathrm{ft}$, and the main breaks are associated with the nearest address. The laterals, which connect mains and service connections (with addresses), were used to adjust the locations of the main breaks to overlap the mains. The main breaks were then assigned to the most likely nearest pipe based on matching pipe characteristics, for example, the estimated installation year and pipe material. If no suitable pipe was found for a given main break, a main was inferred based on the estimated pipe characteristics found in the main break records.

Installation practices vary based on the type of cover (e.g., asphalt versus soil only) and the type of road (e.g., highway versus local street) under which the pipe is located. To reduce variability among the pipes considered for modeling, only pipes adjacent to or buried under local streets were used. Local streets were defined as roads bearing a Caltrans Functional Classification System value of 7. This group of pipes is the largest in EBMUD’s distribution network (76% by length). Additionally, only main breaks assigned the highest priority levels of 4 and 5 (91% of all main breaks) were modeled. EBMUD assigns priority from 1 to 5 to each break based on severity (with priority 5 breaks most severe). This is intended to mitigate the error associated with breaks that are reported long after they actually occur, which may be more likely for breaks of lower severity. For example, a catastrophic main break (likely priority five) is more likely to be immediately noticed and reported by the public than a very small leak (likely a much lower priority). The lag time between the date a break is reported and the date the break actually occurred is of great interest to utilities and could be included in the model by adjusting break times as necessary. However, this is a difficult problem to solve, and since no solution has been found, it must be accepted as contributing to the error of any pipe break model.

Several common deficiencies of distribution network and break history data must be resolved for modeling to proceed. Two widespread issues are left truncation and survival bias (Scheidegger et al. 2015). The first consists of missing break history, which is an inescapable flaw of main break data collected by most water utilities that began operating long before the digital age. This is especially relevant for EBMUD given that pipes were installed in its service territory as early as 1866, while digital break records only go back to 1997. This means that of the approximately 155 years over which main breaks have occurred, records only exist for approximately 23 years, or 15%.

While Lin and Yuan (2019) simulated the missing history, pipe break modelers often choose to ignore this problem by assuming no dependence on past breaks or by including a variable accounting for the number of known breaks. The relevant question is whether or not the missing break history is important and, if so, to what extent. For the SPP model, instead of including the total number of breaks experienced by a pipe at any given time, the effect of past break history is described by the triggering function [Eq. (3)]. Problems arise when events are triggered by an event outside of the temporal boundary, but the magnitude of this problem can be assessed though the triggering parameters and estimating the duration that the temporary increase in risk due to a break is significant. The data shown in Fig. 2 suggest that the self-excitation effect is relatively short, on the order of approximately 5 to 15 years depending on the material. This suggests that the missing break history beyond a few years near the boundary will have minimal impact on the model and can be safely neglected. The error due to the missing break history near the temporal boundary is accepted as a minor source of error. This is a common approach in earthquake and crime modeling, where an additional assumption must also be made about events being triggered from outside the spatial boundaries.

The second issue (survival bias) is the absence of information about abandoned or replaced pipes in the database of pipes, which means that only pipes that have survived are observed. Because the surviving pipes may not be representative of the entire cohort, the results of analyses using only surviving pipes may not be generalizable. Fortunately, EBMUD’s records contain information on many abandoned or removed pipes. The remaining error due to this flaw in the data is partially mitigated by inferring the presence of missing pipes using the break data. The algorithm that assigns each break to a single pipe segment also flags breaks with no satisfactory pipe assignment. Using the pipe characteristics associated with the flagged break, the missing pipe is generated and is assumed to be collocated with the nearest pipe segment. However, this procedure does not mitigate the error associated with missing data on pipes that were abandoned before digital records began or pipes that were abandoned after digital records began but did not fail during that time period and were not kept in the pipe database. For more information on the break-pipe assignment algorithm, see the Supplementary Information.

Air temperature and precipitation data were obtained from the National Oceanic and Atmospheric Administration Global Historical Climatology Network (Menne et al. 2018) and the stations used for this analysis are all those within 16.1 km (10 mi) of the pressure zones operated by EBMUD. There are 15 temperature stations within 16.1 km (10 mi) of the service territory. Only nine of those have full coverage of the study period, and the others are fairly sparse. The missing data were filled with that station’s nearest neighbor data. The same approach was used for the 88 stations with precipitation data within 16.1 km (10 mi) of EBMUD’s service territory, many of which do not have complete records. Each pipe was assigned the weather data of the nearest station.

GIS data for California’s roads were obtained through the freely available Functional Classification data set published by Caltrans (2021). Each pipe was assigned the type of road to which it is nearest, which consists of a value of 1 (interstate), 2 (other freeway or expressway), 3 (principal arterial), 4 (minor arterial), 5 (major collector), 6 (minor collector), or 7 (local). Additionally, monthly geospatial traffic volume (estimated trip counts) data were obtained from StreetLight Data (San Francisco, California, USA) spanning from 2018-07-01 through 2020-01-31 (the date format used in this work is YYYY-MM-DD in accordance with ISO 8601). Owing to the number of pipe segments studied, obtaining unique traffic data for each pipe is not feasible. Thus, the traffic data were aggregated to spatiotemporal grid cells and normalized by the area of the cell. The cells were designed to minimize misrepresentation of the traffic experienced by each pipe by isolating high-traffic roads within their own polygons. The polygons were generated using Voronoi tessellation and manual modification as necessary to include all pipes and to meet the maximum polygon limit imposed by Streetlight. The traffic data contain distinctions for pass-through, end-in, and start-in traffic. The pass-through traffic is distinct from the other types and may be attributed to arterials and freeways rather than local roads. Accordingly, pipes under local roads were assigned the traffic volume values from the start-in and end-in categories. In addition to the time-varying traffic covariate, a per-pipe aggregated value was also included under the assumption that the value is valid beyond the dates for which traffic data are available. All traffic covariates are normalized by the area of the grid cell to which they belong.

Soil data were obtained from the USGS gNATSGO (Soil Survey Staff 2021) data set and assigned to each pipe using ArcMap (Esri, Redlands, California, USA). Unfortunately, soil properties cannot be assigned to all pipes due to the incomplete coverage of EBMUD’s service territory. Soil clay content, electrical conductivity, and pH are available for 95%, 71%, and 76% of the mains, respectively. When these variables were included, the pipes without data were ignored.

Monthly water demand data were provided by EBMUD for all of its customers from 2005-01-01 through 2021-01-31. The addresses were geolocated and used to assign a local demand surrogate for each pipe. To capture all pipes using this method, the 50 nearest customer taps within a 16.1 km (10 mi) radius of each pipe were used to calculate the time-varying demand surrogate. Additionally, a single, all-time average per-pipe surrogate demand value was included as a time-invariant covariate to characterize the general demand trends for each pipe and to isolate the fluctuations of demand over time. A very small group of pipes (0.34% of all pipes) has much greater surrogate water demand than most other pipes. These pipes are considered outliers or not representative and were excluded from modeling.

The estimated model parameters are shown in Table S7. The coefficients may be interpreted relative to each other within each material type because the covariates were all scaled to the same range (between 1 and 2). The self-excitation parameters were statistically significant (95% confidence level) in all cases, suggesting that the observed temporal clustering could be described by self-excitation. Traffic volume was evaluated as a time-varying covariate to the extent permitted by data availability, which was limited to 2018-07-01 to 2020-01-31. Traffic volume failed to produce a significant parameter estimate at the 95% confidence level, with the exception of traffic for CI main breaks. Traffic appears to be significant for CI main breaks at the 95% confidence level, but the resulting coefficient is negative. This implies that increased traffic is associated with reduced break risk. Although unintuitive, this result is consistent with findings by Eisenbeis (1994). To facilitate comparison between material types and to take advantage of the longer time span available for other covariates, traffic was removed from consideration, and the study proceeded with the full set of covariates available from 2013-01-01 through 2020-12-31, as shown in Table S7. Some results are intuitive, including the negative coefficients for diameter (larger-diameter pipes are stronger due to their larger area moment of inertia), positive coefficients for length (longer pipes present more opportunities to break), the positive coefficients for soil clay content (higher clay contents result in more extreme soil expansion and contraction), and the positive coefficient for air temperature for AC mains (in the study area, higher temperatures coincide with arid conditions in which differential soil settlement occurs due to variations in moisture content). Contrastingly, the negative coefficient for the age of the pipe for CI and AC mains is unintuitive when considering only the conceptual failure model described earlier. However, the empirical failure risk shown in Fig. 5 suggests that any single coefficient for age may be inadequate, as failure risk does not necessarily increase monotonically with age. For an extended discussion of the estimation and inference results, see the Supplemental Materials.

The network-wide predictions resulting from the SPP model fitted to CI, AC, and ST main breaks are compared to the observed data in Fig. 6. The models for CI and ST capture the average break behavior, while the model for AC better describes the seasonal variation in breaks. As suggested by the different seasonal trends shown in Fig. S5, CI breaks are likely driven by many different factors. The models presented here do not distinguish between the failure modes, and thus the relationships between covariates and main breaks may be different for different types of breaks. This may explain why the simulated seasonal behavior is damped for the CI model. Conversely, the majority of AC main breaks appear to be correlated with a single seasonal variable, as shown in Fig. S6. Accordingly, the simulations show good agreement with the observed seasonal patterns of AC main breaks. ST main breaks do not display obvious seasonal patterns and appear to be much more random than their CI or AC counterparts, and the fitted model for ST breaks reflects this lack of seasonality (Fig. S7). The network-wide main break predictions resulting from Scheidegger et al.’s (2013), Lin and Yuan’s (2019), and Kleiner and Rajani’s (2010) models are shown in Figs. S8–S10. Kleiner and Rajani’s (2010) model resulted in the lowest mean square error (MSE) for CI, followed closely by the SPP model and Lin and Yuan’s (2010) model. For AC, the SPP model provided the lowest MSE, while Lin and Yuan’s (2019) model provided the lowest MSE for ST. It appears that including self-excitation resulted in higher MSEs in some cases while reducing MSEs in others (Fig. 7).

The PR curves for all six models (three with and three without triggering functions to model self-excitation) are shown in Fig. S11. The AUC-PR for all models are summarized in the top row of Fig. 7. A larger AUC-PR indicates better classification performance. The data suggest that only Scheidegger et al.’s (2013) model is superior without triggering and that predictions of CI main breaks are superior to predictions of AC and ST main breaks. It appears that although self-excitation is observed in the empirical data, including it in the model definition does not always improve failure predictions. Moreover, the second-highest performing model across all materials (based on Fig. 7), Scheidegger et al.’s (2013) original model (without self-excitation), performs nearly as well as Lin and Yuan’s (2019) model with self-excitation. The significant additional computational effort to include triggering may not be worth the small difference in classification performance between the two models. Interestingly, the models that proved superior at classifying individual failures were not necessarily better at predicting network-wide behavior. For example, as shown in Fig. 7, Scheidegger et al.’s (2013) model applied to CI, AC, and ST failures was nearly the best at classifying individual breaks (according to the AUC-PR) but in all cases produced larger MSEs for network-wide predictions compared to the other models.