Calibration of intermittent water supply systems hydraulic models under data scarcity

About 1.3 billion people mainly in South Asia, Latin America, and Africa receive water through intermittent water supply systems (IWSSs) (Farmani et al. 2021). IWSSs are unable to meet customer demands due to water scarcity from the sources or due to economic or/and technical scarcity. Due to water scarcity, these piped water supply systems (WSSs) deliver water to users for fewer than 24 h a day or less than 7 days a week (Kumpel & Nelson 2016). IWSSs have traditionally been designed as continuous water supply systems (CWSSs). Consequently, their operation in a way other than their design concept presents several technical challenges (Klingel & Nestmann 2014). These include increased pipe bursts and background leaks due to pressure surges and fluctuations (Zyoud 2022). The high leakage levels contribute to the high non-revenue water (NRW) prevalent in IWSSs. Moreover, there are high flow rates caused by high load factors during short supply periods (Zyoud 2022). The high flow rates induce excessive friction head losses leading to pressure deficiencies that are the major causes of water quality and inequitable water distribution concerns in these systems (Vairavamoorthy et al. 2007; Ingeduld et al. 2008). Due to their wide application and complexity, IWSSs are a topical research subject (IWA 2023).

According to Zyoud (2022), IWS research can be divided into three major categories, which are (i) causal factors and transition from IWS to continuous water supply (CWS), (ii) water contamination and health risks associated with IWS, and (iii) design and optimisation techniques applied to these systems. Several studies focus on developing, using hydraulic models, optimal designs and operation strategies for equitable water distribution (Ameyaw et al. 2013; Ilaya-Ayza et al. 2017; Solgi et al. 2020; Ceita et al. 2023). The challenge with modelling IWSSs has been the lack of appropriate modelling techniques. Demand-driven analysis (DDA) models, which are usually used with CWSSs, are not suitable for IWSSs. They lead to excessive simulated flows because they indicate that demands are met even at nodes where pressure is too low to supply water. This inaccurate representation of IWSSs by DDA models creates the need for pressure-driven analysis (PDA) hydraulic solvers. EPANET 2.2 solver has both DDA and PDA capabilities (Rossman et al. 2020).

In the EPANET 2.2 PDA option, the actual node consumption demand outflow delivered at the node is a function of the pressure available and the base demand (Rossman et al. 2020). To realistically simulate the hydraulic behaviour of a WSS, NRW must be incorporated in the model because the total water demand comprises consumption water demand and NRW (Butler 2004). However, to the best knowledge of the authors, there are no mathematical models for modelling NRW in the literature, while there are such models for leakage (Germanopoulos 1985; Gupta et al. 2016; Lambert et al. 2017; Simukonda 2021). Consequently, in this study, leakage is assumed to represent NRW because of the following:

• (a)

  The general definition of NRW is the difference between the amount of water put into the distribution system (system input volume) and the amount of water billed to consumers (Billed Authorised Consumption) (Kingdom et al. 2006; Leakssuite Library Ltd 2020). This definition of NRW, considers the total volume of water lost as NRW without necessarily considering how the water is lost (real (physical) or apparent water losses). Moreover, in both the Standard IWA Water Balance and the IWA Water Balance with Additions (Leakssuite Library Ltd 2020), the only apparent water loss component that does not necessarily involve physical water flows through pipes is customer metering inaccuracies and data handling/billing errors. All the other apparent water loss constituents involve water flows through pipes and affect their hydraulic capacity just like real losses. The definition of NRW and this view of apparent water losses from the pipe hydraulics perspective link NRW and leakage in agreement with Anand (2017) who defines leakages as flows of water in the city pipes that are not fully authorised, controlled, and unknown by the city authorities.

• (b)

  The scarcity of data on IWSSs (Hayuti et al. 2007; Klingel & Nestmann 2014; Sarisen et al. 2022) makes it difficult to find instances where NRW for these systems is differentiated into real and apparent water losses. This is critical especially regarding billing errors. According to Mastaller & Klingel (2018), due to less than 100% metering, the determination of the actual apparent water losses is a challenge because of the possibility that those on fixed charge, but supplied with excess water, consume more water than they are billed and those supplied with less water consume less water than billed.

• (c)

  Galaitsi et al. (2016) indicated that estimations show that half of the water losses are due to leaks with the other half being apparent losses caused mainly by illegal connections. Illegal connections may include connecting households to the pipe network without the knowledge of the utility company or may be part of unauthorised consumption aspects which include consumers who reconnect themselves after being disconnected by the water utility, disconnected consumers resorting to getting water from neighbours who are on fixed charge resulting in one bill for the household with the fixed charge connection but the water consumed could be for two or more households, and consumers with missing data in the billing system, who may be legally connected but because there are no records of their property in the utility billing system, the consumers use water for free (Simukonda 2021).

As highlighted in part (a), all these constituents of apparent water losses involve physical flows of water out of the pipes. Consequently, for hydraulic modelling, which is the focus of the current study, discussions presented in parts (a) to (c) give the basis for assuming the synonymity of NRW and leakage. From the law of conservation of mass, the district metered area (DMA) system input volume is equal to its total demand which constitutes leakage (NRW) and consumption demand. Thus, the incorporation of NRW in IWSSs hydraulic models is important because, for these systems, NRW constitutes a significant part of pipe flows. For instance, in India, NRW varies from 15 to 70% of the system input volume (Vishe 2019). In Zambian utilities, NRW levels range from 33 to 74% (NWASCO 2021). However, modelling leakage, which is representing NRW, is not directly provided for in EPANET. Including leakage in EPANET for PDA simulations is complex because the relationship between leakage and pressure must be defined and a decision must be made whether to place leaks along the pipe or at the nodes (WDSA/CCWI Committee 2022).

There are three commonly used methods for including leakage in WSSs hydraulic models: the additional nodal demand, the artificial reservoir/tank, and the artificial emitter methods (Taylor et al. 2019; Muranho et al. 2014; Cobacho et al. 2015; Sophocleous et al. 2019; Marsili et al. 2023; Mottahedin et al. 2023). Traditionally, modelling leakage involves adding additional demands to nodes. This means each node has a consumption and a leakage demand. The leakage demand is commonly tied to the consumption demand pattern (Cobacho et al. 2015) or kept constant over the simulation duration (Sebbagh et al. 2018). The drawback of this method is that leakage is treated as being independent of pressure ignoring the well-known consumption demand–pressure–leakage relationship. For the artificial reservoir/tank methods, the drawback is that they are uncommon and according to Muranho et al. (2014) reservoir/tank-based methods are computationally more complex than the emitter-based approaches. The artificial emitter-based methods use the fixed and variable area discharge concept (Lambert et al. 2017). This concept uses the emitter equation where the leakage flow rate is a function of pressure at the node, the emitter coefficient, and the emitter exponent (Sebbagh et al. 2018; Wu et al. 2022). The emitter can be added as a node at a selected location along a pipe or it can be added as an artificial element connected to a demand node (WDSA/CCWI Committee 2022).

For hydraulic models to be accurately calibrated, including leakage modelling is a vital step (Berardi & Giustolisi 2021). However, IWSS models that include leakage are uncommon in the literature.

Battermann & Macke (2001) developed a PDA method for IWSS that included leakage modelling. The limitation of this method is that it has numerical convergence problems when applied to large systems and its accuracy is affected by the difficulty in determining the correct sizes of the domestic water consumption tanks and the distribution of the leakage tanks in the hydraulic model (Battermann & Macke 2001). Burrows et al. (2003) presented a method for dynamic leakage representation in network modelling studies using EPANET. However, it is not clear how leakage was modelled using an emitter (PDA) together with consumption demand using DDA at the same node. Ingeduld et al. (2008) developed an IWSS modelling approach that considered both consumption demand and leakage, but limited details about their modelling and calibration approaches were provided. Mohapatra et al. (2014) developed a hydraulic model for assessing both intermittent and continuous WSSs using EPANET considering both consumption demand and leakage. Consumption demand was modelled using artificial reservoirs, but it was not stated how leakage was modelled. Rathi et al. (2020) used an IWSS to demonstrate challenges in the calibration of WSS hydraulic models. However, they first converted the water supply, for the case study zone, from intermittent to CWS and modelled both consumption and leakage demand using the DDA method. It is noteworthy that even for CWSSs, modelling leakage using the conventional DDA methods is unrealistic (Hayuti et al. 2007; Sebbagh et al. 2018). Simukonda (2021) developed an artificial emitter-based method for modelling IWSSs considering the differences between consumption demand and leakage. The limitation of the method is that it was developed before the universal publicity of EPANET 2.2 which takes care of the PDA consumption demand aspects. Moreover, the method used manual calibration which is cumbersome when applied to large systems with several measured data points.

The few studies above and Sarisen et al. (2022) demonstrate that there is limited literature on the comprehensive calibration methodology of IWSSs models in general and those with leakage modelling in particular. This is due to the lack of information for model development and the scarcity of field-measured data for model calibration (Hayuti et al. 2007; Sarisen et al. 2022). However, like any model, calibration of IWSSs models is inevitable because it ensures that simulation results match closely with the measured data (Sophocleous et al. 2017). This study proposes a methodology for calibrating IWSS hydraulic models that include leakage modelling. This is a major contribution to the modelling of IWSSs because in these systems, leakage is normally high and the proposed approach will ensure that hydraulic models represent these systems more realistically (Charalambous & Laspidou 2017). The methodology also includes distinct steps to mitigate the problem of data scarcity. Moreover, the proposed calibration approach selects the emitters' exponent and coefficients automatically through optimisation. This simplifies the task of determining leakage emitter coefficients to the setting of their lower and upper bounds based on the minimum and maximum pressure for each DMA. This is vital because while the bounds for the emitter exponents are well documented in the literature such as Fuchs-Hanusch et al. (2016), those for emitter coefficients are not. Thus, for the emitter coefficients, this study presents an approach on how to estimate their bounds. This is also a significant contribution because these coefficients naturally should be developed iteratively for systems with high leakage levels (Burrows et al. 2003).

This section presents the theoretical background on the PDA modelling of consumption water demand and leakage. It includes a step-by-step IWSS hydraulic model calibration process and the case study to which the methodology is applied.

The emitter (discharge) coefficient, for each leak, increases as the leakage volume from it increases with time (WDSA/CCWI Committee 2022). This makes it difficult to know the range of the emitter coefficient values in advance.

The emitter exponent governs the pressure sensitivity of the leakage flow rate (Zaman et al. 2022). It mainly depends on the predominant pipe material in the water distribution system (WDS) (Cobacho et al. 2015). Its values have been found to be close to 0.5 for round holes in all materials. This makes 0.5 the commonly used value for the leakage emitter exponent (α) and the default value in EPANET (Mahmoud et al. 2017).

The length of the artificial pipe with the CV is short and has a wide diameter to ensure minimal head loss across the artificial elements (Mahmoud et al. 2017) thereby leading to similar pressure values at the demand and emitter nodes. Consequently, the outflow from the emitter depends on the pressure at the demand node. This inclusion of these artificial elements at demand nodes and models leaks as dissipating from the emitters.

The calibration of a WSS hydraulic model can be summarised in seven steps (Ormsbee & Lingireddy 1997): (i) identification of the intended use of the model, (ii) model parameter estimations, (iii) collection of calibration data, (iv) evaluation of the model results, (v) performing macro calibration, (vi) performing sensitivity analysis (SA), and (vii) performing micro calibration.

Regarding the identification of the intended use of the model, the methodology presented (in this study) is for hydraulic calibration of models for long-range planning, designing, and WSSs operation. The estimations of initial model parameters (e.g., pipe roughness and node demand) are done during the initial phase of hydraulic model development. Calibration data can be collected from the utility company records where there are consistent flow and pressure measurements for monitoring the system performance. Data can also be collected from consultant or commissioned research reports. However, for IWS data that can be used for calibration is scarce. This makes it inevitable that, where the data are not available, the modeller should collect the hydraulic model calibration data through field measurements as exemplified in the literature (Battermann & Macke 2001; Ingeduld et al. 2008; O'Donnell 2010; Alves et al. 2014; Klingel & Nestmann 2014).

The first three calibration steps (i, ii, iii) apply to both IWSSs and CWSSs. The difference is in step (iv) (evaluation of the model results). The initial evaluation checks the alignment of DDA and PDA results. In high-pressure areas, the simulation results between DDA and PDA should be identical (Mahmoud et al. 2017). However, if there are significant differences in the results in high-pressure areas, the most likely problem could be the PDA settings for the emitter exponent, minimum pressure, and required pressure. For the emitter exponent and minimum pressure, the default values of 0.5 and 0 are used, respectively. Therefore, the required pressure is the only parameter that needs adjustments to ensure that the DDA and PDA simulation results in high-pressure areas match. The macro calibration, SA, and micro calibration (v, vi, and vii) steps are further explored below.

The premise of macro calibration is the analysis of the initially developed hydraulic model. If the model results greatly differ from the measured data (30% difference or more), there is a need for initial rough tuning of the hydraulic model following the procedure in Ormsbee & Lingireddy (1997). The procedure outlines the possible causes of the large discrepancies and how to identify them. Identifying the main cause of the large errors and reducing the errors is important because it affects the effectiveness of micro calibration techniques (Ormsbee & Lingireddy 1997; Simukonda 2021). This is because through the identification of the causes of large simulation-measured data discrepancies, macro calibration reduces the possibility of compensating errors. These errors arise when values for one parameter are adjusted to force the model simulation results to agree with measured values and yet the cause of the disagreement between simulated and measured values is a different parameter (U. S. Environmental Protection Agency 2005; Walski 2017).

After a reasonable correlation between measured and simulated results is achieved (within an error range of less than 30%) through macro calibration, micro calibration is done (Ormsbee & Lingireddy 1997). After macro calibration, the accuracy of the model is usually not yet good enough, and there are still many parameters (variables) that may cause disagreement between simulated and measured values. To identify the most likely source of error, before performing micro calibration, SA is done to the hydraulic model with accepted settings from macro calibration.

SA is conducted to investigate a system's causality between inputs and outputs. SA benefits include the identification of sensitive modelling parameters that heavily influence the model outputs (Razavi et al. 2021). The identification of these parameters is vital to ensuring the efficient hydraulic model micro calibration (Sophocleous et al. 2017).

SA can be done using local or global methods. Local or one factor at a time (OAT) involves changing values of one independent variable within the range of its variation and evaluating the extent of change in the model output. Global methods involve the changing of values of all independent variables in the model at once (Saltelli & Annoni 2011). The OAT method is the most common in the literature (Saltelli et al. 2004) and it is used in this study.

Pipe roughness and nodal consumption demand are the most uncertain parameters in WSS hydraulic models (Ormsbee & Lingireddy 1997). Demand has a fluctuating behaviour compared to pipe roughness, making real-time calibration an important means of achieving updated models. The combination of a lack of measurement data and high uncertainty poses a challenge to the calibration of WSSs hydraulic models (Sanz & Pérez 2015). SA using the OAT method involves iteratively adjusting the pipe roughness and then demand values or vice versa while evaluating the sensitivity of the model output to parameter variations (U. S. Environmental Protection Agency 2005) using the MAE (Equation (3)).

For models that include leakage, leakage values are also unknown and the addition of uncalibrated leakage parameters through artificial elements influences the model fit significantly; hence these must be calibrated during micro calibration to simulate actual leakage behaviour. Consequently, SA involves the PDA simulations before the addition of artificial emitters for leakage modelling.

Micro calibration focuses on the fine-tuning of sensitive parameters to ensure a well-fitted model. There are three micro calibration approaches. The first is the manual (trial and error) approach. This is laborious and difficult to achieve high accuracy calibration (Savic & Walters 1995; Simukonda 2021). The second is the explicit (analytical) approach. Its major limitation is that it can only be conveniently applied to small networks. For large networks, its application requires a lot of network skeletonisation (Savic & Walters 1995). The third is the implicit (optimisation), which is the automated approach. It uses algorithms searching for parameter values to form the optimal solution. Objective functions comparing the measured and simulated results are used to evaluate the fitness of each solution. The drawbacks of such methods occur when the number of parameters increases and the demand for accuracy is high. The solutions become complex and the power of the search method diminishes (Savic & Walters 1995). This study applies the implicit approach using genetic algorithm (GA) due to its automated nature and suitability for complex WSSs optimisation (Farmani et al. 2007). The parameters to be calibrated during micro calibration of the hydraulic model with leakage are (i) consumption base demands, (ii) the emitter exponent for the consumption demand, (iii) the emitter coefficients for the leakage emitters, (iv) the emitter exponent for leakage, and (v) pipe roughness.

In Equation (6), the error is calculated from the known system input volume and the simulated values to ensure that the simulated inflow to the zone does not exceed the actual allocation. This is in line with the assertion that knowledge of system inflows can be used as a system constraint (Branisavljević et al. 2009). Minimum pressure and maximum unit head loss can also be used as constraints as these quantities have practical limitations. Constraints should consist of flexible margins to broaden the range of feasible optimisation solutions.

In Equation (7), is 1 if a single multiplication factor is used for the whole network. The value of N is 2 because one decision variable is for the consumption demand emitter exponent, and the other one is for the leakage emitter exponent. This arises from the treatment of the emitter exponent as a global value in the current version of EPANET. This means at any given instance; all consumption demand emitters in the network have the same value and all leakage emitters have their own same value. For the emitter coefficients, all the nodes in a DMA have the same coefficient value. Hence, the number of leakage emitter coefficients is equal to the number of DMAs.

In EPANET 2.2, for the PDA simulation option, the consumption demand exponent (called the pressure exponent) has a default value of 0.5. However, it varies between 0.5 and 2.5 (Mahmoud et al. 2017) and the actual value can be determined through calibration (Cheung et al. 2005). This range should be used to set the lower and upper bounds during optimisation. For the leakage emitter exponent, derived values from the field and experimental studies identify that its range varies from 0.36 to 2.5 (Fuchs-Hanusch et al. 2016). Thus, this range is used as the bound for the emitter exponent during optimisation. The EPANET-MATLAB-toolkit is utilised to link the hydraulic solver EPANET and available functions in MATLAB such as the GA.

Consequently, pressure analysis of the WSS is essential to use Equation (8) to determine ranges forming bounds for the leakage emitter coefficient for each DMA. The pressure analysis is conducted using a code written using the EPANET-MATLAB-toolkit. To get pressure values for each DMA, the first step is the splitting of the IWSS into DMAs based on water supply schedules. Due to different supply hours and DMA locations, the pressure in each DMA will vary. As leakage varies with pressure, emitter coefficient values vary for the different DMAs. The next step is to run a PDA simulation of the WSS consisting of demand nodes without artificial emitter elements. This simulation run aims to estimate the pressure values at each node. For each DMA, the minimum (P_min) and maximum (P_max) pressure values should be noted. If measured pressure data is available, it should be used instead.

The calibration data used in the current study was collected by (O'Donnell 2010) at different times over 4 months (between April 2009 to July 2009). The measurements were done as part of the consultancy on behalf of the Japan International Cooperation Agency (JICA). The measurement data was, and still is, the only systematic data on record for the LWSS and the Chelstone zone. Its quality was compromised because during 4 months, the hydraulic conditions most likely changed significantly (Simukonda 2021) and the data comprised only flow measurements. Considering the uncertainties arising from the compromised calibration data, the error threshold was relaxed from 10 (U. S. Environmental Protection Agency 2005) to 11%. Thus, if the simulated solution of the isolated zone had a sum of flows in the two pipes greater than of the simulated value for the whole LWSS model, a penalty function was applied to the RMSE calculation as shown in Equation (6).

The first step of the approach was the comparison of DDA and PDA simulation flows in high-pressure areas. This was meant to check if the ‘required pressure’ setting for the PDA simulation was correct. The deviation was 2.1%. This was an acceptable agreement between the analysis options results and no ‘required pressure’ parameter setting was modified.

The macro calibration was done by altering the valves along pipes that influenced the simulated flows. The original valve settings gave the MAE and RMSE between simulated and measured flow values of 0.906 and 1.756 L/s, respectively. Several changes in valve settings gave poorer MAE and RMSE values. This implied that the hydraulic model was already optimally macro calibrated from the rough tuning it underwent during earlier studies (Simukonda 2021).

For SA, there were no alterations to base demand values because of limited information and they were assumed to be correct (Bhave 1988). The model results were not sensitive to pipe roughness changes implying that the hydraulic model was adequately calibrated by the developers – BCHOD consulting engineering (O'Donnell 2010). The model was not sensitive to the demand emitter exponent values. Therefore, the default value of 0.5 was maintained for the consumption demand emitter exponent. Then, artificial leakage emitters (ALEs) were added and the default emitter exponent value of 0.5 was applied. Emitter coefficients were also given the value of 0.5 for this trial run. The RMSE and MAE produced were 45.18 and 9.42 L/s, respectively. These showed significant shifts from the macro calibration values before the addition of the emitters. This led to the identification of the leakage emitter coefficients and exponent as the most influential parameters and the only focus for micro calibration.

The maximum and minimum pressure for each DMA were obtained using the EPANET-MATLAB-toolkit. From the pressure values, the upper and lower emitter coefficient bounds were calculated using Equations (10) and (11). The average NRW values per node were determined from the NRW percentages from the LWSS survey data in O'Donnell (2010) which also contained the average nodal base demands. The data used to determine the emitter coefficient bounds for each DMA are summarised in Table 1. For the LWSS, a minimum pressure of 7 m was required at the demand nodes (O'Donnell 2010). This was the required pressure used in the hydraulic model under calibration. Out of 16, 6 DMAs could not provide adequate pressure during their scheduled water supply periods (Table 1). The table also shows the upper and lower emitter coefficient values for each DMA which were used during the GA optimisation process. The emitter coefficient values for some DMAs such as DMA 13 and DMA 15 are very large compared to others. This is because these DMAs are represented by single nodes and have high levels of leakage, while the other DMAs have multiple nodes into which the average DMA leakage was divided. Conversely, the very small emitter coefficient values are a result of the relatively low leakage values that were divided across DMAs with many demand nodes.

The constraints for the optimisation were based on the water inflow rates to the Chelstone zone. Two pipes that connect Chelstone to the other zones in the complete LWSS hydraulic model were investigated under PDA. The average flows for each of the two pipes were determined from the simulation results of the whole LWSS hydraulic model. The average pipe flows were 38.9 and 372.9 L/s. These were used as comparative values in the GA. Previous calibration work on the whole LWSS model was conducted using manual calibration. In the previous work (Simukonda 2021), the calibration parameters were the leakage emitter exponent (α) and the total base demand (Bd_Totalj). Using such a laborious method could not guarantee good calibration for the large hydraulic model involved. Moreover, the quality of the calibration was also affected by the compromised quality of the available calibration data used.

The MATLAB-provided single-objective optimisation GA was used, and its settings were varied by altering the crossover fraction, population size, and the number of generations. The default mutation function was used. The default crossover fraction setting in the MATLAB-provided GA is 0.8. Using the crossover fraction functionality in MATLAB, the crossover fraction was altered to 0.5, 0.7, 0.75, 0.8, and 0.9. To compare the optimal solutions produced using the different crossover values, the generation and population sizes were kept constant at 30 and 10, respectively. Each crossover fraction alternative was run five times to examine the variations in output for the same inputs. The results of these runs are summarised in Table 2. The table shows that the small generation and population sizes were able to give acceptable optimisation results.

The best runs using each crossover fraction are shown in bold numbers in Table 2. The first run, using a crossover fraction of 0.9, gave the lowest RMSE and MAE values of 4.696 and 0.979 L/s, respectively. Ayad et al. (2013) recommend that the range for the crossover fraction is between 0.7 and 0.9. Thus, this study's findings align with other studies in the field. Table 3 displays a summary of GA runs using 0.9 as the crossover fraction for different population and generation sizes. The table shows that a population size of 50 over 100 generations gave the lowest RMSE and MAE. It also shows how the RMSE and MAE results improved with increased population size and generations.

In Figure 7, the elements are the artificial consumption demand emitter, the artificial flow control valve, the artificial check valve 1 and 2 (ACV1 and ACV2), the artificial node, and the ALE.

The emitter exponent proposed for the optimal solution with a population size of 50 and 100 generations is 0.364. The global value for the emitter exponent can introduce an error as the dominant failure type of the pipes in DMAs can vary. The dominant pipe materials in the Chelstone zone network are asbestos cement and steel/galvanised iron. Emitter exponent values range from 0.36 to 0.7 for circular holes in steel pipes and from 0.5 to 0.55 for circumferential cracks in asbestos cement pipes (Fuchs-Hanusch et al. 2016). Consequently, the proposed emitter exponent is within the lower boundary for this parameter. The emitter coefficients for the various DMAs are shown in Table 4 (rounded off to 5 decimal places). Table 4 shows some coefficients that are very low while others are large. These differences stem from the assumption that all demand nodes in the network will have leakage. This implies that leakage is distributed across all demand nodes in the WSS. Consequently, DMAs with many nodes and relatively low leakage quantities have lower emitter coefficient values than DMAs that are represented by fewer demand nodes (due to network skeletonisation).

Using 0.364 as the global emitter exponent and the DMA emitter coefficients in Table 4, the optimal GA output was produced with a population size of 50 and a crossover fraction of 0.9 and 100 generations. The absolute errors of each simulated pipe flow compared to measured data are summarised in Table 5.

The large absolute error values are associated with pipes P-2864, P-6480, and P-2965. The common feature of these pipes is that they have relatively low flow values. A relatively small change in flow produces a greater relative change in absolute error. Other than these pipes, most simulated pipe flows fall within the 11% range of the measured value indicating that the micro calibration was achieved. This assertion is based on the relaxed error threshold from the guidelines for steady-state simulations contained in U. S. Environmental Protection Agency (2005). It is also based on the understanding that the quality of calibration results depends on the quality and comprehensiveness of the calibration data used. In this case, the quality of the data was compromised as already alluded. It is unlikely that the data collected over different system operation conditions can be accurately represented by simulations presuming that all the data were collected on the same day. Furthermore, the calibration data quality type of good, bad, and useless (Walski 2000) could not be ascertained as this requires knowledge of the data collection procedure. Moreover, comprehensiveness was compromised by the lack of measured pressure data.

A methodology for calibrating IWSS hydraulic models that incorporate NRW through leakage modelling was developed. Inclusion of leakage modelling was done by adding artificial emitters to demand nodes and the model calibration involved four major steps: checking the alignment of the EPANET 2.2 hydraulic simulation options results in high-pressure areas, performing macro calibration followed by SA, and fine-tuning the model through micro calibration using the MATLAB-provided single-objective GA. A penalty value was applied to make solutions outside the feasible range unfavourable. The methodology was applied to the real-world case study IWSS hydraulic model. The data used for micro calibration constituted only pipe flow measurements done by the original model developers. The calibration results gave low RMSE and MAE of 4.527 and 0.944 L/s, respectively. Moreover, the absolute errors for individual sets of measured and simulated pipe flows of the calibrated case study hydraulic model were mostly less than 10% which is the threshold for pipe flows. These statistics reflect good accuracy in simulated data especially since the quality and comprehensiveness of the measured calibration data were uncertain.

The major contributions of the study are the step-by-step approach to the calibration of IWSSs hydraulic models that include leakage modelling under data scarcity conditions, the determination of emitters' exponent and coefficients through optimisation thereby avoiding the trial-and-error procedure of determining the coefficients and the estimation of the lower and upper bounds of the emitters' coefficients based on each DMA average pressure and level of NRW. One of the limitations of the proposed calibration methodology is that it could not be validated due to the scarcity of the measured data. Thus, the micro calibration was deemed sufficient. Moreover, the high sensitivity of the calibration results to the values for the emitters' coefficient and exponent means the accuracy of the values of these parameters, which has a bearing on the convergence of the optimisation solutions, is not yet determined. To validate the proposed methodology, other hydraulic models with sufficient monitoring data for leakage, pressure, and flows retrieved under relatively constant boundary conditions should be used. With adequate high-quality measured data, calibration can be done using a specified level of accuracy for the coefficients and validation can be used to validate them.

The authors gratefully acknowledge the Commonwealth Scholarships Commission for sponsoring (ZMCA-2015-141) the second author's PhD studies.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.