Pre-proof WaterROUTE : a model for cost optimization of industrial water supply networks when using water resources with varying salinity

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

and Cresswell, 2011; Sieber and Purkey, 2015; Sonaje and Joshi, 2015. 126 In this study we present WaterROUTE (Water Route Optimization Utility Tool & Evaluation), a model that 127 adds new functionality to the previously developed WSN model (Willet et al., 2020). The original WSN 128 model generates regional water supply networks only based on water quantity requirements. In the work 129 presented here, the previously developed WSN model is extended to include water quality, specifically in 130 terms of salinity. The addition of water quality as a design criterion for water supply networks is crucial 131 to design regional decentralized water supply networks. The inclusion of water quality makes the delivery 132 of water at the desired quality possible by mixing. In this study WaterROUTE is used to demonstrate how 133 brackish/saline groundwater resources, exploited at sustainable yields, can serve as potential alternative 134 water resources for industrial use. Brackish water resources can ensure a sustainable water supply when 135 combined with optimal network layouts and desalination (Caldera and Breyer, 2017;Reddy and 136 Ghaffour, 2007). For the first time, to our knowledge, we present and apply a modeling approach to 137 create water supply network layouts with optimal pipeline routing at a high spatial resolution, connecting 138 supply sources with different salinities, which also accounts for pipeline placement costs. 139 WaterROUTE optimizes water supply network configurations according to site-specific demands for water 140 quality and quantity with water supply sources that have different and variable water qualities. With this 141 functionality we connect regional hydrological modeling with planning of water supply infrastructure. The 142 model generates the optimal network configuration and quantity of water needed from each supply 143 source to satisfy the (industrial) demand without trespassing sustainable limits for water extractions. 144 WaterROUTE is a valuable tool for IWRM and regional planning in areas where maximum sustainable 145 yields of aquifers need to be enforced. Areas of particular interest are freshwater scarce areas with 146 intensive industrial activities for which lower quality water can be used and where alternative 147 (ground)water resources are available. We present and example simulation with regional brackish 148 groundwater resources as the alternative water source for an industrial site. 149 2 Methodology 150 WaterROUTE is an optimization and visualization model which calculates optimal water supply network 151 configurations. The optimization model mixes water streams with different qualities to supply a single 152 demand location with a desired water quality. Mixing of water is a new and essential functionality to 153 design decentralized water supply networks that use alternative water resources with different qualities. 154 In WaterROUTE, water supply and demand sites are represented as vertices and pipeline connections are 155 represented as edges. The vertex and edge representation of (water) transport networks is commonly 156 used for optimization (Mala-Jetmarova et al., 2017) and was previously used for network design without 157 mixing in Willet et al. (2020). 158 WaterROUTE requires two inputs: the available water sources in a region and a preliminary network from 159 which the optimal network configuration is selected. The preliminary network is created by determining 160 the lowest cost routes between demand and supply locations using geographic information systems 161 (GIS) methods. The inputs are processed by the WaterROUTE optimization model to yield the network 162 configuration with the lowest cost for a specific water demand at the demand location (Section 2.1 -163 2.6). The outputs are then visualized with GIS software for evaluation and decision making. An overall 164 representation of WaterROUTE 1 is shown in Figure 1.   (Hirsch and Dantzig, 1968;Kim and Hooker, 2002). In this study, we alter the original FCNFP formulation 169 to include water quality parameters as a constraint. The water quality parameter included in this study is 170 the salinity (the chloride concentration) of groundwater. Water may mix throughout the network, yet the 171 water reaching the demand location must not exceed the maximum salinity defined by the user. 172 The WaterROUTE optimization problem is described as a planar mathematical system represented by 173 vertices ( ) and edges ( ). Vertices represent supply locations, demand locations, and transport 174 hubs/junctions. Edges represent the possible pipeline connections between vertices. Each vertex ( ) has 175 a chloride concentration ( ) and an associated water supply or demand ( water ( ) and flow of product ( ), the product is the amount of chloride (in mgClday -1 ). The total 184 product is determined based on the concentration and amount of water extracted from each supply 185 Figure 1 The model framework of WaterROUTE vertex. Each edge has an associated cost per unit flow per km ( ), and a length in km ( ). Additional 186 parameters are the maximum flow ( ) and maximum product capacity ( ) over each edge. We define 187 a maximum allowed concentration ( ) that is used to constrain the final quality of the supplied water. 188 Table 1 gives an overview of the parameters in the WaterROUTE optimization problem formulation. 189 The interaction between the available pipeline diameters, flow requirements, and flow velocity leads to 198 investment costs which increase with a stepwise pattern (see Supplementary Information 1). A stepwise 199 increase in costs is referred to as a stairwise arc cost function (Bornstein and Rust, 1988;Du and 200 Pardalos, 1993;Holmberg, 1994). In this study, pipeline diameters with increments of 100 mm and a 201 maximum flow velocity of 1.5 m s -1 are used. The steps in the cost function were determined for a flow 202 range between 0 and 5.5 Mm 3 year -1 by increasing the flow with steps of 0.1 Mm 3 year -1 with a peak 203 factor of 1.5 ( Supplementary Information 1). The stepwise behavior is incorporated in the objective    Table 2 (2) The amount of product (mass) extracted from a supply vertex must be equal to the amount of water 251 (volume) extracted from that vertex times the concentration (mass/volume) at that vertex if the vertex 252 is a source ( ). We ensure this with 253 254 The constraint in Eq. (7) ensures that the water extracted from a supply vertex ∑ ∑ times the 255 concentration at the vertex is equal to the product extracted ∑ ∑ .

256
The amount of product extracted from any supply vertex must be lower than or equal to the amount of 257 product extractable from that vertex 258

259
The product available at a supply vertex is determined by multiplying the concentration at the vertex by 260 the amount of water available 261 Similar to the water flows for a supply site functioning as a transport hub, the sum of the product flows 262 towards the vertex must be lower than or equal to the sum of the product flows exiting the vertex. This 263 is achieved with 264 265 where ∑ is the outgoing product flow and ∑ is the incoming product flow. 266 The product flow (mgClday -1 ) towards the demand site ( ) divided by the water volume (m 3 day -1 ) 267 towards the demand vertex must be lower or equal to the maximum allowed concentration. We ensure 268 this with 269 which is similar to Eq. (7), but the equality condition is replaced by an inequality condition and Eq. Physical bounds: constraints (4), (11) Product conservation: constraints (7), (8), (9), (10).
Special case: minimum salinity network determination 2.6 286 When the desired water quality is set to the minimum salinity achievable for a certain demand (see 287 Supplementary Information 3) the supply sources can be determined before the network configuration 288 optimization. Supply sites are ordered by increasing salinity and the cumulative water availability and 289 associated salinity are calculated. The set of clusters which can supply the demand at the minimum 290 salinity is determined from the cumulative water and salinity list. Clusters not in the set are removed 291 from the WaterROUTE optimization model inputs and the optimal network is determined by omitting the 292 water quality constraints. This procedure reduces calculation time considerably for large networks. 293 year -1 . The inputs for the example simulation are the available local groundwater sources (Section 3.1) 300 and the preliminary network layout between the demand and supply locations (Section 3.2) 301

WaterROUTE example simulation inputs
Groundwater salinity and availability 3.1 302 The groundwater in the example simulation comes from several well clusters identified based on the 303 fresh-salt groundwater interface as well as the transmissivity, which affects the possibility to extract 304 water, of the groundwater system in the region (see Willet et al., 2020). The regional groundwater 305 system has been extensively monitored, mapped and modelled in the past and shows the presence of 306 fresh groundwater resources on top of groundwater with a higher salinity (Delsman et al., 2018). PEST, the most widely used calibration software for groundwater in the world (Doherty, 2005). 332 Parameters that have been changed during the calibration are the horizontal hydraulic conductivity of the 333 aquifers, the vertical hydraulic conductivity of the aquitard, the hydraulic resistance from/to the surface 334 water system and finally the groundwater recharge. The results for Zeeuws-Vlaanderen are as follows: 335 the median of the difference between the calculated minus the measured heads changes from 0.18 m to 336 -0.009 and the average absolute difference between the calculated minus the measured heads changes 337 from 0.29 m to 0.24 m. We believe these differences are good enough calibration results. Validation of 338 the model has not been performed as the entire dataset was believed to be needed for the calibration. 339 In this study, the 3D groundwater model simulates the effect of multiple brackish groundwater 340 extractions (used as the alternative water supply source) over the well clusters on the groundwater 341 salinity over time and the piezometric heads in the vicinity of well clusters. In Willet et al. (2020), 342 analytical equations were used to estimate the upconing of the interface between fresh and saline 343 groundwater (Dagan and Bear, 1968) and the drawdown of the phreatic groundwater level (Bruggeman, 344 1999). The numerical 3D groundwater model incorporates hydro(geo)logical details of the local setting 345 (the heterogeneous salinity distribution, interaction with the surface water system, geology), includes 346 changes in groundwater salinity due to extraction wells, and thus produces more accurate results than The surface water boundary is modelled with a fixed salinity concentration and does not change over the 353 entire simulation period. There is not enough surface water salinity data to insert a seasonal varying 354 surface water salinity boundary condition (though the model can do so; a seasonal varying surface water 355 head boundary was modeled). Several surface water features in the Zeeuws-Vlaanderen region are 356 draining from the groundwater system or are not active in summer. The fresh groundwater recharge is 357 likely a dominant source of fresh water that enters the wells, given that small water courses and ditches 358 are the main surface water phenomena in this region. 359 To meet environmental targets (e.g. Natura2000) and to limit drought effects, the maximum drawdown 360 of the phreatic groundwater level is set to 50 mm (  with the lowest total costs for a specific demand at the demand site. The selected subset is the optimal 391 WSN configuration for a specific scenario. The preliminary network layout in this study represents a 392 water supply network which still needs to be built but the same methodology can be applied for an 393 existing (to be expanded) water supply network. The preliminary network for the Zeeuws-Vlaanderen 394 region was generated following the steps as outlined in Willet  The preliminary network has a total of 408 pipeline segments and 243 transport hubs to connect the 25 406 groundwater supply locations and the single demand location (see Supplementary Information 7). 407

408
WaterROUTE is used to generate the optimal water supply network configuration for five demand 409 The extraction rate from well cluster 14 is capped at 97% corresponding to a flow of 548 m 3 day -1 with a 443 pipeline diameter of 100 mm (see Table 2). Increasing the flow of this cluster would require a pipeline 444 diameter of 200 mm, leading to higher costs. The amount of water that can be extracted from well 445 cluster 10 is flexible and can be increased from 98% to 100% without increasing or decreasing the 446 network investment costs. This flexibility can be used to provide slightly more water but leads to water 447 with higher salinity at the demand location.  supply network. The general trend is that the length, complexity, and costs of the supply network 454 increase when groundwater with a lower salinity is required at the demand location. This trend is the 455 most pronounced for the 2110 simulation; the optimal network for a demand quality of 425 mgCl -L -1 is 456 17.2 km shorter than for a demand of 375 mgCl -L -1 (see Figure 8 and Table 3). The salinization of well 457 clusters leads to longer networks and increasing costs over time. 458 The networks A2, A3, and B3 (see Figure 8)  The networks A1, B2, and C3 (see Figure 8) also share the same configuration. This network configuration 464 of 51.6 km is 0.1% more expensive than the 51.8 km network (A2, A3, and B3). A shorter network can 465 have higher costs depending on the specific pipeline diameters which need to be used. The salinity of the 466 groundwater supplied by this network is lower than the required salinity of the demand site for any of the 467 periods shown in Figure 8. For example, the chloride concentration of groundwater provided by network 468 A1 is 337 mgCl -L -1 instead of the maximum allowed concentration of 375 mgCl -L -1 . This is possible due to 469 the constraint in Eq. (10) which ensures that the salinity of groundwater reaching the demand location is 470 lower than or equal to the demand salinity. A network which supplies groundwater with a lower salinity 471 than the demand salinity, the case for A1, B2, and C3, only occurs when the lowest cost network happens 472 to yield a lower salinity. 473 The networks B1, and C2 (see Figure 8)  As the maximum allowed salinity at the demand site increases, the length and costs of the WSN decrease 496 (Table 3, from top to bottom). Increasing the maximum allowed salinity increases the number of well 497 clusters which can be used in the network. A larger number of usable well clusters increases the 498 probability that well clusters located close to the demand location can be used. The possibility to choose 499 well clusters close to the demand location results in shorter networks. Shorter networks generally have 500 lower costs, with some exceptions (see Section 4.3). 501 In general, the WSN costs increase when the same water quality needs to be supplied further in the 502 future (Table 3, from left to right, except for the minimum salinity network). This is the result of the 503 salinization of the well clusters in the study area. Salinization of well clusters results in fewer clusters 504 which can contribute to the network for a specific demand salinity. As a result, longer networks which 505 transport fresh water from further away are needed. 506 The optimal network configuration for a specific region depends on the local groundwater availability and 507 the groundwater salinization/freshening dynamics. Salinization and freshening of specific well clusters 508 can lead to unexpected WSN costs. For the Zeeuws-Vlaanderen simulation this is reflected in the optimal 509 network configuration for the minimum possible salinity at the demand location when only using 510 groundwater. Based on the modeled changes in groundwater salinity in Zeeuws-Vlaanderen the minimum 511 salinity network for 2110 has lower costs compared to 2045 and 2030. 512 The modelling approach presented in this study expands the functionality of the WSN model (Willet et 518 al., 2020) with the possibility to mix water and further expands the modelling toolbox on which 519 Integrated Water Resources Management is reliant (Srdjevic et al., 2004). Determining the most cost-520 effective network for a specific quality at the demand site needs to consider different water qualities and 521 water quantities at the supply sites. Input data on water quantity and water quality is supplied by 522 existing and tested external hydrological models. WaterROUTE processes these inputs and makes it 523 possible to explore water supply network options when the water quality of regional supply sources 524 changes over time. It shows how small changes in the maximum allowed salinity of water reaching the 525 demand location cause significant changes in the configuration of the water supply network. This 526 knowledge is useful for regional planning purposes. 527 WaterROUTE can also be used to plan network expansion by using the characteristics of an existing 528 supply network as inputs. For existing networks, the capacity of the existing pipelines is fixed but using 529 these pipelines does not require new investments. Other characteristics of existing networks can also be 530 incorporated. For example, if existing networks contain segments with iron pipelines a maximum salinity 531 constraint for these pipeline sections can prevent corrosion when using saline/brackish water resources. 532 The possibility to include several demand locations, instead of a single demand location, for decentralized 533 water supply network design and regional planning is relevant for areas where multiple water users 534 compete for the same water resources. The addition of multiple demand locations, with different water 535 demand quantities and qualities, introduces non-convex quadratic constraints to the optimization model 536 and requires a problem formulation where several water flows of different qualities can flow over the 537 same trajectory in parallel pipelines. Developing an effective problem formulation for multiple demand 538 sites is suggested for future research. 539 Alternative water sources for industrial use 5.2 540 Decentralized supply networks making use of alternative water sources can be a solution to cope with 541 future changes in water availability around the world. Decentralization of water supply can enhance 542 water reuse possibilities (Leflaive, 2009) and can have advantages over centralized systems (Domènech, 543 2011;Leflaive, 2009;Piratla and Goverdhanam, 2015). Supplying industrial sites with alternative 544 regional water resources requires data on the availability of alternative sources, now and in the future. 545 WaterROUTE is a tool that can evaluate the feasibility of using these alternative sources and their 546 corresponding decentralized supply networks at a high spatial resolution. Modeled brackish groundwater 547 is used as the alternative water supply in the Zeeuws-Vlaanderen example simulation. Other 548 alternatives, such as treated wastewater, rainwater, desalinated seawater, or surface water, can also be 549 evaluated with WaterROUTE. 550 The formulation of the optimization problem in WaterROUTE is based on an overall mass balance of water 551 and a product. The product used in the Zeeuws-Vlaanderen simulation is chloride. Other water quality 552 parameters than chloride can also be used. Another possibility is to investigate multiple products 553 simultaneously by adding new variables and constraints for each of the additional products to the basic 554 model framework. This functionality is useful when evaluating other local alternative water sources such 555 as rainwater and treated wastewater from industries, urban areas, and agriculture. When adding 556 additional quality parameters non-linear and non-additive relationships between products should be 557 accounted for. For example, two water steams originally free of microbial activity, the first due to a lack 558 of nutrients, the second due to a lack of organic carbon, can lead to bacterial growth when mixed. The 559 addition of these complex interactions is only possible when they can be accurately predicted 560 mathematically but can lead to computational problems if relationships are non-linear. The fields of 561 industrial ecology (Hond, 1999)  water resources. WaterROUTE is suitable for designing water supply networks which respect sustainable 569 extraction rates. This functionality is needed for regional planning that aims to anticipate on the expected 570 changes in water availability (Hanasaki et al., 2013), salinization of (ground)water resources (H.D. the local/global carrying capacity (Bakshi et al., 2015). Within industrial water use the connection 573 between local carrying capacity and evaluation methods for water use is still lacking (Willet et al., 574 2019). WaterROUTE provides a link between the physical (hydrological) modeling of water resources and 575 regional planning of water supply networks. Through this link the costs for mismanagement of scarce 576 water resources, e.g. overextraction leading to salinization requiring longer supply networks, becomes 577 apparent. 578 In this study, the maximum groundwater extraction rates are made dependent on a maximum drawdown 579 of the phreatic groundwater level for the complete region. It is proposed to replace regional values for The results of WaterROUTE show that in most scenarios not all well clusters are used, or well clusters are 586 exploited below their maximum capacity. WaterROUTE does not yet consider the effects of partial 587 extractions on the complete groundwater system. The simulations performed for this study suggest that 588 interference between well clusters can be neglected for the Zeeuws-Vlaanderen area because well 589 clusters are far enough apart. In other areas interference may occur and simulating the effects of partial 590 extractions on groundwater salinity and drawdown to verify the feasibility of the network design is 591 suggested. Simulating the effects of a network design on groundwater and subsequently updating the 592 network design creates a dynamic interaction between the optimization model and the groundwater 593 model. Such a dynamic interaction is relevant in areas where water extractions at one well cluster can 594 affect other well clusters but is currently computationally infeasible. 595 The WaterROUTE model can also assist in designing regional water supply networks which counteract 596 saltwater intrusion from the sea by using the WSN to recharge aquifers when fresh surface water is 597 abundant. Smart groundwater extractions can lead to freshening of groundwater resources by attracting 598 fresh water from the surface water systems instead of saline groundwater from below the extraction 599 point. Coupling the operation of decentralized water supply networks with locations where this form of 600 freshening is possible can lead to regional benefits besides water supply. Fresh water resources can be 601 stored during the wet season to be retrieved at a later moment with Aquifer Storage and Recovery (ASR) 602 (Maliva et al., 2006). Correct timing of extractions can make the stored water available for use without 603 affecting the fresh-salt groundwater interface in the subsoil while being a cost effective option compared WaterROUTE is a valuable tool for planning and design of water supply networks using local alternative 626 water sources. WaterROUTE designs water supply networks that deliver water at the specified quality 627 and quantity of the user based on the modeled or known availability of water resources in a region. The 628 model is used in an example simulation to show how the dynamics of groundwater resources can be 629 connected to the regional design and planning of water supply networks. Long-term scenarios can be 630 generated which help to anticipate on changes in (fresh)water availability. WaterROUTE is demonstrated 631 with a simulation for Zeeuws-Vlaanderen, the Netherlands, and shows that a small decrease in demand 632 quality (a chloride concentration increase from 375 mgCl -L -1 to 400 mgCl -L -1 in 2110) results in a 633 decrease of the supply network placement costs by 7% for a demand of 2.5 Mm 3 year -1 . Delivering 634 higher quality water leads to higher costs because longer networks are needed. The length of the water supply network for the Zeeuws-Vlaanderen simulation varies between 46.9 km and 122.1 km based on 636 the water quality required at the demand location. The WaterROUTE model shows that costs can be up to 637 68% higher to supply water with the lowest possible salinity compared to a demand with no salinity 638 constraint in the Zeeuws-Vlaanderen simulation. The best network configuration depends on the specific 639 water quality demand of the user, the local water availability, and the time horizon over which planning 640 occurs. As water quality requirements become more stringent, optimal network selection becomes more 641 complex and modeling tools such as WaterROUTE are needed to assist decision makers in designing cost-642 effective decentralized water supply networks. WaterROUTE can, in future work, be expanded, and can 643 be used to determine the optimal balance between water transport and water treatment/desalination, 644 the use of aquifer storage and recovery within decentralized networks, and the creation of decentralized 645 water supply networks based on the exchange of water between urban, industrial, and agricultural areas. 646 Through these applications WaterROUTE can assist in coping with regional water scarcity over time by 647 connecting demand sites with local supply sources. 648

Declaration of competing interest 649
The authors confirm that there are no known conflicts of interest associated with this publication and 650 there has been no significant financial support for this work that could have influenced its outcome. Ministry of Infrastructure and Water Management and partners of the Dutch Water Nexus consortium. 663 The regional partners: DOW Benelux, Evides, KWR, and Shell contributed to this project by providing the 664 context for the example simulation, providing data on regional water demand, and practical expertise on 665 the water supply system. 666