Calibration of Design Models for Leakage Management of Water Distribution Networks

Water losses in urban water distribution networks (WDN) accelerate the deterioration of such infrastructures. The enhanced hydraulic modelling provides a phenomenological representation of WDN hydraulics, including the modelling of leakages as function of pipe average pressure and deterioration. The methodological use of such models on real WDN was demonstrated to support the planning of leakage management actions. Nonetheless, many water utilities are still in the process of designing flow/pressure monitoring, thus data available are not enough to perform detailed calibration of such models. This work presents a physically based approach for the calibration of WDN hydraulic models aimed at supporting leakage management plans since early stages. The proposed procedure leverages the key role of mass balance in enhanced hydraulic models and the technical insight on pipe deterioration mechanisms for various quantity and quality of available data. Two calibration studies of real WDNs demonstrate the feasibility of the approach and show that the distribution of leakages in the WDN does not much influence the pressure values, which confirms the need for flow measurements at monitoring districts for leakage and asset management.

1 Introduction was introduced in the 1980s (e.g. Gemanopoulos 1985), while hydraulic modelling including leakages at pipe level depending on average pressure and deterioration was introduced more than one decade ago .
Enhanced hydraulic models, born to support the management of WDN, entail a phenomenological representation of its hydraulics. This means that they capture the emerging hydraulic behaviour of the network, also in consequence of leakage management actions, e.g. pressure control or pipeline rehabilitation.
The hydraulic models calibration consists of determining parameters that, when input into a hydraulic simulation model, will yield a reasonable match between measured and predicted pressures and flows in the network.
In classic DDA models, the main parameters to calibrate were pipe hydraulic resistances, consistently with the original purpose of verifying adequate pressure, especially in firefighting scenarios. Therefore, energy balance along pipes was of primary importance to maximize the matching between simulated and measured pressure at some nodes (e.g. Walski 1983;Kapelan et al. 2007); while mass balance, i.e. matching observed and simulated pipe flows over time, was aimed at adjusting fixed nodal demands only, possibly including leakages patterns as additional fixed outflows.
In enhanced PDA hydraulic modelling pressure directly affects leakage outflows and, consequently, discharge through pipes and head losses. This means that mass-balance also depends on pressure because of pressure-dependent demand components. Therefore, beside pipe hydraulic resistance, in PDA enhanced models there are additional parameters of pressure-demand models to calibrate, including those of leakage model. Rokstad and Ugarelli (2017), reviewed existing methods for estimating leakage model parameters and reported that accurate representation of leakage in hydraulic models will often be required to distinguish between multiple groups of pipes that have different propensity to leak. This, in turns asks for detailed monitoring of flow and pressure through the WDN, i.e. at district metering area (DMA) level.
Unfortunately, in many regions worldwide water utilities are still in the process of planning flow and pressure monitoring systems. In such circumstances available data usually rely on network topology (i.e. pipe-node connectivity from existing surveys or models), asset information (e.g. roughness, diameters, age of pipes, number of connections), water consumption data based on billing accounts and total network inflow records for annual water balance purposes. In some cases, pressure records might be available at pressure/flow control devices, if any. This paper aims at answering two new technical-scientific questions which are of direct relevance at earlier stages of planning leakage management: (i) how to realistically distribute losses spatially at the pipe level in enhanced hydraulic models using limited available data; (ii) understand, for a given global leakage level, how far pressure measurements can drive the spatial distribution of water losses in the model, i.e. how different spatial distribution of losses might affect pressures in the WDN.
The innovative approach proposed herein aims at distributing water losses consistently with deterioration processes of pipelines reported in technical literature. It leverages the separation between leakages and water consumption outflows while performing mass balance in enhanced hydraulic models, and exploits the technical insight on the relationship between pipe deterioration (i.e. propensity to leak) and asset features. While doing so, it is flexible to use under various quantity and quality of available data. Discussion on two real WDNs demonstrates that pressure measurements do not provide enough information to define the spatial distribution of leakages while flow measurements are of strategic importance, giving rise to technical remarks on planning WDN monitoring systems.
2 Enhanced Hydraulic Modelling: Design model as Methodological Support for Leakage/Asset Management The design model is meant to perform the enhanced hydraulic modelling of a WDN, aimed at representing the phenomenological hydraulic behaviour of the network, which can be changed by alternative management actions. It encompasses all water demand components as function of pressure (Giustolisi and Walski 2012) including pressure-deficient conditions for water supply to consumers (Wagner et al. 1989), leakages at pipe level , private storage tanks (e.g. Giustolisi et al. 2014) or free orifices, and implements pressure control valves (PCV) (e.g. Berardi et al. 2018) or variable speed pumps controlled by any node of the network. The WDNetXL system, which is used herein on the case studies (WDNetXL-WDNetGIS 2020), integrates all such features that are mandatory to cope with the complexity of real WDN to support planning for asset management.
One main assumption of WDN design model is the steady-state simulation in each time interval ΔT, over an operating cycle (e.g. 24 h), neglecting unsteady conditions within ΔT (Todini 2011;Giustolisi and Walski 2012). At earlier stage of planning ΔT = 1 h is usually adopted as it is a technically sound compromise between the average variability of demands over 24-h operating cycle observed during one year and the uncertainties on day-by-day changes in water consumptions. The following equations reports, in matrix form, the mass balance and energy balance equations (Giustolisi 2020) behind the enhanced hydraulic model.
For each simulation snapshot t, Q p is the column vector of unknown pipe flow rates; H n is the column vector of unknown nodal heads; H 0 is the column vector of known nodal heads and V n is the column vector of volume outflows during ΔT lumped at nodes. A pn , A np and A p0 are topological incidence sub-matrices of the general topological matrix, link-node, of the network. The subscript p and n relate to the number of pipes and nodes (unknown heads), while the subscript "0" refers to the number of reservoirs (known heads). The second Eq. (1) entails mass-balance at model nodes with V n = d n •ΔT where d n are the stationary demand components (e.g. Giustolisi and Walski 2012) generally depending on pressure through H n . The following equation explicitly reports all demand components that are included in V n .
V n cons (t, H n (t)) is the water demand supplied to consumers directly connected to the water network (Wagner 1989). It equals the statistical water requests if nodal pressure is higher than the minimum value to supply service to consumers (i.e. as established by contract). V n priv-tank (t, H n (t)) represents the volume of water feeding private storage tanks in ΔT, as a common water supply scheme in many areas worldwide (e.g. Giustolisi et al. 2014).
V n orif (t, H n (t)) is the water volume from uncontrolled free orifices, e.g. from hydrants or known/assumed "bursts orifices". It has to remark that known/detected burst orifices are not included in a design model since can be reasonably assumed as repaired at the time of model development, while future bursts are unpredictable random events.
V n tank (t, H n (t)) represents the water volume to/from tanks whose consistent and stable simulation affects the global water balance in the system. The solving algorithm of model in (1), including the water level variation in tank among the unknowns can be found in .
V n leak (t, H n (t)) represents leakage volume due to pressure-dependent outflows from holes, cracks, joints and fittings. It includes both low discharge background leakages and losses from undetected/unreported pipe bursts, whose location is unknown (Farley and Trow 2003). From WDN management perspective, undetected bursts and background leakages represent the highest volume of lost water on annual basis since they run for longer time before repair; for such reason V n leak (t, H n (t)) is referred also as volumetric leakages (Berardi et al. 2018). Two main mathematical formulations for V n leak are based on conceptualization of leakages as free orifices: the Germanopoulos' (e.g. Germanopoulos 1985) and the FAVAD models (May 1994;Van Zyl and Cassa 2014). The WDN hydraulic model used herein assumes the Germanopoulos' formulation in Eq. (3) ) in terms of leakage outlet volume along the kth pipe. Since the exact location of leaks is unknown, Eq. (3) assumes that leakages depend on average pipe pressure P k,avg (t), as the mean of pressure computed at pipe ending nodes.
Each element of the vector V n leak (t, H n (t)) in Eq.
(2) equals the sum of half of V k leak in Eq.
(3) of all pipes joined at the ith node. From calibration perspective, Gemanopouos' leakage model in Eq.
(3) introduces parameters α k and β k . Exponent α k generally relates to the stiffness of pipe material. In absence of detailed monitoring and filed data, in earlier leakage management stages, α k = 1.0 can be assumed for all pipes in the network (Rokstad and Ugarelli 2017).
The coefficient β k is a deterioration parameter, related to the propensity to leak of the kth pipe due to a combination of factors including ageing, laying conditions, number of connections to private properties, diameter, material, external and internal stresses. It was recognized that differentiating β k among pipes improves the physically consistent distribution of leakage within the WDN and affects the global flow of water volumes (masses) within the network. This aspect is of key relevance while using enhanced hydraulic modelling to figure out the expected impact of various actions like, for instance, operating pressure control devices, closing gate valves (e.g. at DMA borders) or pipe rehabilitation.

Remarks on Global Leakage Assessment in WDN to Plan Investments
Technical-scientific literature reports few performance indicators of WDN asset management based on water losses (e.g. Alegre et al. 2013). The linear leakage index, which is the volume of water losses [m 3 ] per day and per km of pipeline is known to be an effective indicator of water losses. Such index uses the global estimates of water losses figures at annual scale and can be easily applicable also in those system with scarce monitoring and system knowledge.
As from previous sections, leakage outflow q k = V k leak /ΔT from a pipe with length length L k can be computed using Germanoupulos' model or FAVAD formulation: where C v is the outflow coefficient and A 0 is the area of the orifice representing the leakage in 1 m pipeline and m the head-area slope. Using the average pipe pressure P k,avg (t) is a technical assumption similar to the internal diameter and mean pipe velocities allowing to calculate pipe hydraulic resistance in place of the actual variation of the internal diameter for aged pipes. Based on Eq. (4), β k can be considered the propensity of pipe deterioration in terms of size of leak orifices per unit length. Assuming the total volume of water losses per day from the WDN is V WDN leak = ΣV k leak [m 3 ] (i.e. sum over all pipes and time steps of a daily operating cycle), it can write: As for β k , β WDN can be considered as the propensity of deterioration per unit length (km in Eq. (6)) and its value depends on the daily volume of leakages per km, M1 WDN , and on average system pressure, P WDN . For M1 WDN = (P WDN ) α , β WDN = 1.16 × 10 −8 [m 2-α /s] can be considered a reference valus of leak propensity. Therefore, given M1 WDN = 30 m 3 /km/day and α = 1, without impairing the generality of the discussion as demonstrated later in the text, if the P WDN = 15 m, β WDN = 2.32 × 10 −8 [m 2−α /s], meaning that leakages are related, on average, to system deterioration; if P WDN = 60 m, β WDN = 0.58 × 10 −8 [m 2−α /s] meaning that leakages, on average, are related to system pressure. In the earlier case, higher β WDN indicates that the rehabilitation plans play a relevant role in leakage management, in the latter case lower β WDN indicates that the pressure control could be more effective. The formulation in (6) has technical relevance for water companies because it allows comparing different aqueducts, enabling immediate understanding of the most effective leakage management options and assessment of relevant costs. It is worth noting that there is a relationship between β WDN (or β k ) and the laboratory studies about the hydraulic and mechanic phenomenology of leakages. In fact, based on the FAVAD formulation in Eq. (4) the parameter β 1,WDN is the deterioration indicator of the network and can write as: where L N is the leakage number (Van Zyl and Cassa 2014), accounting for the system material and average pressure, which is related to the exponent α as: Based on Eq. (8), assuming α = 1 as first approximation to calibrate the design model is consistent with the assumption that β 2,WDN is about one order of magnitude lower that β 1,WDN with average pressure in urban WDN usually lower than 100 m of water column (Laucelli and Meniconi 2015), meaning L N ≈ 1. WDN needing improvements (e.g. for the Italian regulation) have M1 WDN in the range 15÷140 m 3 /km/day; assuming P WDN in a wide range 5÷100 m, Eq. (6) returns β WDN is in the range 1.74 × 10 −9 ÷3.25 × 10 −7 [m 2−α /s]; it represents also the order of magnitude for pipe level β k values.
Next section faces the issue of determination a specific value of the β k for each single pipe depending on asset features. This means to spread β WDN , which can be determined iteratively using Eqs. (6) or (7) based on M1 WDN .

Introducing Pipe Deterioration Model to Calibrate Leakage Model Parameters
Two main approaches are reported in literature for calibrating parameters β k and α k of pressure-leakage model in Eq. (3) for each pipe. The first option accounts for the lack of monitoring data and assume α k = 1 and the same deterioration parameter (β k , k = 1, …, p) for all pipes in the system (e.g. Berardi et al. 2018). This means to spread uniformly the leakage propensity over the entire WDN and let pressure drive the spatial distribution of leakages in the system.
In the second option β k and α k can be differentiated among pipes using all available information on asset features and resorting to pressure and flow monitoring (Rokstad and Ugarelli 2017). In this case pipes are grouped into homogeneous cohorts based on known features, like diameters, age, material, etc., while water balance and pressure readings at DMA allow local assessment of β k and α k .
The approach proposed herein aims at providing a flexible calibration framework using various quality and quantity of information, ranging from data available at earliest planning stages up to flow/pressure measurements that will be progressively available as soon as monitoring and management actions will be implemented.
Technical literature pointed out that the main factors influencing pipe failure rate are the same affecting leakages in the WDN and data collected over time revealed that the rate of pipe breaks increases with leakage rate (e.g. Girard and Stewart 2007).
Few literature studies investigated the probability of pipe failure events, aimed at establishing functional relationships between the occurrence of reported failures and some covariates representing asset features (e.g. pipe diameter, material, number of connections, length) or surrounding conditions (e.g. burying depth, soil) (e.g. Kleiner and Rajani 2001) using statistical approaches or data-driven modelling techniques. Eq. (9) reports one of these formulas (e.g. Berardi et al. 2008) where bursts occurrence for each pipe BR k depends on total length Lt k of pipeline, length weighted age Ae k and diameter De k of homogeneous cohorts of pipes.
Such model confirmed the technical evidence that smaller pipes are likely to experience higher failure and leakages. Values of n and m, were found between 0.5 and 2; while for most analysed networks it was found exponent l = 1 consistently with the increased probability of failure for longer pipe cohorts. Accordingly, the pipe length L k was used to scale the aggregate propensity to fail at the single pipe level. The coefficient a WDN was estimated for each network to account for other local factors/variables. It is assumed herein that β k for each pipe is proportional to the failure probability based on BR k as in Eq. (9).
A avg and D avg are mean age and mean diameter weighted by pipe length, respectively. The coefficient b WDN , similarly to a WDN in Eq. (9), is supposed to include the effects of other internal and external factors on leakage propensity in the WDN. If maintenance records on failure events are available the failure model, i.e. exponents n and m, can be developed for the peculiar network; otherwise, exponents n and m can be borrowed by literature models. From WDN model calibration perspective, such approach allows assigning each pipe with a different β k value, although only coefficient b WDN need to be estimated in order to match the global WDN water loss volume, i.e. M1 WDN .
The methodology is also scalable since can be used at DMA level, as soon as flow and pressure measurements are available and water balance at DMA provides the volume of water losses (i.e. M1 DMA ) to estimate coefficient b DMA (in place of b WDN ).

Calibration of WDN Design Model in Two Real WDN
The proposed approach was adopted on two real networks located in Southern Italy. The design models were developed as part of a procedure for planning district monitoring areas (DMA) integrated with leakage control actions. The information available were: the records of the water level and water flow rates from reservoirs over the last year; the annual water consumptions at consumers; existing WDN hydraulic models, developed 10 years earlier, containing updated information on topology and pipe asset (i.e. diameters, roughness). Flow and/or pressure measurements were only available at some pressure control valves located upstream of the main distribution network. The water utility provided information about the status of gate valves and the setting patterns of pressure reduction valves, if any.
Some pressure measurements were carried on in the network although not synchronous with inflow data and consumers' consumption data. In addition, they were collected in a specific day, thus were affected by a single demand scenario and possible local disturbances. As such, they were used to validate the calibrated model, without forcing the model to exactly reproduce the observed values.
The model of the first network, named "Apulian-B", comprises 994 nodes and 1108 pipes, with total pipeline length of about 38 km; the system is fed by one reservoir. A pressure control valve (PCV) on the main pipe feeding the network modulates pressure according to the known pressure set pattern reported in Fig. 1; a partially closed valve upstream of the PCV reduces pressure in the newest southern part of the system. Figure 2 shows that altimetry of the WDN gradually decreases from 188 m a.s.l. (red) to 119 m a.s.l. (blue). The average daily volume of water supplied to consumers is about 1589 m 3 and an estimated leakage volume of 2269 m 3 , meaning a linear leakage index of about 59 m 3 /km/day with average pressure of about 18 m. The second network model, named "Apulian-MSA" (Fig. 3) includes 1111 nodes and 1307 pipes, with total pipeline length of about 43 km. Two pipes feed the distribution part of the network from one reservoir and the elevation of the WDN is more variable than the first case, ranging from about 863 to 651 m a.s.l,.
The Pressure and flow measurements in both cases were carried on to calibrate the hydraulic resistance of pipelines feeding each network from the reservoirs as they are of preeminent importance to get correct pressure values at nodes upstream the distribution part of the network, affecting volumetric leakages.
The hydraulic resistance of pipes in the distribution part of the WDN were taken from existing models of the WDN, upon verifying their consistency with technical literature values and correcting some spurious high values, likely adopted to force calibration of the existing model based on classic calibration approach. In fact, due to multiple alternative water paths in the looped distribution part of a WDN, single pipe hydraulic resistances are actually not observable and would require detailed pressure/flow monitoring which are not available. The patterns of customers' demand were assumed as priors of the calibration problem. This is a flexible assumption encompassing billing data for each consumer based on the type of contact (e.g. households, business, industrial, etc.), local water balance at DMA level or smart metering at single consumer, according to the available information.
The information about pipe age was not completely reliable in Apulian-B and not available in Apulian-MSA; the list pipes replaced in the last 5 years was available for both systems. In order to test the proposed procedure accounting for age also, pipes in Apulian-B were assigned with age between 70 (in the old city centre) and 20 years based on the distance from the centre and rough information on age of buildings; for Apulia-MSA WDN, the information on age was neglected. Values of β k were obtained by estimating the coefficient b WDN , such that: pressure and flows resulting from models were consistent with field measurements; total leakage volume returned by the model matched the value provided by the water utility; and pressure distribution within the networks was consistent with that reported by the personnel of water utilities.
The lack of pipe failure records prevented to develop statistical models like in Eq. (9). Accordingly, the calibration procedure was repeated under various combinations of exponents n of age factor (A k /A avg ) and m of diameter factor (D k /D avg ) in Eq. (10). For Apulian-B n was assumed to vary in the set {0, 0.5, 1, 1.5, 2} and m in the set {0.5, 1, 1.5, 2}; for Apulian-MSA WDN m was assumed to vary in the set {0.5, 1, 1.5, 2}, while n = 0, i.e. no information on pipe age. In addition, new pipes (i.e. replaced during the last 5 years), were assumed with a β k equal to value 1/10 of the average β k value of other pipes.

Results of Calibration and Technical Remarks
In Apulian-B, Fig. 4 shows the total outlet water volume for each time step of a 24-h cycle, identifying the water supplied to consumers (cyan) and the volume of leakages (blue), as obtained from the abovementioned leakage model. The volume of water losses changes over time depending on both pressure control pattern (PCV) (see Fig. 1) and consumers' water requests. Quite similar graphs can be obtained for other combinations of exponents (not reported here). Figure 4 also reports the distribution of linear leakage indicator values from each pipe, as colorscale dots in the middle of pipes, which are computed as the average of volumetric leakages V k leak (t) at each simulation time step, based on β k and on average pipe pressure P k,avg (t), divided by pipe length.
As the calibration procedure was repeated by assuming n = 0 and m in the set {0.5, 1, 1.5, 2}, results show that the maximum difference in nodal pressure under various combinations of n and m is about 1.3 m at peak hour (8.00 a.m) and 0.6 m at minimum consumption hour.
Performing the same analysis for n in the set {0, 0.5, 1, 1.5, 2} and m in the set {0.5, 1, 1.5, 2}, the range of variation of simulated pressures increases up to about 2.6 m, although in most nodes the maximum pressure variation across all combinations of n and m is below 1.5 m. Figure 5 shows on the left the variation of nodal pressures for n = 0, (no pipe age information -top-left figure) and for n is in the set {0, 0.5, 1, 1.5, 2} (bottom-left figure) at 8.00 a.m. Pressure variations due to different assumptions on exponents n and m are comparable with the uncertainties of pressure measurements due to the inaccuracies of meters, the unknown minor losses nearby the measurement points and the local effects of consumers' demand variation, as pulses in the networks, as well as the seasonal average increase/decrease of water consumptions.
Top-right and bottom-right plots in Fig. 5 show, as colorscale dots in the middle of pipes, the maximum variation of daily leakage volumes from pipes among different deterioration models: for n = 0 (top-right) and including pipe age information (bottom-right). It can be noted that changing the deterioration model results into a remarkable different distribution of volumetric leakages from pipes. In addition, including pipe age information in the deterioration model increases such differences.
The analyse in Apulian-MSA, confirm similar conclusions as for the first network. The maximum pressure variation under various combinations of n and m is quite below 1 m, thus technically negligible also in face of larger differences in elevation and without pressure control in Apulian-MSA. Figure 6 shows in detail the maximum variation of nodal pressure (at 8.00 a.m.) and maximum variation of daily leakage outflow from pipes among all calibrated design models.
Results demonstrate that, irrespectively on the pressure regime in the network, pressure monitoring does not provide enough information to spatially distribute leakages within the WDN. Increasing the number of pressure measurements would not result into increased information to calibrate β k . Viceversa, flow measurements at DMA boundaries are able to capture such distribution of water volumes through the WDN and drive the calibration of β k . In WDN supplied by multiple water sources (e.g. tanks/reservoir/pumps) with as many inflow meters, the proposed methodology is expected to provide more accurate model calibration since it allows identifying flow paths that are consistent with the monitored water inflow, giving rise to more realistic distribution of pressure (and leakages) within the network.
It also hints that the selection of pipes to install flow meters should be driven, besides topological rationales, by metrological criteria aimed at minimizing the disturbances (e.g. due to possible inversion of flow) or inaccuracies (e.g. due to low flow rates) of flow measurements in order to optimize the efficiency of monitoring.

Concluding Remarks
The proposed calibration methodology leverages the physically-based information provided by enhanced hydraulic modelling where mass balance equations include pressure-dependent components of water demands. In addition, it exploits the technical observation that pipe failure rate increases as leakage rate. Deterioration parameters are differentiated among pipes using a statistical failure model from literature or, if available, maintenance records. The approach is flexible and scalable to incorporate field measurements and data that will be available in the future, thus promoting the iterative refinement of the model towards the calibration of the WDN model as-built.
The application on two real networks and the experience carried on many real WDNs hint some technical remarks on pressure and flow monitoring to calibrate a design model to support early stage leakage management.
Pressure measurements are mandatory to assess hydraulic resistance of pipelines feeding the network from sources (i.e. reservoirs, tanks, pumping stations) which have direct impact on pressure and leakages in the WDN. Nonetheless, they would hardly provide enough information to calibrate the looped part of the WDN.
Pressure monitoring provides the information on average pressure status to allow the global leakage assessment to allocate investments (e.g. as per Eqs. (6) or (7)). Pressure gauges installed in some points of the WDN (e.g. at DMA boundaries) allow verifying the consistency of modelled pressure with real hydraulic status, thus avoiding systematic errors in leakage modelling. The distribution of volumetric leakages in the WDN involves water mass movements activating many water paths in the WDN that cannot be identified by pressure monitoring only. WDN subdivision into DMA and flow measurements at DMA boundaries enables the observability of leakage distribution driving the calibration of pipe deterioration parameters. Therefore increasing the number of pressure gauges in a WDN cannot surrogate the information on mass balance that is provided by flow measurements.
Acknowledgements Data on case studies were provided by IA.Ing s.r.l. Pictures from WDNetXL awas provided by IDEA-RT s.r.l. (www.idea-rt.com).
Availability of Data and Material The data of WDNs are confidential in nature and can be provided only with restrictions. Some data are available from the corresponding author by request.

Code Availability
The system WDNetXL is available by request.
Funding Open access funding provided by Università degli Studi G. D'Annunzio Chieti Pescara within the CRUI-CARE Agreement.

Conflicts of Interest/Competing Interests None.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.