Synthesis of Heat-Integrated Water Allocation Networks: A Meta-Analysis of Solution Strategies and Network Features

: Industries consume large quantities of energy and water in their processes which are often considered to be peripheral to the process operation. Energy is used to heat or cool water for process use; additionally, water is frequently used in production support or utility networks as steam or cooling water. This enunciates the interconnectedness of water and energy and illustrates the necessity of their simultaneous treatment to improve energy and resource efﬁciency in industrial processes. Since the seminal work of Savulescu and Smith in 1998 introducing a graphical approach, many authors have contributed to this ﬁeld by proposing graphically- or optimization-based methodologies. The latter encourages development of mathematical superstructures encompassing all possible interconnections. While a large body of research has focused on improving the superstructure development, solution strategies to tackle such optimization problems have also received signiﬁcant attention. The goal of the current article is to study the proposed methodologies with special focus on mathematical approaches, their key features and solution strategies. Following the convention of Je˙zowski, solution strategies are categorized into: decomposition, sequential, simultaneous, meta-heuristics and a more novel strategy of relaxation/transformation. A detailed, feature-based review of all the main contributions has also been provided in two tables. Several gaps have been highlighted as future research directions. different sections by extending the boundaries to incorporate all aspects in an industrial plant. Several authors aimed at integrating non-water thermal streams [40,42–45], cooling utilities [4], and hot utilities (steam cycle) [40] in their methodologies and found that application of holistic approaches can bring economical and environmental beneﬁts to all parties involved. The topic remains under-addressed in the literature.


Introduction
This paper addresses heat-integrated water allocation networks. Due to the similarities between water and other mass streams [1] such as hydrogen networks [2], property-based networks [3], and more generally resource conservation networks, the terminology used in this paper is based on heat-integrated mass allocation network (HIMAN). This is to emphasize the fact that most of the methodologies presented in the literature, and in this paper, can be easily applied to other resources. In HIMAN problems involving water, integration of cooling water becomes especially important [4] as it should be considered in combination with process water to satisfy industrial demands. An example of this from the pulp and paper industry was presented by Suhr et al. [5]. Figure 1 illustrates typical water pathways in industrial pulp and paper plants and the strong interconnectivity among different water users. HIMANs have been extensively studied in the literature since their emergence in 1990 [6][7][8] with more than 100 articles covering different aspects and proposing various methodologies. Due to the growing interest in this domain, several review papers have been published that directly or indirectly study different features of HIMAN problems, in particular heat exchanger network (HEN) synthesis, heat-integrated water minimization, and wastewater treatment. Bagajewicz [9] provided an overview of conceptual (i.e., insight-based) and mathematical (i.e., optimization-based) methodologies in water and wastewater minimization [10,11] with main focus on the authors' main contributions [12][13][14]. Later, Foo [15] provided a comprehensive overview of conceptual approaches in water network design (i.e., water pinch analysis) in the 21st century covering single-contaminant fixed flowrate and fixed mass load problems applied to water regeneration, treatment and total water network design. Soon after, Jeżowski [1] published an annotated exhaustive literature review on water networks, analyzing formulations, approaches, and solution strategies from 1980 to 2010. The classification of solution strategies is modified and incorporated in the current work. Use of multi-objective optimization techniques to improve optimization and controllability of the processes together with summary of methodologies related to heat, mass and work exchange networks were provided by Chen and Wang [16]. Klemeš [17] studied recent advances in water footprinting and life cycle assessment, wastewater minimization, and heat-integrated water allocation networks. More recently, Ahmetović et al. [18] carried out a comprehensive literature review specific to HIMAN and its features covering research published until 2015. The aforementioned review papers span over 20 years of research and development; therefore, some early research directions/gaps that were highlighted could/should have already been addressed. Most of these gaps (Table 1) are addressed for mass allocation networks, however they can easily be extended to HIMANs.
The remainder of this article focuses on major features of HIMAN methodologies which are reviewed with special focus on mathematical approaches proposed after 2015. Table 2 at the end of this section provides a comprehensive overview of recent publications following the same approach employed in our previous publication [4]. For a complete review of all the related papers, the reader is referred to the published review papers [1,9,[15][16][17][18]. Table 1. Highlighted gaps in the literature on water allocation networks and their specificities.

Highlighted Gaps Remarks/Literature
Fixed concentration problems As opposed to fixed mass load problems in which outlet concentration is limited (e.g., solubility) and hence mass load becomes variable. (not extensively addressed in the literature) Multi-contaminant problems Extensively addressed by mathematical methodologies with use of nonlinear programming techniques [19] Rigorous modeling Rigorous water and treatment unit models for retrofit problems in particular (not extensively addressed in the literature) Batch-wise processes Seminal work by Wang and Smith [20], and several prominent works covering water allocation network synthesis problem for batch processes [21][22][23][24][25][26][27] (not extensively addressed in the literature of HIMANs) Non-water processes In particular hydrogen networks [28]. Topics on "resource conservation network" and "property-based resource conservation networks" are dedicated to address this particularity [29]. Retrofitting Developing methodologies for plant retrofitting considering technical and geographical constraints to find feasible and practical solutions (not extensively addressed in the literature).

Multi-period operations
Considering variations of operating condition, e.g., temperature of freshwater, over multiple time horizons (to some extent, this has been addressed by literature on batch-wise operations).

Heat integration
Extensively studied under HIMAN methodologies and is the main focus of the current article.

Interplant operations
Extensively addressed by Chew et al. [36], Zhou et al. [37], Zhou and Li [38], Ibrić et al. [39], Kermani et al. [40] in HIMAN problems. Improving solution strategies Improving deterministic approaches, application of stochastic or hybrid (combined heuristic and mathematics) approaches. Jeżowski [1] highlighted the use of sequential-decomposition techniques or combination of several meta-heuristic (i.e., stochastic such as genetic algorithm (GA)) approaches as potential directions. Holistic approaches [41] Considering synergies among different sections by extending the boundaries to incorporate all aspects in an industrial plant. Several authors aimed at integrating non-water thermal streams [40,[42][43][44][45], cooling utilities [4], and hot utilities (steam cycle) [40] in their methodologies and found that application of holistic approaches can bring economical and environmental benefits to all parties involved. The topic remains under-addressed in the literature.

Approaches
There are two main approaches in HIMAN synthesis problems: conceptual and mathematical. Conceptual approaches make use of graphical techniques and expert insight. Several conceptual approaches have been proposed in the past with focus on single-contaminant problems including, but not limited to: two-dimensional grid diagram [6,46,47], heat surplus diagram [48], water energy balance diagram [49], superimposed mass and energy curve [50,51], temperature vs. concentration diagram [52,53], and enthalpy difference vs. flow chart approach [54]. Only two conceptual approaches (concentration order and temperature composite curve [55,56] and single-temperature-peak design principle [57]) have been proposed to handle multi-contaminant problems.
Mathematical approaches, conversely, are based on superstructure derivation and optimization which take into account many interconnection possibilities in the network design. The mathematical formulation is generally non-convex mixed integer non-linear programming (MINLP) and the objective function is mainly defined as minimization of total annualized cost (TAC) of the system, including both operating and investment cost. Solving a rigorous superstructure is very complex and hence requires innovative solution strategies. Mathematical approaches and their solution strategies are analyzed in more detail in Section 3.
There exists a third approach combining the synergies of conceptual and mathematical approaches. Hybrid methodologies [58][59][60][61][62][63][64][65] were first highlighted by Bagajewicz [9] as the most effective alternative to their individual applications. Such approaches allow the use of insight-based heuristics in formulating the mathematical models and hence aid mathematical approaches in representing practical and realistic alternatives in their superstructure. Moreover, expert insight can be incorporated in the methodologies to evaluate the solutions at each stage of the solution strategy, similar to the methodology proposed by Kermani et al. [4]. Conceptual approaches can also be used as techniques for initialization of large MINLP superstructures. The research direction is mainly focused on mathematical approaches. From the optimization perspective, mathematical approaches are guaranteed to provide optimal solutions (or near-optimal in non-convex formulations) to the problem, however the feasibility of such solution(s) in practice is not guaranteed. For this reason, hybrid approaches must be the main focus for future research.

Interconnectivity of Heat and Water
The proposed methodologies (being categorized as conceptual or mathematical approaches) can be categorized into three groups considering the interconnectivity of heat and water as "separate", "sequential", and "simultaneous". In "separate" methodologies, fresh water consumption is minimized in the first step, while the water network is designed without considering temperature constraints of the network. Knowing these two, in the second step, thermal streams will be extracted for heat integration [64,[66][67][68]. "Sequential" methodologies are similar to "separate" methodologies in the fact that fresh water consumption is minimized in the first stage; however, this target is incorporated in the second step, where heat integration and water network design are performed simultaneously [6,14,42,69,70]. "Simultaneous" methodologies, on the other hand, consider all the aforementioned steps simultaneously by taking into account the trade-offs between water consumption and thermal utility consumptions.

Water Network Specificities
Single vs. multiple contaminants: With regard to the constituent of water streams, the problem can be formulated as single-contaminant or multi-contaminant. It should be highlighted that the research focusing on property (e.g., toxicity, viscosity, or acidity) integration in resource conservation approaches can also be categorized under this classification. The mathematical formulations dealing with contaminations are generally nonlinear due to the existence of bilinear terms of typeṁ u C u at the inlet of mixers, whereṁ u and C u are unknown mass flowrate and contamination, respectively. Savelski and Bagajewicz [12] showed that for single-contaminant problems, the contaminant will always reach its highest limit at the outlet of a water unit operation. Therefore, the nonlinear equality constraint at the inlet of a mixer can be formulated as a linear inequality constraint with outlet contamination fixed at its maximum value. Using the necessary condition of optimality proposed by Savelski and Bagajewicz [13] and the maximum driving force [11], Yang and Grossmann [71] formulated a linear model for targeting fresh water consumption in multi-contaminant problems by relaxing the equality constraint of a mixer to an inequality constraint. The direction of relaxation was achieved by applying the KKT (Karush-Kuhn-Tucker) conditions of optimality. They stated that this formulation will result in the exact target under certain conditions and otherwise provides a tight upper bound to the problem.
Fixed-load vs. fixed-flow problems: Water minimization problem formulations can be categorized into two groups of fixed-load (FL) problems and fixed-flow (FF) problems [15,57,72]. In fixed-load problems, water is essentially a mass transfer medium with the goal of removing a fixed amount of mass load (e.g., contamination) from a process. Cleaning processes are considered as this type of problem. Since water can enter and leave a process u at any level of contamination (C in u and C out u , respectively), water flowrate (F u ) through each process varies according to Equation (1): where WUP is the set of water unit processes. This type of problem implies equal flowrates at the inlet and outlet of each water unit processes. Nonetheless, water loss or gain can be modeled as well with additional modifications. The limiting composite curve approach [10] and mass problem table (similar to problem table algorithm in heat cascade) [73] are among the well-known conceptual approaches based on fixed-load problems. For fixed-flow problems, the flowrate through each process is fixed while the water unit process is modeled as two separate units, i.e., source and sink. Conversion of the fixed-load problem into a fixed-flow problem for single-contaminant processes is completed by fixing the flow to the limiting flowrate using Equation (1) (By setting C in i = C in,max i and C out i = C out,max i ).

Heat Exchanger Network Synthesis
HEN synthesis problem is among the most researched topics in process integration dealing with developing more rigorous superstructures while providing more efficient solution strategies. An early review paper by Gundersen and Naess [74] presents more than 200 publications on this topic while Furman and Sahinidis [75] provides a comprehensive overview of major solution strategies and studies conducted in HEN synthesis until the end of the 20th century. HEN design in HIMAN synthesis problems is different from the classical HEN synthesis problems due to the possibility of stream mixing and splitting within HENs. All conceptual approaches use the classical pinch design method after having maximized the indirect heat exchanges (non-isothermal mixing (NIM)). To understand the implication of HEN synthesis in HIMAN problems using mathematical approaches, it is vital to provide a brief summary of HEN synthesis methodologies. The complete HEN superstructure is an MINLP model with nonlinear terms in both the objective function and constraints. Proposed solution strategies are directly affected by the proposed superstructure (modification of the original superstructure by relaxation, linearization, etc.) and can be categorized mainly as sequential vs. simultaneous solution strategies: • Sequential approaches: The HEN synthesis problem can be broken down into several subproblems which is then solved successively for the minimum total HEN cost. The general approach is a three-step sequential technique. The first step minimizes the utility consumption through either conceptual techniques such as pinch design method [76] or mathematical techniques by constructing mixed integer linear programming (MILP) models [77,78] . Having the utility targets, an MILP model is formulated in the second step to minimize the number of matches between hot and cold streams which is known as the heat load distribution (HLD) problem [77,79,80]. This step can further be divided into subproblems for each pinch interval which effectively minimizes the number of heat exchanger units instead of matches. In the last step, a non-linear programming (NLP) model [81] can be solved for minimum cost of heat exchanger network subject to results of the two previous steps. Floudas et al. [81] showed that every solution of the second step corresponds to a feasible HEN design in the third step. Floudas and Ciric [82] proposed a decomposition solution strategy for solving the NLP model of HEN synthesis to global optimality using generalized Benders decomposition (GBD) given the HLD matches and the utility targets and a fixed heat recovery approach temperature. Nonetheless, the solutions of the second step will only provide a feasible match with minimum number of matches and cannot guarantee a globally optimum HEN in the third step. Many techniques exist to direct the second step toward better matching results. Implementing integer cut constraints [83] to generate many solutions in the second step with minimum number of matches or using penalty (i.e., ranking) costs for each match in the objective function of the second step are among these techniques. Spaghetti design (i.e., vertical heat transfer model) [84] can also be incorporated into an MILP model (proposed by Gundersen and Grossmann [85] and extended by Gundersen et al. [86]) for targeting and ranking matches which may result in lower capital cost.
• Simultaneous approaches: The goal is to design the HEN at once. The two major contributions in this category are the works of Floudas and Ciric [82] and Yee and Grossmann [87] which are based on MINLP modeling. The former is indeed a combination of the MILP model of Papoulias and Grossmann [77] for heat load distribution and the NLP model of Floudas et al. [81] while the latter is based on a stage-wise representation approach [88]. Several assumption are incorporated in the stage-wise approach which results in a linear set of constraints, while the nonlinearity only arises in the objective function due to logarithmic mean temperature difference formulation. However, as is discussed below, this is not the case in HIMAN due to the presence of NIM.
For mathematical approaches, the main difficulty arises in modeling the heat duty of water streams which are not known a priori. This is important as the HEN formulation is generally constructed by knowing the set of hot and cold streams in advance. Superstructure-based methodologies are often suggested to solve the problem by incorporating all possible interconnections but the computational burden for a comprehensive superstructure is often cited to be an issue [89]. To address this, a subset of water streams can be integrated with HEN. The survey of literature shows that fresh water and wastewater streams are dominantly modeled as a succession of heat exchangers and splitters, and heat exchangers and mixers, respectively. The rest of water streams, i.e., inlet streams, outlet streams, and recycling streams may or may not be included in HEN synthesis superstructure. In mathematical approaches, the dominant HEN superstructures are the modified state-wise superstructure of Yee and Grossmann [87] and the state-space superstructure of Bagajewicz et al. [90]. The former is modified by including stream splitting, stream mixing, and non-isothermal mixing options which consequently makes the problem non-convex with nonlinearities arising in both objective function and constraints. The superstructure formulation of Papoulias and Grossmann [77] has also been used by many authors [4,58] to generate feasible heat exchange matches while the final HEN design is carried out using pinch design method. The HEN hyperstructure of Floudas and Ciric [82] has also been applied by Leewongtanawit and Kim [91] within a decomposition solution strategy.

Wastewater Regeneration and Treatment
Water regeneration implies removing impurities using treatment techniques which can later be reused or recycled in the system. Generally, regeneration units are categorized as fixed outlet concentration (provides linear models [1]) or fixed removal ratio approaches. This should not be confused with treatment units which remove impurities in disposed waste due to environmental regulations; though nevertheless, the classification remains the same. Methodologies for optimization of wastewater regeneration and treatment networks can be broadly categorized into conceptual and mathematical approaches. Conceptual approaches for wastewater treatment are limited to non-heat-integrated networks; nevertheless, a short summary of these methodologies is included here as they bring insights into optimal integration of treatment units with processes. Several conceptual techniques have been developed: • Fixed-load problems: Early work on regeneration targeting in this category is based on limiting composite curve approaches [10,11,92]. However, as stated by Foo [15], these techniques could not handle all different cases that could arise. In particular, there are cases where implementing the regeneration process changed the pinch point [93] and hence cannot correctly define the minimum fresh water target. Later, several studies proposed using graphical and sequential approaches to overcome this issue, known as revised targeting techniques [94,95]. They showed that the inlet concentration of a regeneration unit is not always the same as the pinch concentration (assumption that was made in previous work). In each case, a fixed outlet concentration for regeneration units were considered.
• Fixed-flowrate problems: Hallale [96] presented a guideline for placement of regeneration units in fixed-flowrate problems. Analogous to the placement of heat pumps in thermal processes, they indicated that, to reduce the fresh water intake (analogous to reducing hot utility in conventional pinch analysis), a regeneration unit should be placed across the pinch by regenerating water with higher concentration from above the pinch (having excess water, analogous to excess heat below the pinch in conventional pinch analysis) to the lower concentration region below the pinch (water deficit, analogous to heat deficit above the pinch in conventional pinch analysis). The main conceptual methods include ultimate flow targeting, source composite curve, and automated targeting techniques. Nonetheless, separate analysis of wastewater treatment networks and water networks forbids any potential reduction in fresh water consumption.
Mathematical models of wastewater treatment and regeneration are nonlinear in nature due to the existence of bilinear terms and hence resulting in non-convex NLP formulations. Quesada and Grossmann [97] highlighted two formulation approaches in modeling general multi-component problems: considering mass flow and composition components unknown which makes the mass balance constraints of mixers nonlinear, while another approach is to model the individual flows of components which makes the mass balance constraints of splitters nonlinear. The former approach is the dominant one in water allocation networks. Quesada and Grossmann [97] proposed a linearization formulation using McCormick relaxation [98] within a branch and bound procedure to solve the problem to global optimality. Karuppiah and Grossmann [19] were the first to address the advantage of optimization of integrated water networks including wastewater treatment. Similar to Quesada and Grossmann [97], they incorporated the McCormick formulation for convex relaxation of bilinear terms in the original NLP model. They later extended their superstructure to address uncertainty within contamination generation and treatment removal ratio by proposing non-convex MINLP models and solve the problem to optimality. The detailed description of mathematical techniques in solving non-convex NLP and MINLP models are beyond the scope of this review and the readers are encourage to refer to [19,97,[99][100][101][102][103].
Dong et al. [104] were the first to address total heat-integrated water allocation networks incorporating wastewater regeneration units and their interconnections within an MINLP superstructure. They showed that the combination of deterministic and stochastic search techniques typically reached the global optimum. Yang and Grossmann [71] have proposed a linear programming (LP) targeting model for water networks involving treatment units using a similar approach to the HEN stage-wise superstructure of Yee and Grossmann [87]. Their proposed model includes as many stages as the number of treatment technologies and one pathway for each unit operation, implying no mixer at the inlet of a treatment unit. The resulting targeting superstructure does not necessarily provide an upper bound, yet does provide approximations of the optimal value of the objective function for the original NLP model. Sharma and Rangaiah [67] applied the same formulation of Bogataj and Bagajewicz [105] for regeneration units while using multi-objective optimization (MOO) through a GA minimizing the total fresh water intake and total regenerated water. Besides two studies [69,106] that have modeled the treatment units using the fixed outlet contamination approach, others have used the fixed removal ratio approach.

Superstructure Generation and Solution Strategies
In general, the synthesis problem of HIMANs is formulated as an MINLP problem. This is due to the presence of binary variables (existence of heat or mass exchange matches) and continuous variables (operating conditions, e.g., temperature and contamination levels) which are complex to solve. This necessitates the development of robust and efficient solution strategies. Several solution strategies can be applied to HIMAN superstructures depending on the interconnectivity, complexity and completeness of the superstructure. As stated by Jeżowski [1], they can be categorized into linearization, initialization, sequential, decomposition, meta-heuristics and simultaneous techniques. Simultaneous solution strategies for solving MINLP problems may exhibit decomposition or sequential techniques intrinsic to the solver being used. They are, however, categorized under simultaneous strategies. It should also be highlighted that the categorization can overlap to some degree, i.e., several techniques can be combined in a solution strategy. Table 3 provides a feature-based representation of all mathematical approaches in HIMAN synthesis problems addressing their objective function(s), mathematical formulation and solution strategy.  [37] • • MINLP • Zhou et al. [118] • • MINLP • Tan et al. [120] • • MINLP • Sahu and Bandyopadhyay [61] • • LP (min freshwater) → LP (min thermal utility) • • Renard et al. [43] • • MINLP → pinch design method •  [136,137] • • Generalized disjunctive programming (GDP) MINLP • Ghazouani et al. [138] • • MILP • Yan et al. [139] • • NLP (relaxing integers with fractional continuous variables) • Torkfar and Avami [140] • • MINLP (including pressure drops in water network) N/A

Decomposition
Two terminologies of "decomposition" and "sequential" must be clarified here. In both cases, the problem is divided into two or more steps; however, the former consists of a finite number of iterations between the steps given defined termination criteria, while the latter is a uni-directional solution strategy with no iteration. More importantly, in decomposition solution strategies (such as GBD) the results of one step provide inputs for the subsequent steps while in sequential solution strategies, no such interactions exist. Papalexandri and Pistikopoulos [7] proposed an MINLP superstructure for heat and mass exchange networks and solved it using GBD [148] by decomposing the problem into a master MILP model optimizing the network configuration and a primal NLP model optimizing the operating conditions. The MILP model provided a non-decreasing lower bound to the objective function while the NLP model gave a non-increasing upper bound. The stopping criterion was the convergence of the objective functions in the two steps below a predefined threshold. They have allocated the binary and continuous variables to master and primal problems, respectively. Similarly, Leewongtanawit and Kim [91] decomposed their MINLP model into MILP and NLP models and solved them iteratively until no further improvement (beyond a threshold) was observed in the objective function of the NLP model. They decomposed the variables similar to the approach proposed by Floudas and Ciric [82] in HEN synthesis.

Sequential
Several sequential solution strategies have been proposed which commonly optimize water and heat targets within the first or second step using LP [14,61], MILP [4,44,58], or NLP [59,60] models. Having these targets (with or without the design of water network), an MILP model can be formulated to minimize the number of heat exchange matches [77,79]. In all cases, the HEN is designed using the pinch design method. Liao et al. [111] proposed a two-step approach for targeting and design by formulating two MINLP models. The first model obtained the water and thermal utility targets with number of stream splits, while the second model minimized the number of heat exchange matches. Dong et al. [104] proposed an iterative sequential solution strategy by first solving an MINLP model using random initial guesses and later improving the results (solving the MINLP model at each step) by iterative heuristic perturbations in both continuous and binary variables. Liao et al. [115] solved an MILP model minimizing the operating cost together with number of matches (similar to HLD) and used the targets and matches as initialization for the MINLP model of HIMAN. Most recently, Ibrić et al. [142] proposed an iterative sequential solution strategy consisting of two steps. An NLP water network model was solved for minimum heat recovery approach temperature (HRAT) to provide a lower bound on the problem. In the second step, an upper bound is assigned to HRAT using a finite number of iterations and a sequence of NLP-MINLP models were solved for each value with the final solution selected as the best among all solutions. Jagannath and Almansoori [141] proposed a sequential solution strategy by introducing three MINLP models: Model A (water network with NIM), Model C (HEN synthesis [87]) and Model B (combined Models A and C). The problem was solved sequentially by solving Model A and a simplified version of Model A to find the water and energy targets. Depending on the results of the two versions, a relaxed version of Model B was solved. At this point, the flow and concentration variables in Model A were fixed and the problem was solved to generate set of solutions using techniques similar to integer cuts. Model C was applied for each solution and the final optimal solution was thus selected as the best of all solutions. The authors mentioned that their sequential approach is computationally exhaustive, yet the solutions are similar to those obtained by simultaneous approaches. More recently, Hong et al. [146] extended their targeting approach [145] by addressing multi-contaminant as well was treatment problems using a sequential solution strategy. An NLP model was formulated to minimize the fresh water consumption in a first step, providing initial values on flow rates and concentrations. An MINLP model was then solved with relaxed/linear TAC in the second step followed by optimization of the original MINLP in the third step.

Simultaneous with or Without Initialization
For simultaneous approaches, models are mainly formulated as MINLP which are solved using commercial software such as DICOPT [149], BARON [150,151], SBB or LINDO [152]. It should be noted that the algorithms underlying these solvers (e.g., outer approximation [153] in DICOPT) are based on decomposition techniques. Nevertheless the original mathematical formulation was not decomposed before solving [110,113,120,143]). Generally, an MINLP model requires good initialization which can be achieved by solving a relaxed instance of the original model either through fixing continuous variables (resulting in an MILP initialization model), fixing binary variables (mainly by excluding the HEN matches, and hence resulting in an NLP model [124,126,[129][130][131]134]), solving an MINLP model of water network with no heat integration [133], or by generating random values [38,135]. Boondarik Leewongwanawit [108] used an MILP model for initialization of MINLP superstructure by fixing the contamination loading at the outlet of water units and linearizing the HEN cost formulation. Dong et al. [104] used an NLP formulation of the water network for targeting the utility consumption and further labeling thermal streams (i.e., hot or cold) which was then solved simultaneously.

Meta-Heuristics
Du et al. [42] have employed GA combined with SA for optimizing water network (MINLP superstructure) and HEN (MINLP superstructure) in an iterative manner. They, however, neglected temperature effects in the water superstructure (first step) and hence construct (i.e., extract) the required thermal streams based on the optimized water network for the second step. Liu et al. [63] proposed a hybrid methodology combining GA-SA with mass pinch and pseudo-T-H-diagram [154]. Other meta-heuristics approaches such as particle swarm optimization were incorporated in the work of Li [123].

Relaxation/Transformation
Under this classification, the original MINLP model is transformed/simplified/relaxed by redefining/removing/adding extra constraints. Zhou et al. [135] formulated a HIMAN superstructure using MPEC and applied complementarity formulations [155] to model binary variables. They solved the problem by transforming the model into NLP. GDP has been used by Liu et al. [137] to formulate discrete and continuous variables of HIMAN by incorporating logical propositions, disjunctions and algebraic constraints. They solved the problem using BARON by transforming the GDP model into MINLP. However, their formulation does not address HEN design. Several authors [139,147] proposed HIMAN superstructures using NLP techniques. Yan et al. [139] adapted the MINLP superstructure of Ahmetović and Kravanja [126] and avoided the use of binary variables by introducing continuous variables in the form of y = f /( f + ζ) where y indicates the existence of the unit of size f and ζ is a very small parameter (∼ 10 −5 f ). Nevertheless, the superstructure a priori treated fresh water and wastewater streams as the only cold and hot streams, respectively, involved in heat integration. More recently, Hong et al. [145] proposed a targeting approach using an MILP formulation in which the HEN was designed in a single step. They adapted the HEN transshipment model [77,87] addressing stream splitting and NIM within the HEN superstructure. One other possibility is by using a discretization approach in which known variables are discretized into a set of known values which will result in an MILP model [62,156]. Several authors have highlighted the challenges in solving medium and large MINLP problems and proposed a reduction strategy and a reduced superstructure to solve heat-integrated water allocation networks more easily [65,142]. Wang et al. [65] proposed several heuristics related to contamination monotonicity (only applicable to single-contaminant problems) and pinch principles together with rational NIM in order to simplify the superstructure of Liu et al. [136]. Ibrić et al. [142] also applied several rules to eliminate infeasible and impractical variables and connections from the superstructure. They showed that these simplifications can reduce the computational efforts and hence increase the solving efficiency.
In summary, it was observed that several major superstructures have been used/proposed by previous authors for synthesizing HIMAN problems, addressed specific regions of the solution space and required specific solution strategies. One approach was the HEN hyperstructure formulation of Floudas and Ciric [82] combined with the water allocation network problem such as in Leewongtanawit and Kim [91], which was solved via a decomposition strategy. Alternatively, Papalexandri and Pistikopoulos [7] proposed a hyperstructure for heat and mass allocation networks by integrating the HEN hyperstructure [82] and an analogous mass exchange hyperstructure [157] to construct the HIMAN hyperstructure. This approach encompassed all possible interactions between the two, where each stream could be split and directed to all heat and mass exchangers while bypass streams were also included. This comprehensive model was formulated as an MINLP problem and solved using GBD. The stage-wise HEN superstructure of Yee and Grossmann [87] was the major formulation used in literature. As discussed in Section 2.4, the major assumption of the original formulation is the isothermal mixing of streams which consequently forbids many promising alternatives in the HIMAN synthesis problem. This was modified by many authors to address non-isothermal mixing and splitting. This superstructure and the one proposed by Floudas and Ciric [82] ("simultaneous match-network optimization") are categorized under simultaneous approaches. Their usage in HIMAN synthesis problems requires additional assumptions and simplifications (e.g., not considering all water stream participating in heat exchange) in the hope of alleviating the computational burden of solving the superstructure of all possible opportunities. As an alternative to superstructure representation, the state-space representation of Bagajewicz et al. [90] was addressed in synthesizing HIMANs by several authors [37,38,104,112,118]. The state-space representation contains the superstructure representation of HIMAN as a special case. Nonetheless, as stated by Bagajewicz et al. [90], this representation can alleviate some of the difficulties arising in superstructure optimization. The representation is based on the definition of a set of input and output variables (e.g., input and output temperatures of thermal streams) and their relations. A state variable, such as heat exchanger inlet temperature or mass exchanger inlet concentration, is defined as a variable which allows calculation of the output variables by relation with the inputs. The state-space is referred to as the set over which the state variables take their values. To this end, the overall input-output relations can be solved via two operators: one dealing with mixing and splitting, and the other with mass or heat exchange. It was shown that use of a particular operator-the assignment operator-can lead to NLP formulations which are better suited for solving large-scale problems.

Superstructure Extension
Despite the importance of holistic approaches to capture the trade-offs between heat and water, a limited amount of research has included non-water thermal streams in their methodologies. This potential was first addressed by Renard et al. [43], although no case study was presented. Later, Kermani et al. [44] extended their superstructure and presented a simplified Kraft mill process by incorporating non-water thermal streams and showed the large potential in their combination. This has been later applied to a real Kraft mill [4,40], also addressing interplant operations. Zhou et al. [37,118] developed a multi-scale, stage-wise superstructure addressing interplant operations for fixed-load as well as fixed-flowrate problems using MINLP. Most recently, Ibrić et al. [39] extended their superstructure [142] to address interplant operations by use of additional binary parameters to allow or forbid interplant connections. Heat transfer coefficient calculations have also been considered in the work of Torkfar and Avami [140] as a function of stream velocity. They further included pressure drop calculations in HEN and water network design. Their superstructure was formulated using MINLP which is stated as a modified and improved version of the superstructure by Jiménez-Gutiérrez et al. [128]. Use of live steam was first investigated by Savulescu and Smith [6]. They argued that use of live steam in water networks is doubly beneficial since it can reduce the steam consumption and also eliminate the use of heaters, hence reducing the capital cost. To the best of the authors' knowledge, the methodology presented by De-León Almaraz et al. [64] is the only work considering the use of vapor-state water in the water network. For implementation reasons, they modeled the phase change by only considering the sensible heat in the HEN design, while the latent heat is added to the solution via addition of a corresponding hot utility.

Physical Improvements
Importance of NIM and its effects on network performance were extensively highlighted throughout the literature [58,107,119,158]. Recently, Martinez-Patiño et al. [159] analyzed direct and indirect heat exchanges by incorporating an exergy component. They concluded that networks having the same water and energy targets may exhibit different exergy losses (due to NIM) which negatively impact the cost of heat exchanger area.
Features of non-water thermal streams and use of live steam in water networks must be studied further as they represent a more realistic approach for application of HIMAN methodologies in industrial applications.

Water-Energy Nexus
Heat-integrated water allocation network synthesis problems have been treated as a separate research field over the past decades; however, they should be regarded as a special case in the field research related to the water-energy nexus which, in turn, is part of a broader water-energy-food nexus. The water-energy nexus was first mentioned in an annual review paper by Gleick [160], highlighting the interconnectedness of energy and water. The reasoning followed that water is used in extracting and producing fuels and producing electricity via steam while energy is used to produce, transport, and purify water. In addition to the research discussed in this paper, work in water-energy nexus domain encompasses developments in the field of desalination technologies, membrane systems, water use in biorefineries, and in shale gas production. The water-energy-food nexus considers the intertwined nature of water, energy, and food by highlighting that water is used to produce food and varieties of crops which, in turn, can be used to produce biofuels. Garcia and You [161] highlighted several research opportunities related to the water-energy nexus: energy and water use in households, novel water sources such as rainwater or water being produced from extraction of fossil fuels, hydropower plants, climate studies, policy planning, and holistic approaches in design of industrial wastewater treatment networks. For comprehensive review of recent contributions and future directions related to the water-energy nexus, the reader is encouraged to refer to review papers by Garcia and You [161], Martinez-Hernandez and Samsatli [162], Lee et al. [163], Albrecht et al. [164], and Dai et al. [165].

Benchmarking Analysis
Similar to the previous review on HIMANs [18], a benchmarking analysis is carried out to illustrate the main features of different methodologies. A well-known single-contaminant case study (Table 4) originally proposed by Savulescu and Smith [6] was selected. Over 35 articles have evaluated their proposed methodologies using this case study. The results are provided in Table 5 while selected key features are plotted in Figure 2 using parallel coordinates for 33 out of 35 articles (two of them lack data to be visualized).
Several network indicators have been selected for the analysis: Number of thermal streams including thermal utilities (N th s ), number of heat exchangers (N HE ), total area of heat exchangers (A total HEN ), total number of mixing points (N mixer ), number of non-isothermal mixing points (N N I M mixer ), and number of mass streams in the water network excluding the thermal ones (N m s ) are among them. Figure 2 shows that the number of non-isothermal mixing points as well as number of mass streams (N m s ) are inversely proportional to the HEN cost, i.e., increasing either of the two will decrease the HEN cost.
This case study further illustrates the necessity of considering more indicators while optimizing heat-integrated water allocation networks. Among the 33 visualized results in Figure 2, 29 results possess 3-5 heat exchangers (the maximum is 10). This may indicate good compromise in HEN investment costs; however, for the same solutions, the number of mixers, non-isothermal mixers and total heat exchange area vary 3-13, 0-10, and 3500-6300 m 2 , respectively. These "neglected" indicators should be somehow addressed in HIMAN synthesis problems which necessitate the application of multi-criteria decision making approaches. In addition, no single solution can possess the optimal value for all indicators, which correspondingly requires the generation of a set of potential solutions that can be analyzed, instead of a single "optimal" solution.  )   2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34 (Table 5) using parallel coordinates [166]. Table 5. Benchmarking of several HIMAN methodologies using the four-water-process test case of Savulescu and Smith [6] (1) .

ID
Network Indicators (3) Economic Indicators (4) M/C ( (1) All methodologies reached the fresh water target of 90 kg/s. All but two of them reached the thermal utility targets of 3780 kW of hot utility. Savulescu et al. [47] and Martínez-Patiño et al. [52] reported 485 and 406 kW of cold utility and 4265 and 4186 kW of hot utility, respectively. Where information was not enough to calculate the indicators, "NA" is indicated. (2) "M" indicates mathematical approach, while "C" denotes conceptual approach.

Concluding Remarks and Future Directions
This work presented a meta-analysis of literature on heat-integrated water allocation networks. Key features of the proposed methodologies have been analyzed with special focus on mathematical programming approaches including HEN synthesis. Developing more rigorous mathematical superstructures necessitates proposing novel solution strategies. The proposed solution strategies have been categorized into decomposition, sequential, simultaneous (with or without initialization), meta-heuristics, and relaxation/transformation strategies. A benchmarking analysis was presented comparing the results of different proposed methodologies from the water and heat exchanger networks perspective. It illustrated how methodologies can produce different results considering other indicators than the typical TAC. Following this review, several gaps have been identified (Table 6 summarizes the main gaps): As mentioned previously, despite the importance of addressing synergies among various elements in a typical industrial plant, holistic approaches have rarely been addressed in HIMAN synthesis problems. Apart from a limited number of specific publications [4,40,[43][44][45], non-water thermal streams have not been combined in HIMAN synthesis. Future research directions should therefore focus on this aspect by proposing more rigorous and efficient superstructures. In addition, use of live steam should be investigated using improved formulations. As water is subject to heating and cooling duties, water loops have a role in recovering heat within and between processes. This feature is even more sensible when considering inter-plant operations. Moreover, following the observed gap in holistic approaches, HIMAN synthesis problems should be considered in conjunction with other heat recovery technologies including organic Rankine cycles (ORC)s and heat pumps. The literature lacks multi-period operations of HIMANs. This is an important feature considering daily and seasonal variations of operating conditions of an industrial plant, including the temperature of freshwater. Thermal storage tanks must be combined within HIMAN problems to provide a flexible heat transfer medium over time. Uncertainty analysis of HIMANs must be addressed to find resilient networks given the uncertainties in the system including costs and operating conditions. Following the benchmarking analysis, multi-criteria decision making approaches must be incorporated in HIMAN synthesis problems to find sets of promising optimal or near-optimal solutions considering diverse economic, environmental, and practical indicators. The application of stochastic optimization and hybrid approaches should be favored in this direction. Upon the survey of the literature, only one article mentioned large-scale industrial applications [4], yet the methodology is limited to the targeting step. Most of the proposed mathematical methodologies are highly complex and their applications to industrial cases may face computational challenges. Hence, research toward efficient solution strategies must be the future trend thus shifting the focus toward reaching practical and good solutions, not necessarily the global optimum. Following the highlighted gaps in Table 1, batch processes and retrofitting remain largely untreated which necessitates further research. Table 6. Summary of identified gap in HIMAN synthesis problem.

Unaddressed literature gaps
Fixed concentration problems Problems with variable mass load Rigorous modeling Water and waste treatment models Multi-period operation Considering the dynamic nature of systems Retrofitting Methods covering partial system retrofit and redesign instead of design

Newly identified gaps
Better treatment of thermal streams Considering non-water thermal streams and potential for live steam as part of the problem definition Utility integration Considering HIMAN with utility selection and integration concepts Sensitivity analysis Generation of multiple or resilient solutions in lieu of global optima Multi-criteria optimization Methods which address multiple criteria for decision-making which extend beyond minimization of cost or fresh water consumption, considering an expanded system. Large-scale problems Developing approaches to adapt formulations to larger scale problems or reformulation to encourage solution generation for problems on the industrial scale Author Contributions: This study was done as part of M.K.'s doctoral studies supervised by F.M. I.D.K. has contributed to this paper through reviewing, editing, and presentation and discussion of results.
Acknowledgments: This research project is financially supported by the Swiss Innovation Agency Innosuisse and is part of the Swiss Competence Center for Energy Research SCCER EIP.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: