Simulation of a Water Distribution Network with Key Performance Indicators for Spatio-Temporal Analysis and Operation of Highly Stressed Water Infrastructure

: An annual and lumped water balance assessment of a water distribution network is recommended by the International Water Association as a ﬁrst step and prerequisite for improving the performance of the network by minimizing real / physical water losses, burst incidents, water theft, nonrevenue water, and energy consumption, among others. The current work suggests a modeling approach for developing the water balance of a network spatio-temporarily, in hour time-scale and neighborhood granularity. It exploits already established key performance indicators and introduces some new ones to highlight the potential in improving the management of a water distribution network when having a detailed spatio-temporal supervision, especially when the spatial and temporal conditions are variable. The methodology is applied in a seasonally touristic and hilly case study. Additionally, a pressure management scheme is applied to further exploit the potential of such a toolkit. For the investigated case study, the town of Skiathos, the annual real losses are estimated equal to 50.9–52.2% of the system input volume, while apparent losses are estimated to be about 5.6–6.6%. Losses depict intense seasonal variability. Real losses range from 38.8–39.6% in summer months to 63.3–64.7% in winter months, while apparent losses range from 8.4–9.3% in summer to 1.3–2.5% in winter. Annual water theft is estimated to be at least 3.6% of system input volume. Spatial variability, which is linked to the elevation and the di ﬀ erent urban land uses is proven to play a signiﬁcant role in the neighborhoods’ water balances and various key performance indicators are suggested and applied for the pressure control scheme. The annual potential savings due to the applied scheme rise up to 51,300 m 3 for leakage and 53,730 m 3 for pressure-driven demand.


Introduction
The United Nations (UN) 2030 Agenda, which culminates in the sustainable development goals (SDGs), sets the preservation of freshwater at the core of sustainable development. Discussions on SDG 6, which refers to clean water and sanitation, focus on the threat of climate change to the availability of water, as well as on the need to manage, in an integrated way, water and energy and their interlinkages to climate change. The European Environmental Agency [1] states that water scarcity is a major threat for Europe, especially for the southern regions. Among the reasons that cause this stress and have triggered drought incidents are the spatial and temporal variability in demand and availability, the overall increasing demand due to the rise in urbanization and changing lifestyles, the unevenly At the same time, Liemberger and Wyatt [7] estimated annual global nonrevenue water (NRW) levels in urban water distribution networks (WDNs) as high as 126 billion m 3 or 39 billion USD, with 9.8 billion m 3 or 3.4 billion USD corresponding to Europe. The International Water Association (IWA) Water Loss Task Force recommends four primary measures for optimized real losses management, namely: (i) management of pipelines and assets, (ii) pressure management, (iii) speed and quality of repairs, and (iv) active leakage control. To enhance these measures, it sets at the core of such a multitargeted task a detailed component analysis of the WDN inflow into revenue and NRW, real and apparent losses, etc, [8].
Such an analysis of the annual system input volume (SIV) has increasingly been implemented for water supply systems (WSSs) around the world with variable precision, depending on the available data, measurements, knowledge of the WDN and available expertise. The analysis typically focuses on completing the iconic IWA water balance (WB) table of SIV components, which are defined by a taxonomy that distinguishes between revenue and nonrevenue water, authorized and unauthorized consumption, real and apparent losses, metering inaccuracies, billed and unbilled consumption, metered and unmetered consumption, and leakage in the transmission/distribution mains, service connections and tank overflows. The task has usually been implemented in an aggregated way for the whole WDN and for annual time step, although there have also been few approaches that focus on the estimation of leakage, rather than on the construction of the whole, or part of the IWA table in a spatio-temporal manner. The difficulty in constructing a precise and trustworthy WB table lies in the difficulty of assessing the real, apparent losses, and unmetered water use.
Significant effort has been made towards suggesting a reliable methodological approach for estimating, and/or modeling leakage (physical losses in mains and service connections), metering inaccuracies and unauthorized consumption (water theft). Almandoz et al. [9] proceeded to a real losses assessment by perceiving the apparent losses as a water consumption pattern dependent component and used four different demand profiles: domestic, industrial, commercial and official. Cabrera and Pellejero [10] and Ismail and Puad [11] conducted bottom-up assessments of physical losses and highlighted that estimating leakage through night flow can minimize the error, while night consumption is minimum and easier to determine. Giustolisi et al. [12] proposed an algorithm for the simultaneous simulation of pressure-driven demand (PDD) and leakage by fully integrating the two At the same time, Liemberger and Wyatt [7] estimated annual global nonrevenue water (NRW) levels in urban water distribution networks (WDNs) as high as 126 billion m 3 or 39 billion USD, with 9.8 billion m 3 or 3.4 billion USD corresponding to Europe. The International Water Association (IWA) Water Loss Task Force recommends four primary measures for optimized real losses management, namely: (i) management of pipelines and assets, (ii) pressure management, (iii) speed and quality of repairs, and (iv) active leakage control. To enhance these measures, it sets at the core of such a multitargeted task a detailed component analysis of the WDN inflow into revenue and NRW, real and apparent losses, etc, [8].
Such an analysis of the annual system input volume (SIV) has increasingly been implemented for water supply systems (WSSs) around the world with variable precision, depending on the available data, measurements, knowledge of the WDN and available expertise. The analysis typically focuses on completing the iconic IWA water balance (WB) table of SIV components, which are defined by a taxonomy that distinguishes between revenue and nonrevenue water, authorized and unauthorized consumption, real and apparent losses, metering inaccuracies, billed and unbilled consumption, metered and unmetered consumption, and leakage in the transmission/distribution mains, service connections and tank overflows. The task has usually been implemented in an aggregated way for the whole WDN and for annual time step, although there have also been few approaches that focus on the estimation of leakage, rather than on the construction of the whole, or part of the IWA table in a spatio-temporal manner. The difficulty in constructing a precise and trustworthy WB table lies in the difficulty of assessing the real, apparent losses, and unmetered water use.
Significant effort has been made towards suggesting a reliable methodological approach for estimating, and/or modeling leakage (physical losses in mains and service connections), metering inaccuracies and unauthorized consumption (water theft). Almandoz et al. [9] proceeded to a real losses assessment by perceiving the apparent losses as a water consumption pattern dependent component and used four different demand profiles: domestic, industrial, commercial and official. Cabrera and Pellejero [10] and Ismail and Puad [11] conducted bottom-up assessments of physical losses and highlighted that estimating leakage through night flow can minimize the error, while night consumption is minimum and easier to determine. Giustolisi et al. [12] proposed an algorithm for the simultaneous simulation of pressure-driven demand (PDD) and leakage by fully integrating the two demands in pipe level. Tabesh et al. [13] introduced a methodology for distinguishing the pressure-dependent components from the independent ones.
Puust et al. [14] distinguished the literature approaches of leakage assessment into top-down approaches where leakage is estimated through WB assessment, and the bottom-up approaches, where the leakage is estimated through the summing of its components. Kanakoudis and Tsitsifli [15] shared a thorough aggregate and annual component analysis for a series of case studies in WDNs in Greece following the established IWA methodology. They concluded that such an analysis is a useful tool for the local water utility, while a simple definition of SIV and billed water can be misleading for the performance of the network. Regarding the assessment of metering inaccuracies, they suggested laboratory investigation, while they distinguished under-registering and over-registering water-meters according to their age. Additionally, the practice of pressure management (PM) and division of the network in district metered areas (DMAs) occurred as a necessary approach for reducing real losses. Knakoudis and Tsitsifli [16] implemented the IWA WB assessment on a semiannual time scale, suggesting that a finer scale than the annual can be revealing for cases that depict seasonal demand peaks, while a few years later they proceeded to an even more detailed, bimonthly WB analysis [17]. Cobacho et al. [18] introduced a methodology for spatio-temporal simulation of leakage, using the EPANET software and the concept of simulating leakage through weighted leakage emitters. Sophocleus et al. [19] introduced a leakage localization methodology based on head pressure and flow measurements tested in two case studies in the U.K. thus reducing the repair time for leakage; according to the article, this could lead to water savings as high as 70% of total loss.
More macroscopic investigation has also been implemented in an attempt to address the need of integrated water management linking the urban water supply to other urban water systems, such as wastewater and stormwater. Makropoulos et al. [20] and Rozos and Makropoulos [21] extended the boundaries of the typical investigation focus on the urban water cycle, involving wastewater, urban stormwater, potable water, grey water, and green water. They also incorporated indicators from all sustainability pillars: environmental, economic, social, and technical, such as operational costs and energy. In their review on urban hydroinformatics, Makropoulos and Savic [22] presented a scheme according to which the modeling is facilitated by the increasingly fluent flows of data and information in the context the developments in information and communication technologies, cloud based information platforms, and remote monitoring. Such a technological landscape in combination to well-advanced methodological frameworks could lead to the next step in water distribution networks modeling, what has been referred as water distribution network digital twins [23].
The importance of assessing the IWA standard WB table in combination with employing a PM scheme and the need to enhance leakage localization and detection of unauthorized consumption raise the need for a methodology that would facilitate the spatio-temporal supervision and would provide the WDN operator with insight on the various components of SIV. The temporal analysis is essential, especially in cases of intense seasonality in the water demand, such as in the highly touristic regions. Seasonality also becomes relevant in the spatial distribution of land uses within the urban landscape, which in turn causes an intense seasonal shift in the demand map.
In this article, a toolkit for the spatiotemporal analysis and simulation of SIV, SIV components and critical key performance indicators (KPIs) throughout the WDN is presented to serve as supervision support for the network and its properties, stresses and potential of improvement, or the basis for a digital twin. The well-established IWA WB table is assessed in neighborhood granularity and hourly time step. The WB components assessed include system input volume, billed consumption, nonrevenue water, unauthorized consumption, real losses, and metering inaccuracies. All the components are identified locally and instantly offering a useful tool to a water utility for creating dynamic hotspot maps that can contribute to the optimization of the WDN management through: (i) the localization of leakage and water theft in neighborhood resolution and the identification of drivers of leakage; (ii) the assessment of metering inaccuracies; (iii) the quantification of tourism impact on water demand; and iv) the application of a PM scheme and the quantification of its beneficial effect. The assessment of the WB table is implemented by distinguishing the components in pressure-dependent and time-pattern-dependent ways and an application of a nested system of two loops, an inner and an outer, which, through iterative runs, close when the simulated pressure best fits the actual network pressure. Spatio-temporal critical KPIs are introduced to facilitate a detailed supervision of the WDN. The KPIs include the energy consumption in the water components (water-energy nexus), a pressure-driven demand indicator, and various expressions of water components such as "per-connection", "per-customer," and "per-network-length" indicators. The aforementioned KPIs are also used for the assessment of the PM scheme performance. The localization of the water components is used to conduct rough conclusions regarding the urban land uses, specifically the residential and the touristic ones.

Case Study Description
Skiathos is a small island municipality in Thessaly, Greece. The island has approximately 6000 registered inhabitants according to the 2011 census [24], while the water utility counts approximately 3500 consumers [25]. The island population shows intense seasonal variability due to high tourist influx. The tourist season, from April to September, has a high peak in August, which often exceeds 90,000 tourists for the whole month with an annual average ascending trend [25]. As expected, this is reflected in the water demand, which depicts seasonal oscillations in SIV, billings, nonrevenue, and nonrevenue proportioning. Indicatively, the August peak summer water withdrawals often exceed a 130% increase compared to the winter minimum withdrawals, which shows a link between the touristic activity and the water demand profiles. The water distribution system (WDS) includes a single drilling, a tank, and the WDN, which is currently under renovation. This research applies to the WDN of Skiathos before its renovation works had started, from 2011 to 2016. During these years the WDN is characterized by an extremely high rate of NRW gradually climbing up from 40%, during the last trimester of 2011, to 70 % of SIV, during the last trimesters of 2015 and 2016 ( Figure 2). The SIV values in Figure 2 from January 2011 to December 2016 are trimester sums of daily SIV readings. The SIV has an accuracy range of ±1% and the NRW has an accuracy range from ±1.5 to ±3.2% (95% confidence limits) depending on the billing period (Table 1). The WDN operates through a single-DMA scheme. The town is characterized by bold relief, with two hills. No PM scheme is applied at the WDN except for a switch of pressure every six months from the tourist to the nontourist period.
Water 2020, 12, x FOR PEER REVIEW 4 of 32 system of two loops, an inner and an outer, which, through iterative runs, close when the simulated pressure best fits the actual network pressure. Spatio-temporal critical KPIs are introduced to facilitate a detailed supervision of the WDN. The KPIs include the energy consumption in the water components (water-energy nexus), a pressure-driven demand indicator, and various expressions of water components such as "per-connection", "per-customer," and "per-network-length" indicators. The aforementioned KPIs are also used for the assessment of the PM scheme performance. The localization of the water components is used to conduct rough conclusions regarding the urban land uses, specifically the residential and the touristic ones.

Case Study Description
Skiathos is a small island municipality in Thessaly, Greece. The island has approximately 6000 registered inhabitants according to the 2011 census [24], while the water utility counts approximately 3500 consumers [25]. The island population shows intense seasonal variability due to high tourist influx. The tourist season, from April to September, has a high peak in August, which often exceeds 90,000 tourists for the whole month with an annual average ascending trend [25]. As expected, this is reflected in the water demand, which depicts seasonal oscillations in SIV, billings, nonrevenue, and nonrevenue proportioning. Indicatively, the August peak summer water withdrawals often exceed a 130% increase compared to the winter minimum withdrawals, which shows a link between the touristic activity and the water demand profiles. The water distribution system (WDS) includes a single drilling, a tank, and the WDN, which is currently under renovation. This research applies to the WDN of Skiathos before its renovation works had started, from 2011 to 2016. During these years the WDN is characterized by an extremely high rate of NRW gradually climbing up from 40%, during the last trimester of 2011, to 70 % of SIV, during the last trimesters of 2015 and 2016 ( Figure 2). The SIV values in Figure 2 from January 2011 to December 2016 are trimester sums of daily SIV readings. The SIV has an accuracy range of ±1% and the NRW has an accuracy range from ±1.5 to ±3.2% (95% confidence limits) depending on the billing period (Table 1). The WDN operates through a single-DMA scheme. The town is characterized by bold relief, with two hills. No PM scheme is applied at the WDN except for a switch of pressure every six months from the tourist to the nontourist period.   The deployment of a more efficient WDN in Skiathos has become an urgent need due to several reasons. The island has been facing seasonal water scarcity, due to the increased water demand driven by tourism. At least two drying-out incidents occurred over the last five years, in August, due to the extreme drop of the water table of the aquifer level at the drilling location. Due to these incidents, the water supply was interrupted in the island, until the aquifer recovered. Naturally, such incidents have multiple social, financial, and environmental impacts, as well as sanitation impacts, since water lies in the core of most everyday activities, especially in a tourist destination. On top of that, the overall aquifer level drop results in seawater intrusion and groundwater salinization, with detrimental consequences, such as the natural mineral mercury release into groundwater. Mercury concentration has repeatedly been recorded to exceed WHO limits, regarding toxicity, even by six times [26]. The water utility of Skiathos (DEYASK) has declared the water in the island as nonpotable. Overall, reducing water losses, optimizing water source withdrawals, investigating the potential of alternative water sources and applying innovative treatment methods are some of the tasks that are of utmost importance in Skiathos.

Infrastructure and Data
The following infrastructure installed in Skiathos is used for retrieving data needed for the simulation of the dynamic characteristics of the WDS: flowrates are being monitored at the inflow point of the WDN, while pressure is being monitored at three points of the WDN. The daily SIV values from January 2011 to December 2016 are readings of a mechanical meter with accuracy of ±1% for each value. The values of SIV are also checked often by the water utility volumetrically by cross-checking with the tank's volume. The values of 2016 in particular are additionally compared to the readings of an electromagnetic flowmeter with accuracy of ±0.25%. For year 2016, there is a pressure-reducing valve (PRV), a CSA XLC 410, DN150 PN16 in particular, that records pressure and flowrate (electromagnetic flowmeter) in the inflow of the WDN, after the tank, with a 15 min time step-currently not regulating pressure automatically-and three cello sensors that record pressure at three points around the network with a 15 min time step. The pressure data are used for the calibration Water 2020, 12, 1149 6 of 31 of the model, as the fitting of the simulated pressure values to the recorded ones is used as the criterion for the reliability of the simulation. Billing data were provided by DEYASK. There are 3500 water meters in the city used for billing purposes, with a billing period of three months, from January to March, April to June, July to September, and October to December, or the first, second, third and fourth period, respectively. The procedure of manually recording the water meters lasts for approximately ten days, while this procedure starts before and finishes after the end of the billing period, which means that there is no systematic hysteresis between the records and the consumption. No fixed minimum consumption is imposed on the billings, so the billings represent what is metered and billed. All water-meters are geolocated with a hand-held GPS device of 4.9-m accuracy, during multiple field trips with the guidance of the utility's staff and matched to corresponding billing data. In order to keep the privacy of water consumers, when these data are presented in public access, the scale is zoomed out to "neighborhood" level and the billed consumptions are summed.

Simulation of the WDN with EPANET
The first step of the methodology includes the simulation of the static and dynamic characteristics of the WDN, through EPANET software [27]. For this purpose, historical maps of the WDN are used to incorporate network geometry and pipe diameter, length, and material. The overview map in Figure 3 presents the main elements of the WDS. The single water supply drilling is found in the northwestern part of the city. Groundwater is pumped from the drilling and pumped up to the water tank located on the hill ( Figure 3). From the tank, water is distributed by gravity to the city area. The main pipeline of the WDN goes back down the hill to the SE direction and then to the E-NE around the city centre ( Figure 3). The pipe diameters are presented in Figure 3 in three clusters and the pipe materials, amiant, metal (cast-iron) and PVC, are presented in Figure 4.
Water 2020, 12, x FOR PEER REVIEW 6 of 32 reducing valve (PRV), a CSA XLC 410, DN150 PN16 in particular, that records pressure and flowrate (electromagnetic flowmeter) in the inflow of the WDN, after the tank, with a 15 min time stepcurrently not regulating pressure automatically-and three cello sensors that record pressure at three points around the network with a 15 min time step. The pressure data are used for the calibration of the model, as the fitting of the simulated pressure values to the recorded ones is used as the criterion for the reliability of the simulation. Billing data were provided by DEYASK. There are 3500 water meters in the city used for billing purposes, with a billing period of three months, from January to March, April to June, July to September, and October to December, or the first, second, third and fourth period, respectively. The procedure of manually recording the water meters lasts for approximately ten days, while this procedure starts before and finishes after the end of the billing period, which means that there is no systematic hysteresis between the records and the consumption.
No fixed minimum consumption is imposed on the billings, so the billings represent what is metered and billed. All water-meters are geolocated with a hand-held GPS device of 4.9-m accuracy, during multiple field trips with the guidance of the utility's staff and matched to corresponding billing data. In order to keep the privacy of water consumers, when these data are presented in public access, the scale is zoomed out to "neighborhood" level and the billed consumptions are summed.

Simulation of the WDN with EPANET
The first step of the methodology includes the simulation of the static and dynamic characteristics of the WDN, through EPANET software [27]. For this purpose, historical maps of the WDN are used to incorporate network geometry and pipe diameter, length, and material. The overview map in Figure 3 presents the main elements of the WDS. The single water supply drilling is found in the northwestern part of the city. Groundwater is pumped from the drilling and pumped up to the water tank located on the hill ( Figure 3). From the tank, water is distributed by gravity to the city area. The main pipeline of the WDN goes back down the hill to the SE direction and then to the E-NE around the city centre ( Figure 3). The pipe diameters are presented in Figure 3 in three clusters and the pipe materials, amiant, metal (cast-iron) and PVC, are presented in Figure 4.

Water Consumption, Leakage and Other SIV Components Spatio-Temporal Distribution Simulation
The SIV components ( Figures 5 and 6), as defined and established by IWA Blue Pages [28] and Manuals of Best Practice [29], are simulated in EPANET either as demand driven or as a pressure driven, according to the distinction made by Cobacho et al. [18]. They decoupled the leakage from the base demand by introducing nodal emitters-open valves to the atmosphere-to better mimic the pressure-dependent behavior of leakage, as it is described by Germanopoulos [30] and Germanopoulos and Jowitt [31]. Cobacho et al. [18] tested their methodology on five variable normalized demand time patterns of different ranges: pattern A and B with a range of 0 to 2.5 and two peaks, pattern C with a range of 0 to 1.5 with two peaks, pattern D with a range of 0 to 2.5 and one peak and pattern E with a range of 0 to 4 and two peaks. The study tests a case of 3500 m 3 /d. Although the recommended methodology does not set any restrictions regarding the demand timepatterns type or the volume of the demand, Skiathos case study demand is quite close to the daily volume demand tested (as shown in Figure 7) as well as to the tested pattern C (as shown in Figure  8).
The components that are simulated as demand driven are the billed authorized consumption (BAC) and the apparent losses, while real losses are simulated as pressure driven ( Figure 5). This does not mean that the components which are simulated as demand driven are not also partially pressure dependent. For the Skiathos case study, BAC consists only of metered BACs. Unbilled authorized consumption (UAC) is considered as negligible. The main demand for UAC would be that of the fire department for fire-fighting; however during the simulation period, no significant fire occurred in the island. Apparent losses, which consist of unauthorized consumption and customer metering inaccuracies, along with reading and handling errors, cannot be perceived as negligible, since water theft in Skiathos has been repeatedly noticed and the recording of water meters-that are older than 20 years-is implemented manually. These last two components are uniformly distributed to the WDN model following the time-pattern of the demand simulated as demand driven. Regarding real losses, "leakage on transmission and/or distribution mains" and "leakage on service connections up to the point of customer metering" are expected to constitute the largest part of NRW, since the WDN during the simulation period is characterized as a rather aging network in a bad condition overall-the construction of the WDN was initiated in the mid-1960s and evolved in the following decades. These two components are simulated as pressure-driven flows. Leakage and overflows at utility's storage tanks and operational losses are considered negligible and occurs before the SIV metering in the specific case study. There is no adequate information on leakage after the point of

Water Consumption, Leakage and Other SIV Components Spatio-Temporal Distribution Simulation
The SIV components ( Figures 5 and 6), as defined and established by IWA Blue Pages [28] and Manuals of Best Practice [29], are simulated in EPANET either as demand driven or as a pressure driven, according to the distinction made by Cobacho et al. [18]. They decoupled the leakage from the base demand by introducing nodal emitters-open valves to the atmosphere-to better mimic the pressure-dependent behavior of leakage, as it is described by Germanopoulos [30] and Germanopoulos and Jowitt [31]. Cobacho et al. [18] tested their methodology on five variable normalized demand time patterns of different ranges: pattern A and B with a range of 0 to 2.5 and two peaks, pattern C with a range of 0 to 1.5 with two peaks, pattern D with a range of 0 to 2.5 and one peak and pattern E with a range of 0 to 4 and two peaks. The study tests a case of 3500 m 3 /d. Although the recommended methodology does not set any restrictions regarding the demand time-patterns type or the volume of the demand, Skiathos case study demand is quite close to the daily volume demand tested (as shown in Figure 7) as well as to the tested pattern C (as shown in Figure 8).
The components that are simulated as demand driven are the billed authorized consumption (BAC) and the apparent losses, while real losses are simulated as pressure driven ( Figure 5). This does not mean that the components which are simulated as demand driven are not also partially pressure dependent. For the Skiathos case study, BAC consists only of metered BACs. Unbilled authorized consumption (UAC) is considered as negligible. The main demand for UAC would be that of the fire department for fire-fighting; however during the simulation period, no significant fire occurred in the island. Apparent losses, which consist of unauthorized consumption and customer metering inaccuracies, along with reading and handling errors, cannot be perceived as negligible, since water theft in Skiathos has been repeatedly noticed and the recording of water meters-that are older than 20 years-is implemented manually. These last two components are uniformly distributed to the WDN model following the time-pattern of the demand simulated as demand driven. Regarding real losses, "leakage on transmission and/or distribution mains" and "leakage on service connections up to the point of customer metering" are expected to constitute the largest part of NRW, since the WDN during the simulation period is characterized as a rather aging network in a bad condition overall-the construction of the WDN was initiated in the mid-1960s and evolved in the following decades. These two components are simulated as pressure-driven flows. Leakage and overflows at utility's storage tanks and operational losses are considered negligible and occurs before the SIV metering in the specific case study. There is no adequate information on leakage after the point of customer metering, which is not considered to be negligible on the one hand and is pressure-driven on the other. Due to this lack of information, this hidden component in BAC is not decoupled by the demand driven components. However, the private property and infrastructure in Skiathos is constantly renovated to meet touristic demands; thus, leakage after the water meter is expected to be much less than the corresponding one before the water meter.
Another specification in the Skiathos case study is the fact that a lot of houses and small enterprises have small tanks in their roofs. These tanks are not used in a regular basis, but only if there is a crisis, i.e., when the utility is disrupting the supply. There are only few houses (less than twenty) at the top of the two hills that use the tanks on a more regular basis, when, due to very high demand during the summer peak, the pressure is very low. These tanks are neither considered to significantly influence the overall consumption pattern nor the water metering in the three-month period of any customer.
Water 2020, 12, x FOR PEER REVIEW 8 of 32 customer metering, which is not considered to be negligible on the one hand and is pressure-driven on the other. Due to this lack of information, this hidden component in BAC is not decoupled by the demand driven components. However, the private property and infrastructure in Skiathos is constantly renovated to meet touristic demands; thus, leakage after the water meter is expected to be much less than the corresponding one before the water meter. Another specification in the Skiathos case study is the fact that a lot of houses and small enterprises have small tanks in their roofs. These tanks are not used in a regular basis, but only if there is a crisis, i.e., when the utility is disrupting the supply. There are only few houses (less than twenty) at the top of the two hills that use the tanks on a more regular basis, when, due to very high demand during the summer peak, the pressure is very low. These tanks are neither considered to significantly influence the overall consumption pattern nor the water metering in the three-month period of any customer.

Time-Pattern Spatio-Temporal Demands
The spatial distribution of BAC is implemented by linking the billing time-series of each water customer with the location of the respective water meter in the Skiathos WDN. Firstly, each customer code that appears in the billings is assigned to the corresponding water meter code, after the anonymization of the billings. Each water meter code is spatio-located with a GPS device. All locations are later imprinted in AutoCad with the water meter and customer codes assigned to them as attributes. Billed water that is not assigned to a water meter (due to customers moving, or changing SIV components as specified for the Skiathos case study and the respective simulation approaches. The spatial distribution of BAC is implemented by linking the billing time-series of each water customer with the location of the respective water meter in the Skiathos WDN. Firstly, each customer code that appears in the billings is assigned to the corresponding water meter code, after the anonymization of the billings. Each water meter code is spatio-located with a GPS device. All locations are later imprinted in AutoCad with the water meter and customer codes assigned to them as attributes. Billed water that is not assigned to a water meter (due to customers moving, or changing services) is introduced to the algorithm as "bulk" demand, which means demand that is not spatio-located but is distributed around the spatio-located water-meters in a weighted proportioning, according to the water-meters' assigned consumptions. According to this distribution of the "bulk" demand, the water meters that showed higher consumptions will be assigned bigger parts of the "bulk" demand and water meters with lower consumptions will be assigned smaller parts of the "bulk" demand in an analogous manner.
On the next step, the 3500 water-meters are grouped into 121 neighborhoods, with an average of 26.6 water-meters/neighborhood. The rest of the billed water-meters are treated as "bulk" demand. The neighborhoods are designed, in a way that clusters service connections that are close to each other. These WDN neighborhoods are referred to as landzones in this article. Each landzone is assigned to a specific demand geonode, a node on the WDN that is considered to supply all the water-meters of the landzone. The geonodes are also simulated in EPANET as demand nodes. The output of these steps is the spatial distribution of three-month-step BAC time-series.
In order to temporarily disaggregate BAC maps from a three-month time-step to its daily values, the SIV daily time-series patterns are followed with the assumption that each landzone follows the same SIV pattern. In Figure 7, the intense seasonal variability can be seen, as well as an ascending annually average SIV trend, also apparent in Figure 2. This procedure is described in detail in Kofinas et al. [32], where the three-month BAC values were simply disaggregated to daily, following the SIV pattern and the aggregate SIV values were disaggregated to the landzones.
The daily BAC is not perceived to follow a uniform distribution within the day, but it is expected to follow hourly patterns that are formed by averages of PRV SIV recordings, after they are modified, by subtracting the leakage, which is estimated using the bottom-up approach of minimum night flows combined to an iterative EPANET based process for the production of average WDN profiles. The description of leakage estimation is described in more detail in Section 2.5. PRV recordings are taken every 15 minutes, however it is decided that hourly time patterns are sufficiently detailed. The seasonality of the island touristic activity alters the daily SIV profiles, since rush hour times and the ratio of day-over night-consumption are driven and co-shaped by touristic water demand. In Figure 8, the twelve hourly SIV profiles are shown, as estimated by averaging PRV recordings, one for each month of 2016. It is apparent that for April to June and July to September billing periods, the relative night to day consumption is lower that the respective January to March and October to December. This does not imply that absolute summer night consumptions are lower than winter; on the contrary, they are much higher as shown in Figure 8. seasonality of the island touristic activity alters the daily SIV profiles, since rush hour times and the ratio of day-over night-consumption are driven and co-shaped by touristic water demand. In Figure  8, the twelve hourly SIV profiles are shown, as estimated by averaging PRV recordings, one for each month of 2016. It is apparent that for April to June and July to September billing periods, the relative night to day consumption is lower that the respective January to March and October to December. This does not imply that absolute summer night consumptions are lower than winter; on the contrary, they are much higher as shown in Figure 8.

Pressure-Dependent Spatio-Temporal Demands
Leakage is spatially distributed as pressure-dependent demand, through the introduction of emitters in every node. This means that every node is assigned with two demands, a time-pattern demand (demand-driven components) as described in Section 2.4.1 and a pressure-dependent flow (pressure-driven components) (Figure 9). The procedure followed by Cobacho et al. [18] is adopted. The emitters are described by a leakage coefficient or emitter coefficient Kj. The flowrate of the leakage demand is given by Equation (1).

Pressure-Dependent Spatio-Temporal Demands
Leakage is spatially distributed as pressure-dependent demand, through the introduction of emitters in every node. This means that every node is assigned with two demands, a time-pattern demand (demand-driven components) as described in Section 2.4.1 and a pressure-dependent flow (pressure-driven components) (Figure 9). The procedure followed by Cobacho et al. [18] is adopted. The emitters are described by a leakage coefficient or emitter coefficient K j . The flowrate of the leakage demand is given by Equation (1).
Leakage is spatially distributed as pressure-dependent demand, through the introduction of emitters in every node. This means that every node is assigned with two demands, a time-pattern demand (demand-driven components) as described in Section 2.4.1 and a pressure-dependent flow (pressure-driven components) (Figure 9). The procedure followed by Cobacho et al. [18] is adopted. The emitters are described by a leakage coefficient or emitter coefficient Kj. The flowrate of the leakage demand is given by Equation (1).

= ×
(1) where Kj is the emitter coefficient at node j, as defined by Equation (2), Pj is the pressure at node j, Figure 9. Scheme of the two types of node demands [18].
where K j is the emitter coefficient at node j, as defined by Equation (2), Pj is the pressure at node j, and N is the pressure exponent, as defined in Table 2 equal to the weighted average value of the pressure exponents of different materials. N = 0.64 × 1.13 + 0.06 × 1.41 + 0.3 × 0.91 = 1.08. Lambert [33] suggests that N is taken equal to 1 in the absence of more information regarding the distribution of materials. However, due to the unbalanced distribution of materials (PVC: cast iron: amiant = 0.64: 0.06:0.3), the pressure exponent's value is estimated in respect of this issue, while the diameters and the age of the pipes depict small variation and balanced distribution.
where K WDN is the WDN leakage coefficient, as defined by Equation (3), and Γ j is the normalized leak variable for node j, as defined by Equation (4) where Q WDN,real is the leakage of the whole WDN, and P WDN is the average WDN pressure, which is defined as described in Section 2.6.
where Γ j, service connections is the normalized leak variable of node j due to the variability of the number of service connections assigned to each node, through the respective landzone (Equation (5)). and Γ j, pipes length is the normalized leak variable of node j due to the variability of the total pipes' length assigned to its node (Equation (6)).
For each node and respective landzone, Equations (5) and (6) are firstly used for the estimation of the two specific Γ j s and then the final Γ j is estimated as their average with Equation (4).
Γ j, service connections = j, service connections N WDN, service connections (5) where N j, service connections is the number of service connections assigned to node j and N WDN,service_connections is the total number of service connections of the WDN Γ j, pipes length = L j L WDN (6) where L j is the length of the pipes that are connected to node j, taken that if a pipe is connected to two nodes, its length is equally distributed to the two nodes, and L WDN is the total length of the WDN pipelines. It is noted that a uniform ratio of 1.2 water-meters per service connection is assumed throughout the whole network as recorded during the field visit, which allows the ratio of Equation (5) to be estimated through the corresponding water-meter number.

Leakage Estimation and SIV Component Analysis
Through the sequence of Equations (1) through (6), the leakage is distributed to nodes according to the nodes' relative importance regarding leakage, as it is expressed by service connections, pipe length partitioning and node pressure. The Q WDN,real initial value is assumed to be equal to the NRW, which includes real losses and apparent losses. This is an initial value hypothesis for a first estimation of the WDN average pressure through EPANET. Once real and apparent losses are better approximated, following the night flow methodology (described in more detail in Section 2.5.2), the Apparent Losses component is transferred to the demand driven components (even though they are partially demand driven as well) and added to the distributed BAC ( Figure 10). This is achieved through two nested loops: an "inner" loop that iterates until convergence, at which time it triggers the "outer" loop that iterates and converges. Figure 11 presents the inner loop. partially demand driven as well) and added to the distributed BAC ( Figure 10). This is achieved through two nested loops: an "inner" loop that iterates until convergence, at which time it triggers the "outer" loop that iterates and converges. Figure 11 presents the inner loop. EPANET initially produces the QWDN,real amount, which is initially equal to the whole nonrevenue component and, as it converges, eventually to the real losses component, following the inner loop shown in Figure 11. This loop ends, when simulated and input leakage converge with an acceptable error (ε) of 0.001. The input leakage is estimated through the bottom up approach of the minimum night flow for each iteration. The average WDN pressure for the initial step can be a random value (preferably a value that is lower than the PRV inlet set pressure).
For decoupling real losses from apparent losses, the night flow leakage estimation approach is implemented, through the outer loop. This estimation assumes that night flow consists of leakage EPANET initially produces the Q WDN,real amount, which is initially equal to the whole nonrevenue component and, as it converges, eventually to the real losses component, following the inner loop shown in Figure 11. This loop ends, when simulated and input leakage converge with an acceptable error (ε) of 0.001. The input leakage is estimated through the bottom up approach of the minimum night flow for each iteration. The average WDN pressure for the initial step can be a random value (preferably a value that is lower than the PRV inlet set pressure).
For decoupling real losses from apparent losses, the night flow leakage estimation approach is implemented, through the outer loop. This estimation assumes that night flow consists of leakage and night consumption. For this purpose, 2019 average SIV quarterly flowrate profiles are constructed for every month of the year, while the minimum value of flowrate recorded for every 15 minutes of the day throughout each month is also recorded. The minimum average flowrate of all months' night flowrates gives the leakage flowrate for the specific WDN pressure, after the estimated night consumption is subtracted. Leakage profiles for every month are constructed based on WDN pressure profiles (an output of the inner loop). The new, decoupled from apparent losses, leakage is then input to the inner loop, while apparent losses are added to the time-pattern dependent demands, to give new WDN pressure profiles. The process iterates until the average actual recorded pressure values at the three points (variable in location and altitude) converge to the simulated ones. This closes both nested loops. Figure 12 presents the whole outer loop and its interactions with the inner loop. It is evident that increasing the number of pressure meters at variable points of the WDN improves the accuracy of the method.   Figure 11. Scheme for the inner loop of the methodology as adapted from the corresponding Cobacho, et al. [18] loop.

Daily Pressure Patterns
The leakage function is inextricably linked to the pressure dynamics. Estimating leakage through the night flow method implies that we know the pressure level at the time when the minimum night flow (MNF) occurs, while creating the daily leakage profile based on the MNF implies that we know the corresponding pressure pattern. With this in mind, twelve daily pressure profiles are created, one for each month. Skiathos intense touristic activity results in a rather intense shift of the residents' consumption patterns, since the composition of the consumer body changes seasonally from winter (permanent) residents to summer tourists. Except for the inflow of tourists, the change in weather causes changes in water consumption routines. The above cause in turn a seasonal change in the average daily water consumption profiles and the average daily pressure profiles. For the construction of the daily pressure profiles, the following steps are implemented ( Figure 12): (a) The recorded pressure time series by the three cellos were processed; outliers are removed and replaced with averaged values of the adjacent time slot records, and the quarterly time step is transformed to hourly time step by averaging the four records of every hour. The recorded data use refers to a time period from February 2015 to March 2017. (b) From the above profiles, average daily pressure profiles are constructed for each one of the three locations, by averaging horizontally the records that correspond to the same hour of the day for the whole month. Days of the same month but different year (e.g., all July days, independent of the year) are used for the construction of the respective month profile. (c) The average WDN pressure value of the monthly period is estimated through an EPANET initial run. This initial run holds the assumption of zero leakage, so the WDN pressure is expected to be overestimated. This is corrected through an iterative process, after the first estimation of leakage. (d) The cello pressure profiles are moved laterally in order to meet the average WDN pressure as estimated in the previous step.
(e) The three produced WDN pressure profiles are averaged, so that a single WDN daily pressure profile is created for every month.
It should be noted that all profiles obtained here are only indicative and are used as initial values for the process.

Estimation of Leakage Using the Night-Flow Approach
For the estimation of leakage, the night-flow approach is implemented; this would imply the use of the MNF of SIV recorded to estimate night leakage through the subtraction of an estimate of the night consumption [34]. MNF = Real losses + Minimum Night Use (7) where MNU= household use + commercial use + special use + after water meter burst [35] (a) For this estimation, the following steps are implemented: where 49.6 m 3 /hr is the minimum of averaged recorded flowrates and it occurs in January 4.15 am time slot (b) The formula L 1 L 2 = P 1 P 2 N 1 i is used to construct the leakage daily profiles (as suggested in [34]) according to the pressure profiles produced as described in Section 2.5.2, the night leakage and the respective night leakage pressure value. N 1 is estimated equal to 1.08. (c) For the twelve produced leakage profiles (Figure 13), the leakage estimated is compared to the minimum recorded night flow. The estimated amount should be less or at most equal, assuming that there might be a quarter of night hours, especially in winter time that for such a small village the consumption is negligible, if not zero.  (d) The leakage is integrated for every one of the 12 months and a percentage of leakage over SIV is estimated. (e) The difference of NRW percentage and leakage percentage is considered as the decoupled apparent losses percentage at this iteration and is added to the billed consumption of each EPANET node for the next EPANET run iteration.
The produced monthly diagrams of PRV recordings and leakage of Figure 13 (which are the final ones, after the end of the iterative process) show the intense seasonality of the case study. The winter months show relatively low SIV (≈90m 3 /sec) during the day and an extended minimum night SIV zone (from 60 m 3 /sec to 70 m 3 /sec) which lasts from 00:00 to 06:00. The level of the leakage during the winter months displays small variability (≈45 m 3 /sec). During the summer months the overall SIV gets higher. The summer profile has two peaks, a morning one around 09:00 and an evening one around 19:00. The evening one is the highest (≈180 m 3 /sec). The summer months do not depict an extended night minimum SIV zone, but a nadir which occurs around 04:00. Summer leakage level is higher in absolute values (≈50 m 3 /sec), while it gets more variable during the day in August (from 30 m 3 /sec to 50 m 3 /sec). The ratio max SIV/min SIV is higher in the summer (≈2) than the winter months (≈1.5), which is indicative of the more intense variability of water demand linked to the touristic activities.

Pressure Management Scheme and Key Indicators for the Performance of the WDN
The utility and value of such a tool, among many other aspects that will be analyzed in the results and discussion chapter, can also lie within the assessment of alternative PM schemes or even a scenario on dividing the network into DMAs. To this purpose, a number of established KPIs such as the water losses per connection, the water losses per network length, the absolute leakage, and the electrical energy costs have been used [36][37][38][39] and some original KPIs such as the absolute leakage difference between two scenarios, the PDD reduction (BAC or BAC + apparent losses) and the sum of leakage and energy cost are introduced and estimated spatio-temporally or in an aggregate value depending on the nature of the KPI.
The absolute leakage and the absolute leakage difference between two scenarios in m 3 would be the most common and useful KPI that the water utility manager would be interested to know. The amount of real losses saved due to an application of a different managerial scenario corresponds to water that is saved for the water resources. It is a variable that can easily be translated into groundwater level difference and help assess the issues of seawater intrusion, due to aquifer level drop.
Leakage per network mains length and the respective difference between two scenarios, a KPI suggested by the IWA best practice as indicative for assessing the performance of a management scheme [40] is also estimated.
An aspect that becomes increasingly interesting regarding the management of resources in general and specifically the management of water in an urban WSS is the resource nexus, which means the interlinkage of resources though a complex system of interrelations that may lead to synergies or trade-offs [41]. The water-energy nexus within the WSS can be found in the impact of leakage on head loss throughout the network [42] and the coupling of energy consumption and water withdrawals, since every cubic meter of SIV has used energy to be pumped, transferred and/or treated. The second energy component is taken into account and estimated for this analysis, according to the user requirements. This way, the two commodities can be managed synergistically in such a system, since saving water, almost proportionally, would save energy. The KPI that is used to express the water-energy relation is defined as the energy saved due to leakage reduction comparing various managerial schemes. The actual value of energy that is spent for a cubic meter of water in kWhr/m 3 is estimated by the actual energy consumption data of DEYASK and is perceived as a flat rate equal to 0.423 kWh/m 3 . This assumption needs to be reconsidered if the variability in short term pumping depths were more significantly influencing the pumping energy.
An insightful KPI for the performance of a PM scheme is the reduction of PDD. This can be estimated for a geonode based on the actual pressure of the node, the PM scheme pressure of the node, the actual consumption and the pressure exponent [43]. Two kinds of PDD volumes can be estimated: the reduction of the BAC and the reduction of the BAC plus the apparent losses, since they are expected to reduce as well even if this is not depicted to revenue. For the purpose of this work, the overall reduction of both components is estimated following the following formula.
In a similar manner, a KPI regarding the economic impact of leakage reduction is introduced. This KPI is expressed in euros and represents the cost of the leakage difference between two scenarios as this is priced by the utility, also including the respective price of energy consumed. For Skiathos, this cost is estimated as 0.112 euros/m 3 Other KPIs that are used to assess the performance of a PM scheme are relevant to the pressure fluctuations, since they are considered to relate to the networks stress and, eventually, bursts. These KPIs are the decrease in pressure fluctuation, the decrease in pressure at nights and the decrease in the number of nodes where a pre-set maximum pressure value is exceeded.
All these metrics are applied for a tested scenario of PM (compared to the actual operation), in which a virtual CP, a different one for every time step, is set to never fall under the value of 5 m. This pressure level is not respectful to the legislative framework in Greece, which requires a pressure value at the ground of at least 4 m*(number of building's floors + 1) for any building but is quite realistic according to the utility. Nevertheless, the purpose of this task is not to suggest a valid PM, but to make a proof of concept regarding the usability of the suggested tool to assess scenario performance. Figure 14 and Table 3 present the fitting of the simulated pressure at three points of the WDN to the actual recorded pressure as measured by pressure sensors for hourly values of the year 2016. The fitting is adequate, with R 2 values equal to 0.6729 (moderate fitting) for the central point, 0.7104 (strong fitting) for the eastern point, and 0.8852 (strong fitting) for the western point [44].

Results and Discussion
In Figures 15 and 16, SIV, BAC, apparent losses, real losses, revenue, and NRW are calculated for each trimester of 2016 and for the whole year, using the developed model. One of the results that emerge is that in the warm trimesters (April to September), the SIV increases by approximately 50% compared to the winter and fall periods. The BAC and revenue water percentages are higher during these trimesters (44-52% of SIV), while the respective percentages for the winter and fall decrease to the levels of 32-34%. Real losses fluctuate from the summer level of 38.8-39.6 % of SIV to the high winter level of 63.3-64.7 % of SIV. Apparent losses fluctuate from 1.3-2.5% (first trimester) to 8.4-9.3% (third trimester) with the highest values noticed from July to December. NRW is inversely proportional to SIV, thus the former is at lowest levels during the summer trimester when SIV is the highest. Regarding the overall annual analysis (Figures 16 and 17), BAC is estimated equal to 42.1-42.3% of SIV, real losses equal to 50.9-52.2%, theft higher than 3.6% and metering accuracies lower than 2.4% (Apparent losses equal to 6%).
Regarding the uncertainties of the measured and estimated values of the WB components, the error of SIV is taken equal to ±0.25% according to the flowmeter specifications, the error of BAC and revenue water is taken equal to 0, since they represent the exact amount of water that is billed, and the real losses error is estimated following the rules of error propagation for the equation that produced their values, as follows: where Error (L 1 ) is the error of real losses, Error (P) = Error (P1) = Error (P2) is the error of the simulated pressure values as estimated for each trimester and for the whole year of simulation, Error (L 2 ) = ±0.25% is the error of the minimum night flow measurement (the night use is a relatively small amount that does not change the error of the measurement even with a ±5% error), and N1 = 1.08.
The highest error estimated for the real losses equals to 1.55% and corresponds to the fourth trimester. The error of NRW is estimated, with 95% confidence limits, from ±0.37% to ±0.52% for the trimester values and equal to ±0.43% for the annual value. The error of the apparent losses is the most significant one with highest value equal to ±30.3% for the first trimester, while it is lower, equal to ±9.46% for the annual WB.
From the above it can be concluded that higher accuracy of the WB assessment can be achieved by increasing the accuracy of the model. The methodology followed for developing the simulation of the WDN uses average values in some specific steps, such as the monthly averaged consumption profiles and the monthly averaged pressure profiles. These average values are used as initial values at the beginning of its iteration, for disaggregation purposes and for the estimation of leakage variation during the day. The error in this specific application is low and leads to acceptable fitting, however if the N 1 value was higher and not close to 1, it would introduce a more significant error to the simulation. For this reason, it is suggested that the methodology developed can perform better with real time data and a denser grid of monitoring devices. This way the tool can evolve to a real digital twin.
Assuming a linear relation between metering inaccuracies and the BAC (the water meters in Skiathos are older than 10 years, thus expected to be under-registering flow with a rather linear relation [15,45,46]), theft can be estimated to take its higher values during the third and fourth trimester, higher than 6% and 5.2% of SIV respectively ( Figure 17). This implies that theft cannot be linked wholly to tourist activity but may be influenced by other social and/or economic factors such as the demand patterns of the local gardens and orchards. The metering inaccuracies are smaller than 3% (or even smaller) in all trimesters, and smaller than 2.4% annually. For this analysis it should not be forgotten that apparent losses as a whole have a relatively high level of uncertainty.
Water 2020, 12, x FOR PEER REVIEW 21 of 32 with real time data and a denser grid of monitoring devices. This way the tool can evolve to a real digital twin. Assuming a linear relation between metering inaccuracies and the BAC (the water meters in Skiathos are older than 10 years, thus expected to be under-registering flow with a rather linear relation [15,45,46]), theft can be estimated to take its higher values during the third and fourth trimester, higher than 6% and 5.2% of SIV respectively ( Figure 17). This implies that theft cannot be linked wholly to tourist activity but may be influenced by other social and/or economic factors such as the demand patterns of the local gardens and orchards. The metering inaccuracies are smaller than 3% (or even smaller) in all trimesters, and smaller than 2.4% annually. For this analysis it should not be forgotten that apparent losses as a whole have a relatively high level of uncertainty.        Except for the percentage of SIV that makes up real losses, the infrastructure leakage index (ILI) is calculated and is found equal to 8.66-larger than 8, which is the limit between technical performance categories C and D. This would set the Skiathos WDN technical performance to the Category D, which implies inefficient use of resources, indicative of poor maintenance and system condition in general [47]. Table 4 presents the parameters used for the estimation of ILI. A validation of the minimum night-flow analysis leakage assessment can be also made through the component analysis of real losses, as suggested by A. Lambert [48], who claimed that leakage assessment should preferably be implemented with more than one methodology to avoid errors. The component analysis of real losses is based on the assessment of the unavoidable annual real losses (UARL) due to leakage in mains and service connections. The component-analysis approach, based on the details of Table 4, gives a leakage of 45.57 m 3 /hr (=0.5 m 3 /d/connection × 187 connections /24 hr), taken into account the suggested rough value of 0.5 m 3 /d/connection for technical performance category D (ILI very close to the benchmark of 8 as estimated in Table 4), developed countries, and average WDN pressure of 50 m [47]. This value is quite close to the minimum night-flow approach estimation that equals 44.99 m 3 /hr.  Except for the percentage of SIV that makes up real losses, the infrastructure leakage index (ILI) is calculated and is found equal to 8.66-larger than 8, which is the limit between technical performance categories C and D. This would set the Skiathos WDN technical performance to the Category D, which implies inefficient use of resources, indicative of poor maintenance and system condition in general [47]. Table 4 presents the parameters used for the estimation of ILI. A validation of the minimum night-flow analysis leakage assessment can be also made through the component analysis of real losses, as suggested by A. Lambert [48], who claimed that leakage assessment should preferably be implemented with more than one methodology to avoid errors. The component analysis of real losses is based on the assessment of the unavoidable annual real losses (UARL) due to leakage in mains and service connections. The component-analysis approach, based on the details of Table 4, gives a leakage of 45.57 m 3 /hr (=0.5 m 3 /d/connection × 187 connections /24 hr), taken into account the suggested rough value of 0.5 m 3 /d/connection for technical performance category D (ILI very close to the benchmark of 8 as estimated in Table 4), developed countries, and average WDN pressure of 50 m [47]. This value is quite close to the minimum night-flow approach estimation that equals 44.99 m 3 /hr. A major function of the produced tool is the spatio-temporal supervision of the WDN regarding its SIV distribution throughout different landzones and different time intervals of selected lengths. The analysis scales down to hourly steps, while the longer the selected period is, the safer the results can be. The same also applies to the spatial scale, meaning that, even though the tool can zoom down to household level, the larger the area of focus is the more accurate the estimations are. This is because for greater amounts of water in WBs of wider areas or longer periods the different uncertainties, such as the uncertainty of the water meter locations and readings, play a less significant role. Additional features that can be supervised are the BAC, which for the case study of Skiathos equals to the sum of revenue water, real losses, and apparent losses.
Discerning apparent losses into unauthorized consumption and metering inaccuracies in smaller time intervals (daily) would require additional information regarding the performance of meters and transparent past reporting of unauthorized consumption incidents [15]. Such a spatio-temporal analysis can enhance the operator's understanding of the WSS and reveal the drivers of NRW and its components. The comparative diagrams of daily profiles on 1 January and 1 August for two landzones of different elevations, 0.5 m and 39.5 m, (Figure 18) can shed light on the intense seasonal and spatial variability of the investigated features: comparing the daily landzone input volumes (LIV) for the two landzones, one can notice that the landzone of higher elevation requires less LIV for both dates. The reason for this might be threefold. Firstly, the two landzones include different number of water-meters and different land-uses, secondly the elevation may have an impact both in leakage and pressure driven demand, and thirdly the infrastructure condition may also vary. Further inspection of the rest of the components can shed more light on the situation. Regarding the LIV increase for the two landzones, it can be estimated that the hilly one has an increase of 70%, while the coastal has an increase of 78%. This difference in the seasonal LIV increases of the two landzones can imply that the coastal landzone attracts much more touristic activity than the hilly one. This is actually true, since the coastal landzone in Skiathos is a part of the highly visited Skiathos promenade that hosts a large number of tourist-related small businesses, such as restaurants, cafes, pubs, and taverns. This is reinforced when revealing that the elevated landzone includes 25 water meters, while the coastal one includes eight. More evidence regarding the land uses of the town of Skiathos would require a spatial cluster analysis that would be implemented with the criterion of water-meter labeling into "households", "hotels", "café-restaurants", etc. An urban land-use clustering analysis would also improve the spatio-temporal assessment of the WB, since it would facilitate the linking of different consumption patterns with different landzones. Such a work is already conducted with the use of the self-organized maps (SOM) methodology [49]. The impacts of the intensity of touristic consumption and pressure driven consumption due to elevation difference are further reinforced, when comparing the respective BAC seasonal increases, which are 111% for the hilly landzone and 203% for the coastal one. The comparison of real losses throughout time and space may prove a very useful tool for the prioritization of the network maintenance with quantifiable evidence. For example, the analysis for the two landzones of Figure 19 offers all the evidence needed to prioritize maintenance works for the coastal landzone. Specifically, the coastal landzone has real losses as high as 70 % of SIV on 1 January that fall down to 44% on 1 August. Respectively, the hilly landzone has 56% on 1 of January, decreasing to 40%. Absolute quantities of water can also be converted in expressions relative to service connections, or main length (Table 5), which-depending on the purpose-are more useful than their expressions as total amounts. For instance, for evaluating the management of the WDN, IWA recommends the annual expression of the real losses in m 3 /connection [50]. However, such an expression would be more helpful for comparing the overall network performance to that of another network, rather than the internal comparison of two landzones, especially in the absence of DMAs. Anyway, these metrics can be indicative of what causes leakage in different landzones (e.g., altitude vs. infrastructure condition). On top of that, they can be indicative of the cause of leakage, since on one hand, leaks and bursts are more frequent on the service connection links, but on the other hand bursts in mains result in much higher flowrates [8]. Taking this into account, comparing the two landzones, one can suspect that the coastal zone might be more stressed due to bursts in mains, while the hilly landzone due to leaks in the service connections, since the coastal zone has much higher values of real losses/ service connection than the hilly, while they have almost the same real losses/ network length. The much higher value in real losses/ service connection could be justified either by more incidents of bursts or Absolute quantities of water can also be converted in expressions relative to service connections, or main length (Table 5), which-depending on the purpose-are more useful than their expressions as total amounts. For instance, for evaluating the management of the WDN, IWA recommends the annual expression of the real losses in m 3 /connection [50]. However, such an expression would be more helpful for comparing the overall network performance to that of another network, rather than the internal comparison of two landzones, especially in the absence of DMAs. Anyway, these metrics can be indicative of what causes leakage in different landzones (e.g., altitude vs. infrastructure condition).
On top of that, they can be indicative of the cause of leakage, since on one hand, leaks and bursts are more frequent on the service connection links, but on the other hand bursts in mains result in much higher flowrates [8]. Taking this into account, comparing the two landzones, one can suspect that the coastal zone might be more stressed due to bursts in mains, while the hilly landzone due to leaks in the service connections, since the coastal zone has much higher values of real losses/ service connection than the hilly, while they have almost the same real losses/ network length. The much higher value in real losses/ service connection could be justified either by more incidents of bursts or by the elevation difference.
Additionally, the per-service-connection expression is more useful for conclusions regarding the level of touristic consumption and pressure-driven consumption, since the need to specify the number of households within a landzone is eliminated. In particular, comparing 0.072 m 3 /service connection to 0.450 m 3 /service connection for the hilly and coastal landzones respectively is more safe for making the conclusion that the coastal zone has more consumption in the nontouristic winter due to pressure-driven demand, while comparing the corresponding amounts in the summer, 0.152 and 1.362 m 3 /service, can shed light on the more intense touristic activity of the coastal landzone.  Table 6 present various KPIs indicative for the performance of the simulated PM scheme. The PM scheme does not involve any DMA division, but only a temporal adjusment of the pressure to satisfy the pressure needs of the node that for each timestep is estimated to have the lowest pressure of the WDN. The pressure at the starting point of the network, at the PRV, is set to keep the pressure at this virtual critical point higher than five meters. This may not be respective to the relevant legislation, that requires pressure higher than 16 meters (=4m × (3 floors + 1), but according to the water utility, it is rather realistic and applies to the current practical requirements of the WDN operation. At this point, it should be noticed that the PM scheme is not suggested, but rather applied to make a proof of concept. Figure 20 shows the effect of the PM sheme in a high temporal scale of four trimesters. The annual potential leakage reduction is estimated to reach 51,300 m 3 , which corresponds to 21,763 kWhs of energy for pumping, transfering, treating and other operational costs, and to 5671 euros of revenue (Table 6). By the trimester analysis, it can be deduced that the highest potential for leakage reduction exists in April to June and July to September, 17.5 % and 16 %, respectively (Figure 20, right). The absolute amounts in cubic meters are presented in the left of Figure 20, where the deep blue shaded area shows the quartely variation of the leakage reduction potential. Table 6 presents the trimester values of energy savings and revenue for each trimester as well as two KPIs that often are interesting to water utilities, the night pressure decrease and the pressure variation decrease. The two KPIs that involve pressure are more relevant to the protection of the WDN against bursts.  Figure 19 proves that a spatial analysis of the critical designing parameters can be an important tool that quantifies needs and prioritizes improvements for the water utility to help in asset management. The constructed tool can quantify the benefits of a PM scheme in time and in space, at the scale of a landzone. By comparing the two indicative landzones (or any other neihboring landzone), the operator can realize that the touristic zone in the promenande of Skiathos has significantly more room for improvement in terms of leakage reduction, than the elevated zone. At the same time, comparing the temporal differences it can be deduced that the amounts of leakage that can be saved through a PM scheme can double on 1 August compared to 1 January for both landzones.
Water 2020, 12, x FOR PEER REVIEW 27 of 32 Figure 19 proves that a spatial analysis of the critical designing parameters can be an important tool that quantifies needs and prioritizes improvements for the water utility to help in asset management. The constructed tool can quantify the benefits of a PM scheme in time and in space, at the scale of a landzone. By comparing the two indicative landzones (or any other neihboring landzone), the operator can realize that the touristic zone in the promenande of Skiathos has significantly more room for improvement in terms of leakage reduction, than the elevated zone. At the same time, comparing the temporal differences it can be deduced that the amounts of leakage that can be saved through a PM scheme can double on 1 August compared to 1 January for both landzones. Figure 19. Beneficial effect of a theoretical PM scheme in two landzones of different elevation presented through different KPIs: hourly comparative curves of leakage for actual pressure scenario and PM scenario, the curve of the respective leakage reduction, daily values of leakage reduction, and daily values of savings in energy and euros, for 1 January and 1 August.
The reduction of PDD is also presented in Table 7, in macro scale and in Figure 19 for specific landzones and dates. For the whole town, a total annual PDD saving potential of 53,730 m 3 is estimated, with much higher savings-more than four times-from April to September (43,430 m 3 ), than from October to March (10,300 m 3 ). Here, it should be noted that the total annual savings of SIV due to the applied PM scheme (leakage reduction plus PDD reduction) are calculated to reach a Figure 19. Beneficial effect of a theoretical PM scheme in two landzones of different elevation presented through different KPIs: hourly comparative curves of leakage for actual pressure scenario and PM scenario, the curve of the respective leakage reduction, daily values of leakage reduction, and daily values of savings in energy and euros, for 1 January and 1 August.
The reduction of PDD is also presented in Table 7, in macro scale and in Figure 19 for specific landzones and dates. For the whole town, a total annual PDD saving potential of 53,730 m 3 is estimated, with much higher savings-more than four times-from April to September (43,430 m 3 ), than from October to March (10,300 m 3 ). Here, it should be noted that the total annual savings of SIV due to the applied PM scheme (leakage reduction plus PDD reduction) are calculated to reach a percentage of 13.3% of the initial SIV. Comparing the PDD reduction in the two landzones, it can be concluded that in the hilly residential area, the summer PDD reduction is 3.5 times the respective winter amount, while in the coastal tourist area the summer reduction is 5.4 times the winter value. This indicates the much higher potential of PDD reduction in the coastal tourist zone. In Table 7, two alternative expressions of PDD reduction are introduced. Expressing PDD in m 3 per m 3 of BAC + apparent losses is more helpful when it comes to drawing conclusions on the impact of altitude in PDD reduction. In these two landzones, it seems that the specific scheme has more impact in higher altitudes than in lower. This is reasonable because the chosen PM scheme, which is lumped regarding DMA division, causes much greater percentage pressure reduction in high attitudes than in lower ones. This fact indicates that a multi-DMA scheme would be beneficial, since with a single DMA relative savings are lower where the potential for savings is higher. The second expression introduced for PDD reduction is in m 3 per service connection. This is a more helpful expression for assessing the impact of urban land uses in PDD reduction. In particular, we can see that the touristic zone, with all the restaurants, bars and cafes, has much higher PDD reduction both in winter and summer than that of the residential zone.
Water 2020, 12, x FOR PEER REVIEW 28 of 32 altitudes than in lower. This is reasonable because the chosen PM scheme, which is lumped regarding DMA division, causes much greater percentage pressure reduction in high attitudes than in lower ones. This fact indicates that a multi-DMA scheme would be beneficial, since with a single DMA relative savings are lower where the potential for savings is higher. The second expression introduced for PDD reduction is in m 3 per service connection. This is a more helpful expression for assessing the impact of urban land uses in PDD reduction. In particular, we can see that the touristic zone, with all the restaurants, bars and cafes, has much higher PDD reduction both in winter and summer than that of the residential zone. Through the results that have been presented, it becomes apparent that a spatio-temporal analysis of the key SIV components and a series of relevant KPIs can make a conclusively insightful tool for the improvement of the WDN operation. The quantification of real losses, apparent losses, authorized consumption, the respective quantities in energy consumption and revenue mapped spatially and temporally can give a full picture regarding the weaknesses of the WDN and help the operator prioritize asset replacement or renovation needs. Through such an analysis, the localization of leakage, metering inaccuracies and theft can be realized at least at the scale of a landzone. Additionally, through such a tool the benefits of a PM scheme can be quantified in spatial and temporal detail.
The potential of this tool can extend to the investigation for an optimum DMA division scheme Through the results that have been presented, it becomes apparent that a spatio-temporal analysis of the key SIV components and a series of relevant KPIs can make a conclusively insightful tool for the improvement of the WDN operation. The quantification of real losses, apparent losses, authorized consumption, the respective quantities in energy consumption and revenue mapped spatially and temporally can give a full picture regarding the weaknesses of the WDN and help the operator prioritize asset replacement or renovation needs. Through such an analysis, the localization of leakage, metering inaccuracies and theft can be realized at least at the scale of a landzone. Additionally, through such a tool the benefits of a PM scheme can be quantified in spatial and temporal detail.
The potential of this tool can extend to the investigation for an optimum DMA division scheme where that can be achieved by a multiobjective optimization algorithm that would involve leakage minimization, energy consumption minimization, revenue maximization, and other KPIs that are presented [51][52][53]. Such an attempt would firstly require the decoupling of energy consumption and revenue from leakage. The two quantities shall alternatively be estimated also as a function of elevation, taking into consideration the dynamic energy needed to transfer water to higher levels, rather than calculate them with use of a per-cubic-meter flat rate throughout the whole WDN. This can suggest potential future work.

Conclusions
A spatio-temporal simulation of a WDN is implemented, with the use of EPANET. In order to decouple the demand-driven and pressure-driven components of SIV, leakage "valves" are used in every EPANET node. The leakage is estimated through the IWA-suggested night-flow approach that is based on the night-flow consumption. The component analysis approach is also used to validate the first estimation. A scheme of iterative processes has been constructed with two nested loops. The inner iterates until leakage converges to a given value. The outer iterates until the simulated pressure at specific points matches the actual known pressure values. The results of such an analysis can offer conclusive insight to the operator of the WDN regarding a series of designing parameters and KPIs of the network. The produced toolkit can serve as supervision support for the WDN or the basis for construction of a digital twin. Leakage, consumption, and theft can be localized at a relatively small scale, while the intense seasonal impact of tourism on the WDN can be quantified, and a series of useful conclusions can be conducted. Some of these conclusions regard the impact of elevation, the change and impact of land uses, and the respective amounts of energy consumption and revenue. The annual and trimester estimated WB tables depict the intense seasonal variability of the components and highlight the need to increase the time scale in the relevant assessments at case studies with seasonally variable demographics or weather, such as the Greek tourist islands. For Skiathos case study the analysis has shown that the annual real losses are from 50.9% to 52.2% and the NRW from 57.2% to 58.00 %. The first trimester shows the highest percentage of real losses, from 63.3% to 64.7%, while the third trimester shows the highest percentage of BAC, from 51.8% to 52.1%. The theft is estimated to be higher than 3.6%, while it increases significantly from July to December, both in absolute values and percentages. A theoretical PM scheme is applied to prove the value of such a tool, for quantifying the relevant benefits spatio-temporarily. Leakage reduction, PDD reduction, energy and economic savings are four KPIs examined for the PM. An annual potential of 51,300 m 3 in leakage reduction is estimated with the application of the PM scheme. The highest leakage reduction in absolute values can happen from April to September. Annual energy savings, economic savings, and PDD reduction potentials are estimated equal to 21,763 kWh, 5761 euro and 53,730 m 3 , respectively The form of the most appropriate expression for each KPI is also investigated, concluding that different expressions help reach conclusions in different aspects. Funding: The work described in this paper has been conducted within the project Water4Cities. This project has received funding from the European Union's Horizon 2020 Research and Innovation Staff Exchange programme under grant agreement number 734409. This paper and the content included in it do not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of its content.