## Abstract

This article aims to present a new methodology for waste load allocation (WLA) in a riverine system considering the uncertainty and achieve the lowest amount of inequity index, cost, and fuzzy risk of standard violation. To find a surface of undominated solutions, a new modified PAWN method, initially designed for sensitivity analysis, was developed and coupled with a simulation-optimization process using multi-objective particle swarm optimization (MOPSO) algorithm, to consider the uncertainty of all affecting variables and parameters by using their probability distribution. The proposed methodology applied to Sefidrood River in the northern part of Iran. Graph model for conflict resolution (GMCR) as a subset of game theory was implemented to attain a compromise on WLA among the stakeholders of a river system's quality in Iran: Department of Environment, Municipal Waste Water, and Private Sector. Some undominated solutions were used in GMCR model and modeling the conflict among decision makers reveals that their preferences and the status quo do not lead to a solely stable equilibrium; thus the intervention of a ruler as arbitrator leads them to reach a compromise on a scenario that has a median FRVS and cost. Sensitivity analysis was done using the PAWN method to assess the sensitivity of three intended objectives to all variables and parameters.

## HIGHLIGHTS

• WLA in a river with three decision makers done by considering minimizing inequity, cost, and fuzzy risk of standard violation.

• The optimum solution found by using a stochastic simulation-optimization model.

• The sensitivity of functions’ outputs to inputs has been determined by using PAWN method.

• Game theory and GMCR+ model used to discover the option agreed upon by the parties.

• Cost and standard violation have the most sensitivity to the upstream flow and inequity has the most to the dischargers' BOD.

### Graphical Abstract

Graphical Abstract
Graphical Abstract

## INTRODUCTION

Increasing human demand for water has made this vital element noteworthy, and its ownership rights to use have become more challenging. Using water bodies, especially rivers, for industrial purposes such as cooling the industries and assimilating the wastewater, is one of the consequences of the industrialization era. Dwindling the ecological capacity of the rivers to treat wastewater, population growth – which leads to producing more wastewater – and industries bloom located along the rivers have reduced the quality of the river's water. This reduction is due to the drainage of the wastewater from the riverside dischargers into it. The requirements to maintain water quality in compliance with the approved standards and preserve its capability to receive the maximum flow of wastewater from industries has led to allocating the shares of wastewater discharge to the dischargers as a controversial and problematic issue. In other words, waste load allocation (WLA) in rivers refers to the determination of required pollutant treatment levels to a group of pollution dischargers to ensure that water quality standards are maintained along the river (Ghosh & Mujumdar 2010). Determining the shares for each of the dischargers, whether industries or municipal wastewater treatment facilities or non-point sources at the river border to discharge the wastewater, can be accomplished using various methods. Optimization (Yandamuri et al. 2006; Cho & Lee 2014), water quality trading (Breetz et al. 2004; Sarang et al. 2008), and game theory (Niksokhan et al. 2009; Nikoo et al. 2012) are among the most common ways to do this.

Application of game theory in water quality management issues is preferred to conventional decision analysis methods because of two reasons. Firstly, unlike multi-criteria decision making (MCDM), this approach addresses decision-maker behavior and finds solutions that will not be refused by them. Finally, this method, unlike MCDM, is capable of expressing and reflecting the various characteristics of water resource problems, such as engineering, socio-economical, and political issues (Mirchi et al. 2010).

Optimal waste load allocation implies that the treatment levels not only maintain the water quality at standard level but also result in the best value for the objective function defined for the water quality management problem in the riverine system (Burn & Lence 1992). Chadderton et al. (1981) introduced and argued various basic theoretical approaches to apportion waste load and brought forward ‘percentage of equal treatments’ as the best option among those 20 methods. Different scholars with the different point of views tried to allocate waste load among dischargers by using various methods and tools: Hathhorn & Tung (1989), Burn & Lence (1992), Kerachian & Karamouz (2005), Saadatpour & Afshar (2007), Qin et al. (2009), Li et al. (2014), Allam et al. (2016), Mohan & Kumar (2016), Xu et al. (2017), Rafiee et al. (2017), and Afshar et al. (2018) are among the extensive range of them.

Burn & Yulianti (2001) presented a methodology that contains three models, two of which are pairwise models (cost versus standard violation, the cost versus equity) to allocate waste in a riverine system. Feizi Ashtiani et al. (2015) studied an approach to meet environmental, economic, and equity objectives simultaneously. The waste loads were allocated according to the discharge permit market. They used two discrete graphs of cost versus violation and cost versus inequity to select four scenarios based on these two models' Pareto Fronts. Saberi & Niksokhan (2017) proposed a methodology using the decision support system (DSS) for the graph model for conflict resolution II (GMCR II), multi-criteria decision making (MCDM) analysis and the multi-objective particle swarm optimization (MOPSO) algorithm while considering cost and standard violation to allocate waste load among dischargers.

The lack of considering the equity as a third objective in the optimization model leads to the setup of a three objectives simulation-optimization model. The difference between the current work and those mentioned above is that ongoing work leads to a surface of the undominated solution, where all points on that surface are optimum in terms of intended objectives. In contrast, the others are two-axis diagrams where the feasible space of the answers is a curve. This study, instead of using a deterministic model that they used, applied a stochastic model, which uses the distribution function of all variables and parameters that take part in the simulation-optimization process.

In this article, a novel methodology proposed to allocate waste load among dischargers located adjacent to a riverine system equitably and cost-effectively, as well as encompassing the lowest level of fuzzy risk of violating the standard (FRVS) for biochemical oxygen demand (BOD). Achieving a precise allocation for removal percent so that three minimizing objectives attained with a porous ratio, taking into account the uncertainty of variables and parameters affecting the simulation model, is required. Three different criteria (objectives) to allocate the waste load to three dischargers, which have controversial preferences, lead to the use of the game theory instead of conventional optimization approaches.

A stochastic simulation-optimization process was run applying MOPSO optimization algorithms, and quality simulation derived from fate and transport equations to find undominated solutions. In the last section, the GMCR+ model was utilized to find a stable equilibrium for the conflict about selecting a specific scenario by three decision makers: Department of Environment (DOE), Municipal Waste Water (MWW), and Private Sector (PS). The impossibility to attain one specific equilibrium through conflict resolution caused the entrance of the judiciary; this intervention leads to achieving one particular equilibrium from six selected scenarios from the simulation-optimization process. A sensitivity analysis (SA) by the PAWN method (Pianosi & Wagener 2015) has been carried out to assess the outputs' sensitivity to the inputs parameters and variables uncertainty in the simulation process.

Using an integrated simulation-optimization which considers all goals simultaneously to find the Pareto Front was not done before since the previous studies usually use two separated simulation-optimization process considering two goals (Burn & Yulianti 2001; Feizi Ashtiani et al. 2015). In this article, an integrated model was set up that considers the uncertainty of all parameters.

How to sum up with the obtained output series from stochastic simulation-optimization reliably and flexibly, whether the median, average, or any other statistics should be used, was one of the difficulties which led authors to use a risk-based index: FRVS. The way sensitivity analysis is used for three intended functions output was another challenge, which forced us to run a separate simulation optimization process that uses standard violation function as a substitute function for FRVS to assess the output sensitivity to the inputs. According to existing decision makers and status, the conflict did not lead to a single stable equilibrium; thus, an arbitrator interferes with the conflict resolution process and is forced to omit some of the available scenarios to achieve a unique stable equilibrium.

The current article includes four main sections: introduction, methodology, results, and conclusion. The introduction section comprises a summary of works done in the field of WLA. The methodology section has eight subsections. The proposed model is the first subsection of this part; four different objective functions are presented and discussed. The optimization algorithm, simulation-optimization model, sensitivity analysis are the next subsections. At the end of this section, the PAWN method – which is used for sensitivity analysis – conflict resolution and the case study are described in detail. The third section, which is named results, covers the simulation-optimization, conflict resolution, and sensitivity analysis results. At the end of this article, a conclusion will be presented.

## METHODOLOGY

In this section, four objective functions, optimization algorithm, uncertainty analysis, simulation-optimization, sensitivity analysis, PAWN method, and conflict resolution, are discussed, and finally, the case study used to apply the methodology on is described. As depicted in Figure 1, the proposed model to WLA is made up of three parts: simulation, optimization, and conflict resolution. In the first section, parameters and variables used to solve the problem are defined. These are collected from official reports and previous studies. It is vital to determine the acceptable range for data by using a literature review and field study for variables. Although considering the seasonal uncertainty and extreme events is noteworthy and helps to have a cost-effective allocation, the lack of precise data series leads to consider the worst situation that covers all plausible conditions. The worst situation for water quality issues is the lowest flow, defined by the 7Q10 method and the highest temperature.

Figure 1

Methodology flowchart made up of three main parts of the proposed methodology: simulation, optimization, and conflict resolution.

Figure 1

Methodology flowchart made up of three main parts of the proposed methodology: simulation, optimization, and conflict resolution.

The river water quality simulation will be run in the deterministic mode. To define sensitive parameters and variables, a sensitivity analysis of water quality simulation was carried out, and sensitive inputs reassessed to be at their correct, acceptable range. In the optimization section, first of all, a stochastic simulation-optimization will be run, and optimal removal percent defined. The undominated surface (Pareto Front) of solutions will be depicted based on the triple objective functions which try to be achieved. The generated Pareto Front presents the feasible space for the problem's answer, and each point on this surface that will be selected has equal worth in comparison to the other points on this surface. In the next step of this section, six scenarios will be chosen from 100 scenarios, which generate Pareto Front, and each of these points contains removal percent for each discharger. Finally, a stochastic simulation model will be run against 100 selected scenarios of removal percent to define the value for violating standards from chosen scenarios.

In the last section of the proposed methodology, a conflict resolution approach will be used to select the most acceptable scenario by stakeholders who are called players in this method. After defining players and their preferences based on the six selected scenarios, the impossible states and status quo will be set for the conflict resolution model. Stable equilibriums will be defined by analyzing the game. If there was no stable equilibrium, the third-party intervenes in the conflict as an arbitrator and imposes its preferences on the game to lead to a sole stable equilibrium and put an end to the conflict.

### Objective function

Three objectives have been considered to allocate waste load for this river in stochastic mode. In the deterministic mode, one of those objectives will be replaced by a fourth function based on the reasons that will be presented in the following parts: (1) minimizing amount of violation from the national standard of water bodies contamination concentration (standard violation) (Equation (1)), (2) minimizing total cost for wastewater treatment and violating standard penalties (Equations (10)–(13)), and (3) minimizing the inequity index (Equation (14)) are used in a deterministic way, and in addition to the first two, (4) minimizing the FRVS (Equation (15)) had been used instead of minimizing the violation of the national standard in the stochastic model. The equations used in the simulation process can be found in the next sections.

### Standard violation

In this article, the objective of shielding water bodies' quality is to retain the BOD level under a predefined standard – the Iranian standard of discharge to surface water. In order to estimate the allowable level of BOD for dischargers, it is essential to use fate and transport equations for BOD to observe predefined criteria under different scenarios. Equations had been applied to estimate the decomposition of BOD and change in the dissolved oxygen (DO) level in the river system. Based on these equations, the standard violation will be estimated by Equations (6) and (9) will be used to calculate DO at the intended location.
(1)
(2)
Subjected to:
(3)
(4)
(5)
(6)
(7)
(8)
(9)
In the abovementioned equations, is the standard violation at point j that is set at the end of the river path. It is a function of waste load removal percent for each discharger, waste load weight , flow , temperature , reaeration and decomposition rates as well as national standards of water quality . Cost is a function of the fee paid by dischargers to treat their wastewater and a penalty relating to violating water quality standards.

is the concentration of oxidizable matter, which here is BOD in mg/L (expressed in oxygen level), is the initial concentration of oxidizable matter in mg/L at the upstream, is the total removal rate in , which is made of decomposition rate (kd) and settling rate (ks). ka is reareation rate (), is the concentration of saturated DO (mg/L), o is the concentration of DO in the river (mg/L), D is the oxygen deficit (mg/L), is the initial oxygen deficit (mg/L), x is the distance from the previous discharger (m), u is contamination velocity (m/s) and and are DO and BOD's standards established by Iran National Standards Organization (INSO) in mg/L. Due to the standard requirement to comply BOD and DO levels at checkpoints simultaneously, the proposed model always checks whether the DO concentration at checkpoint satisfies this requirement or not. Based on INSO, the DO concentration in checkpoints should be more than 5 mg/L.

### Cost function

One of the crucial parts of the WLA in a river system is defining a cost function. Based on previous works carried ou in Iran, in this paper, two distinct cost functions are used for municipal wastewater treatment and industrial wastewater treatment facilities (Jamshidi et al. 2013; Saberi & Niksokhan 2017). They used the following equations to estimate the treatment cost of municipal and industrial wastewater treatment. It has to be considered that for each discharger (PS or MWW), treatment cost and penalty will be calculated and used in the conflict resolution phase. In contrast, the total PS and MWW's costs and penalties have been used in the optimization process.
(10)
(11)
(12)

In the abovementioned equations, is the cost pay for removal percent treatment at the ith discharger, is the penalty applied to this discharger for violating the standard at the checkpoint. The second equation (Equation (11)) is used to estimate the cost for municipal wastewater treatment and the third (Equation (12)) belongs to industries cost to treat their wastewater. In these two equations, ci is the total payment for the treatment of one cubic meter of wastewater in one year (y) in USD.

Estalaki (2010) estimated a penalty function (Equation (13)) based on the replacement cost method and assumed five different combinations of BOD level that at each of which, the BOD level at checkpoint exceeded the allowable amount in the regulation, and eventually estimated the penalty function by using regression. In the current paper, the penalty divides equally among all of the dischargers, similar to Saberi & Niksokhan (2017). In the following equation (Equation (13)) v (mg/L) is a representative for the amount of the violating BOD standard, which based on the Iran national standard equals 5 mg/L in water bodies, and the penalty cost is determined in USD per year.
(13)

### Inequity index

Burn & Yulianti (2001) used the inequity concept, which is a problematic issue in WLA and is the result of a desire for fairness in the distribution of the treatment level and the associated costs among dischargers. They used the inequity index (Equation (14)) to take fairness into account. A decrease in the amount of this index indicates the improvement of fairness among dischargers, and a rise shows the decline of fairness. In this relation, NS is the number of dischargers, and are removal percent and input waste load of the ith discharger, respectively. can be defined by the multiplying flows (q) of the ith discharger by its BOD concentration. and are the average removal level and the average input waste load for all of the NS dischargers.
(14)

### Fuzzy risk of violating standard

Risk assessment is used extensively in water resources management fields such as water shortage (Feng & Huang 2008; Buurman et al. 2020), supply chain (Schaefer et al. 2019), drought (Shahid & Behrawan 2008; Chou et al. 2019), climate change (Dessai & Hulme 2007; Ronco et al. 2017), quality management (Mujumdar & Sasikumar 2002; Pullan et al. 2016; Lin et al. 2018; Liu et al. 2018), and waste load allocation (Rehana & Mujumdar 2009).

In the process of allocating waste load among dischargers in a river system under the deterministic situation, the first objective function was the standard violation, vj. To consider the uncertainty of variables and parameters in stochastic mode, instead of vj, the FRVS concept (Ghosh & Mujumdar 2006) was used for risk minimization in water quality control of the river system at the checkpoint. In their proposed equation (Equation (15)), is the fuzzy membership function of water quality related to the concentration of water quality index or c, here is BOD, and are the minimum and maximum concentration levels of BOD, respectively. indicates the probability of occurrence of concentration of c mg/L at the checkpoint and perch between zero to one. Since the standard for BOD concentration at the checkpoint is 5 mg/L, Figure 2 depicts the proposed fuzzy membership function and (Equation (16)) indicates its equation. If the BOD concentration at checkpoint exceeds than 8 mg/L, the fuzzy membership number is one. If BOD is less than 5 mg/L, then is zero. In the case where the BOD concentration at the checkpoint is in the range of 5–8, the will be between zero and one. The FRVS value declares the risk that BOD or any contamination level proceeds at the predefined level at the checkpoint due to the waste load discharged to the system. FRVS with a low value presents a lower risk for water quality. While zero FRVS is representative of a situation in which the river will not experience any violating from the predefined standard, FRVS equals 1 indicates that the situation at which the BOD level exceeds the allowable level, is inevitable. Therefore, to decrease this risk, either the probability of violating the standard or the level of this violation should be reduced.
(15)
(16)
Figure 2

Fuzzy membership function for water quality: is fuzzy membership number and range zero to one, and BOD is the concentration of the biochemical oxygen demand parameter at the checkpoint.

Figure 2

Fuzzy membership function for water quality: is fuzzy membership number and range zero to one, and BOD is the concentration of the biochemical oxygen demand parameter at the checkpoint.

### Optimization algorithm

The use of evolutionary algorithms for multi-objective optimization problems has been increased and has led to the development of various algorithms; NSGA-II (Deb et al. 2002), PAES (Knowles & Corne 2000), microGA (Coello & Pulido 2001), and MOPSO (Coello et al. 2004) are some of this evolutionary multi-objective optimization algorithms.

Particle swarm optimization (PSO) is a single objective optimization algorithm inspired by the choreography of a bird flock. It was developed for multi-objective issues and has been used in a wide range of problems such as multipurpose multi-reservoir operations (Fallah-Mehdipour et al. 2011), the study of ship's principal parameters (Lei 2011), chemical engineering (Taha et al. 2011), and textiles (Taheri et al. 2013). It is essential to use an optimization algorithm coupled with the simulation model to minimize the total cost of the system as well as another two objectives. The MOPSO model had been used since it is easy to implement and use, and its efficiency had been proven in some of the studies which were carried out in different fields (Coello et al. 2004; Goudos & Sahalos 2006; Durillo et al. 2009). The structure and principle of the MOPSO algorithm can be found in Coello et al. (2004).

### Uncertainty analysis

Outcomes or events which cannot be predicted with certainty are often called risky or uncertain. The term ‘risk’ is usually reserved to describe situations where probabilities are available to explain the likelihood of various events or outcomes. If the probabilities of multiple events or outcomes cannot be quantified, or if the events themselves are unpredictable, some would say the problem is one of uncertainty, and not of risk (Loucks et al. 2005). Uncertainty and associated terms such as error, risk, and ignorance are defined and interpreted by different authors, see Walker et al. (2003) for a review.

Refsgaard et al. (2007) adopted a subjective interpretation of uncertainty based on Klauer & Brown's (2004) definition in which the degree of confidence that a decision maker has about possible outcomes and/or probabilities of these outcomes is the central focus. Thus, according to the definition, a person is uncertain if s/he lacks confidence about the specific outcomes of an event. Reasons for this lack of confidence may include a judgment of the information as incomplete, blurred, inaccurate, unreliable, inconclusive, or potentially false. Similarly, a person is certain if s/he is confident about the outcomes of an event. It is possible that a person feels certain but has misjudged the information (i.e. his/her judgment is wrong).

Loucks et al. (2005) described two important classifications for uncertainty, where the first one divides uncertainty to knowledge uncertainty, decision uncertainty, and natural variability. In the second classification, they defined them as follows:

1. Informational uncertainties, which are made of imprecision in specifying the boundary and initial conditions as well as imprecision in measuring observed output variable values.

2. Model uncertainties, which are made of the uncertain model structure and parameter values, as well as the variability of observed input and output values over a region smaller than the spatial scale of the model. The variability of observed model input and output values within a period shorter than the time scale of the model and errors in linking models of different spatial and temporal scales are also made model uncertainties.

3. Numerical errors or errors in the model solution algorithm are specific sources of uncertainty.

While uncertainty analysis examines the lack of knowledge about the real value of model parameters, sensitivity analysis (SA) attempts to distinguish the ordinal or cardinal rank of variables' effects on changing output functions. Uncertainty can sometimes be reduced through further study, and by collecting additional data. Uncertainty analysis tries to determine the characteristics of different model's inputs and output distributions as well as functions of those random output variables that are performance indicators or measures. There are many tools and methods developed to take uncertainty into account. Refsgaard et al. (2007) listed 14 different methodologies and tools suitable for supporting uncertainty assessment and their guidance to the applicability. Data Uncertainty Engine (DUE) (Refsgaard et al. 2005), Error Propagation Equations (e.g. Mandel (1984)), Expert Elicitation (e.g. Spetzler & Stael von Holstein (1975)), Inverse Modelling (parameter estimation, freeware like PEST (Doherty 2003) and UCODE (Poeter & Hill 1998)), Monte Carlo Analysis (free software packages like @risk (Corporation 2009) and SimLab (Saltelli et al. 2004)), NUSAP (Funtowicz & Ravetz 1990; Van der Sluijs et al. 2004), and Uncertainty Matrix (Janssen et al. 2003; Walker et al. 2003) are among these methodologies and tools.

This article used SAFE Toolbox (Pianosi et al. 2015) and embedded methods for sensitivity analysis to take uncertainty modeling into account. Since the first phase of sensitivity and uncertainty analysis are both generating a large number of input from their acceptable range and their second phase is evaluating the generated inputs against fitness function (objective function), two initial steps of PAWN method code were used to generate input variables and parameters in their acceptable range and evaluate them. At the final phase, the output generated through the objective function calculated by the previous step will be used to generate three objective functions distributions. The average statistic of cost and inequity index distribution, as well as the FRVS as a sole value, will be assessed in the optimization process.

### Simulation-optimization

Based on the suggestion given by Nikoo et al. (2012) to take uncertainty into account in the simulation model, Q, BOD concentration, and temperature of upstream, tributaries and dischargers, as well as Ka and kr, which include 76 parameters and variables, will be considered.

A set of 100 series of random are generated, i.e. 100 different values for , , , …, and in the range 0–1. A set of 76 input parameters and variables will be generated randomly from their acceptable range to calculate three output values for cost, inequity index, and standard violation. Therefore, in the initial phase, there are 100 undominated solutions in the repository, each of which represents a series of (, , , …, and ) and have three values for the triple abovementioned objective function.

In the second phase, a new set of 100 series of is generated randomly and evaluated against objective functions by the MOPSO algorithm. If some of these new series of dominate any of the past 100 series, it will substitute the old and dominated ones. To take the uncertainty of variables and parameters into account in the simulation-optimization process for each set of 100 different , 200 different series of variable and input parameters were generated by the Latin-Hypercube method based on a normal distribution from their allowable space. Thus, for 200 different sets of inputs, which were generated by Latin-Hypercube and 100 sets of that were generated by MOPSO in the previous step, the simulation method gives three sets of 200 output for each of series, so for 100 different there are three distributions, each of which has a set of 200 numbers.

These 200 numbers lead to a probability distribution for each of the triple objectives. To assess each set of 100 with the previous set, which is in the repository, the MOPSO compares the average of cost and inequity index and the FRVS of 100 new with the prior set of in the repository. This process will repeat 500 times to achieve 100 different sets of undominated .

Eventually, at the end of 100,000 runs of the simulation-optimization process, an undominated surface will be generated that is made of 100 sets of undominated . Each of these points has a coordinate of three values: cost, inequity index, and FRVS. The first and second values are the average of a distribution made of 200 values. It is essential to indicate that each of these 100 points on the undominated surface, which is a Pareto Front, has the same worthiness, and all of them are optimal based on three objective functions.

One may use MCDM to select some of these optimal values, but authors use game theory to choose one distinct point which leads to a set of , those twelve dischargers are required to treat their wastewater up to this level and then empty to the river.

### Sensitivity analysis

Sensitivity analysis attempts to determine changes in outcomes values that are caused by changes in inputs. In other words, sensitivity analysis aims to comprehend the extent of change in output related to the change in inputs. It will evaluate the importance of uncertainty or inaccuracy of the model's inputs for the modeling or decision-making process. Sensitivity analysis could determine the change in system's optimal efficiency regarding the change in different parameters and also could determine how an optimal decision would change due to changes in the levels of resource constraints or desired outputs. Various techniques are developed to determine the sensitivity level of a model's outputs regarding a variety of input values. Most methods examine the effect of changing the value of a single parameter or input variable by taking into account the remaining parameters or inputs at a constant value. Sensitivity analysis could also be developed to investigate the composite effects of multiple sources of error (Loucks et al. 2005).

Sensitivity analysis is made up of three fundamental phases, which can be summarized as: (1) generation of input samples through their acceptable range; (2) model assessment for sets of the generated input value; and (3) post-processing to analysis input and output samples to estimate the sensitivity index (Pianosi & Wagener 2015).

Factor prioritization (FP), factor fixing (FF) or screening, and factor mapping (FM) are three primary purposes of the sensitivity analysis. The relative ranking of the model inputs contribution to output uncertainty is the aim of FP, determining the inputs that have no contribution to output uncertainty is the aim of FF, and discovering the regions in the inputs space that create specific output values is the aim of FM (Saltelli et al. 2008). FP and FF in sensitivity analysis are carried out by using the sensitivity index – a synthetic index quantifies the relative contribution to output uncertainty from each input. In a case where a parameter's sensitivity index value is zero, this parameter does not have any contribution to the output uncertainty, and wherever this index obtains a more significant value, expected parameters have a more substantial contribution to the output sensitivity (Pianosi & Wagener 2015).

Four main categories can be expressed for the sensitivity analysis as: (1) local sensitivity analysis and global sensitivity analysis; (2) mathematical, statistical and graphical; (3) screening and refined; and (4) qualitative and quantitative. Song et al. (2015) compiled a summary table of these basic categories for the sensitivity analysis methods, which contain descriptions, characteristics, and application cases of them as well as summary definitions of sensitivity analysis in the different fields. In a regional sensitivity analysis, the uncertainty caused by the change in input values around a specific point such as the average will be defined. General sensitivity analysis will determine the change in input values within the entire feasible space; this definition is also correct if used for outputs, which means it could investigate either entire feasible space or a specific area of an output (Pianosi & Wagener 2015).

Global sensitivity analysis methods are usually recommended in hydrological modeling applications because they have certain advantages compared with the local sensitivity analysis methods. The elementary effects test (EET, or Morris method (Morris 1991), regional sensitivity analysis (RSA, Spear & Hornberger 1980; Wagener & Kollat 2007), variance-based sensitivity analysis (VBSA, or Sobol’ method, e.g. Saltelli et al. (2008)), the Fourier amplitude sensitivity test (FAST by Cukier et al. (1973)), dynamic identifiability analysis (DYNIA by Wagener et al. (2003)) and a novel density-based sensitivity method, PAWN, by Pianosi & Wagener (2015) are embedded in the SAFE Toolbox (Pianosi et al. 2015).

Global, quantitative, model-independent, unconditional on any assumed input value, easy to interpret, easy to compute, stable, and moment independent are those features which a sound sensitivity analysis index should satisfy (Liu & Homma 2009). While variance-based sensitivity indices satisfy all the above features, but the last one, the PAWN method, meets all features (Pianosi & Wagener 2015).

### PAWN method

In this article, the authors use the PAWN method (derived from the names of Pianosi and Wagener) to perform sensitivity analysis. PAWN is a method for global sensitivity analysis and its central idea is to characterize output distributions by their cumulative distribution functions (CDF), which are easier to derive than probability distribution functions (PDFs). This method uses the sensitivity index to define outputs' sensitivity to inputs variability. SAFE toolbox and PAWN can use different methods such as Latin-Hypercube (McKay et al. 1979) and one-factor-at-a-time at the first phase to generate random input series. PAWN makes an N*M matrix of input values within their acceptable space, which contains N different values for M model inputs. In the second phase, a series of N values for M different input parameters and variables will enter the simulation model to calculate a series of N values for Y outputs. In the final phase of the process, the model's output series will be assessed by a sensitivity analysis method such as Morris, Sobol’, PAWN, and so on to analyze the sensitivity (Pianosi et al. 2015).

PAWN measures the distances between conditional and unconditional CDF of outputs for a specific parameter in which this method is trying to estimate its sensitivity index. In the first phase of estimating this distance, selecting all parameter and variables to calculate CDF of the output y will be done in a stochastic manner, and every value in the region of the allowable space of inputs has the chance to be selected and used; this kind of CDF is called unconditional CDF. To create a conditional CDF, parameter x, the parameter which the model tries to measure output sensitivity relative to it must be set to a specific value and selection of the other parameters and variables will be stochastic (. As one of the fundamental properties that a ‘good’ global sensitivity index should satisfy is to be unconditional about any assumed input value (Liu & Homma 2009), PAWN uses n different values for x to estimate n conditional CDFs and address this issue (Pianosi & Wagener 2015).

As is evident in the third box of Figure 3, there is an infinite amount of distance between an unconditional and a conditional CDF; therefore, the PAWN uses the Kolmogorov–Smirnov (Kolmogorov 1933; Smirnov 1939) statistic ) to calculate this length (Equation (17)). Since there are n different conditional CDFs and each one has a specific distance to the conditional CDF, this method uses the maximum of the KS statistic set to declare the final result, which is the x sensitivity index (Equation (18)). Steps in the numerical implementation of the PAWN sensitivity index can be found in Figure 3:
(17)
(18)
Figure 3

Phases in the numerical implementation of the PAWN sensitivity index (Pianosi & Wagener 2015).

Figure 3

Phases in the numerical implementation of the PAWN sensitivity index (Pianosi & Wagener 2015).

### Conflict resolution

The existence of disagreement among stakeholders led to the development of methods to address conflicts in different fields such as social, systems engineering, psychology, and operation research. A mathematical approach to model the behavior of decision makers (DMs) – here called players – on occasions where the preferences and decisions of one player may affect the other's preference and decision, is known by the term ‘game theory’. Although each decision maker has some power over the outcome, no single decision maker has full control (Hipel & Fang 2005). Rationality is the fundamental assumption of this method and means players who are making the decisions are intelligent, hence when they are trying to achieve distinct objectives, they understand and take into account other decision makers' rationality and, accordingly, tune their preferences based on this point (Parrachino 2006).

According to the scope of problems, a variety of game theory methods have been developed and extended. Madani & Hipel (2011) depicted two general categories of games and their subcategories: (1) quantitative procedures, which include normal form, extensive form, and cooperative game theory; (2) nonquantitative approaches which include metagame analysis and this, in turn, includes conflict analysis and drama theory.

Conflict resolution, which is a subset of non-quantitative approaches, has a supportive history of methods such as non-cooperative game theory (Von Neumann & Morgenstern 1947), metagame analysis (Howard 1971), conflict analysis (Fraser & Hipel 1984), drama theory (Howard 1994a, 1994b), theory of moves (Brams 1994), theory of fuzzy moves (Kandel et al. 1998; Li et al. 2001), hypergame analysis and graph model for conflict resolution (Kilgour et al. 1987).

Conflict resolution has been used in different fields such as tribal issues (Pereira & Athparia 2017), air traffic management (Hwang et al. 2007; Hou et al. 2017), automated guided vehicle (Reveliotis 2000), public–private partnerships (Osei-Kyei et al. 2019) and so on. In the case of water resources management, this method was used in both interstate (Eleftheriadou & Mylopoulos 2008; Mianabadi et al. 2015; He et al. 2018; Tayia 2019) and intrastate disputes (Amakali 2005; Schlager & Heikkila 2014), specifically for waste allocation (Karamouz et al. 2006).

Kinsara et al. (2015) introduced the fundamental steps of GMCR as modeling and analysis (Figure 4). Determining the DMs, all possible options, feasible state, allowable state transition, and relative preferences are the first step's aims. On the other hand, ascertaining the possible equilibria or outcomes of the conflict under a common stability definition is the purpose of the analysis step. The difference in concepts of stability definition refers to the difference in human behavior in dealing with strategic risk and their insight into the future. Four common stability definitions are Nash stability (R) (Nash 1950, 1951), general meta-rationality (GMR), symmetric meta-rationality (SMR) (Howard 1971), and sequential stability (SEQ) (Fraser & Hipel 1984). Three models have been developed for GMCR: GMCR I, GMCR II, and GMCR+. While the GMCR I (Hipel et al. 1990) model lacks the interface, it allows for the implementation of the conflict modeling and analysis. GMCR II had an interface and allowed interpreting and analyzing the results. Due to the introduction of the inverse graph model for conflict resolution concept, the GMCR+ model was introduced by Kinsara et al. (2015).

Figure 4

The basic procedure of applying the GMCR methodology to a real-world conflict (Kinsara et al. 2015).

Figure 4

The basic procedure of applying the GMCR methodology to a real-world conflict (Kinsara et al. 2015).

To analyze conflict in the GMCR+ model, infeasible states must be omitted, preferences prioritized, the direct ranking should be done, and the reachability matrix created. Preferences prioritization is a process that can make a list of players' (DMs') preferences based on their payoffs from corresponding states. Weighing is used by the GMCR+ for preferred states of every single DM and also presents manual prioritizing in which users could sort DMs preferred states manually (Kinsara et al. 2015).

Sakamoto et al. (2005) showed that in the cases where a stable equilibrium is not achieved, third-party intervention could alter the game status. It can play different roles: changing preferences and priorities of other decision makers (coordinator), leads other players to exclude some of their preferences (the arbitrator), and arrange for relevant parties to come to the negotiation table. Zanjanian et al. (2018) indicated that the entrance of the judiciary as an arbitrator in Iran could lead to a stable equilibrium.

### Case study

The proposed methodology was applied to the Sefidrood river in Gilan province, in northern Iran, to assess its applicability and efficiency. This river is 114.95 km length from the Sefidrood Dam outlet to the river meadows at the Caspian Sea. The river is divided into nine reaches according to its hydraulic characteristics. Twelve industries and municipal wastewater are discharging to the river. In this study it is supposed that each discharger uses its facilities to treat wastewater. Figure 5 depicts the location of Sefidrood and three tributaries that join the mainstream of this river: Tootakbon, Tarikrood, and Zilakirood.

Figure 5

Case study: Sefidrood river and three tributaries discharging to it are located at Guilan Province in northern Iran.

Figure 5

Case study: Sefidrood river and three tributaries discharging to it are located at Guilan Province in northern Iran.

## RESULTS

This section provides the results of different processes and models that are used in the proposed methodology in three subsections. In the first subsection, simulation-optimization results are presented, then the results of sensitivity analysis are described, and ultimately, the final solution for the WLA is presented through a conflict resolution subsection.

### Simulation-optimization

One hundred undominated solutions, demonstrated in Figure 6(a), were obtained through the stochastic simulation-optimization model, which used MOPSO as an evolutionary optimization algorithm. Figure 6(b)–6(d) display the undominated solution in a 2D view from different perspectives. Six scenarios (Table 1) were selected based on the most distributed state according to the three axes, which are objective functions.

Table 1

Removal percent of twelve dischargers (D1–D12) in six scenarios (S1–S6) and their related outputs for the intended objective function

DischargerRemoval percent in selected scenarios
S1S2S3S4S5S6
D1 0.84 0.16 0.94 0.28 0.38 0.53
D2 0.06 0.81 0.91 0.66 0.28 0.69
D3 0.68 0.48 0.61 0.57 0.49 0.60
D4 0.35 0.75 0.92 0.67 0.19 0.80
D5 0.09 0.43 0.42 0.45 0.24 0.43
D6 0.57 0.39 0.41 0.36 0.39 0.40
D7 0.34 0.37 0.63 0.49 0.47 0.59
D8 0.52 0.29 0.51 0.38 0.57 0.40
D9 0.38 0.81 0.84 0.71 0.18 0.79
D10 0.09 0.35 0.16 0.37 0.51 0.25
D11 0.43 0.18 0.62 0.40 0.59 0.45
D12 0.25 0.51 0.81 0.61 0.62 0.77
Inequity index 17.68 9.63 11.78 11.19 14.65 13.11
Cost *106 ($/y) 1.59 2.36 2.47 1.95 1.30 1.91 FRVS 10.80 14.01 0.04 15.9 11.29 6.72 DischargerRemoval percent in selected scenarios S1S2S3S4S5S6 D1 0.84 0.16 0.94 0.28 0.38 0.53 D2 0.06 0.81 0.91 0.66 0.28 0.69 D3 0.68 0.48 0.61 0.57 0.49 0.60 D4 0.35 0.75 0.92 0.67 0.19 0.80 D5 0.09 0.43 0.42 0.45 0.24 0.43 D6 0.57 0.39 0.41 0.36 0.39 0.40 D7 0.34 0.37 0.63 0.49 0.47 0.59 D8 0.52 0.29 0.51 0.38 0.57 0.40 D9 0.38 0.81 0.84 0.71 0.18 0.79 D10 0.09 0.35 0.16 0.37 0.51 0.25 D11 0.43 0.18 0.62 0.40 0.59 0.45 D12 0.25 0.51 0.81 0.61 0.62 0.77 Inequity index 17.68 9.63 11.78 11.19 14.65 13.11 Cost *106 ($/y) 1.59 2.36 2.47 1.95 1.30 1.91
FRVS 10.80 14.01 0.04 15.9 11.29 6.72
Figure 6

Pareto Front derived from the stochastic simulation-optimization process. (a) 3D view of the position of 100 undominated solutions (red dots) and six selected scenarios (black star named from S1–S6) for conflict resolution. (b) Solutions mentioned above in a 2D view, where its axes are inequity index versus FRVS. At the same time, (c) and (d) are a 2D view of FRVS versus cost and cost versus inequity index functions, respectively. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

Figure 6

Pareto Front derived from the stochastic simulation-optimization process. (a) 3D view of the position of 100 undominated solutions (red dots) and six selected scenarios (black star named from S1–S6) for conflict resolution. (b) Solutions mentioned above in a 2D view, where its axes are inequity index versus FRVS. At the same time, (c) and (d) are a 2D view of FRVS versus cost and cost versus inequity index functions, respectively. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

To obtain a comprehensive view of the alteration of the outputs of objectives functions, Figure 7 depicts the normalized values of the outputs in the range of zero to one. In this chart, the decline in the amount of one objective does not necessarily lead to a decrease in the amount of the other. For the S1 scenario, the inequity index is at the highest level, while FRVS and cost values are around the median in their allowable ranges. Cost function has the highest level in the S3 scenario, while this scenario has the lowest value in its acceptable range for FRVS and a low value for the inequity index. The S4 scenario, which is representative of a highly risky situation, has a low inequity index and a median cost. Nevertheless, the decline of FRVS will lead to an increase in the system's cost.

Figure 7

Normalized values of the outputs for objective functions: Blue pillars represent the inequity index, striped pillars represents cost, and gray pillars show FRVS. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

Figure 7

Normalized values of the outputs for objective functions: Blue pillars represent the inequity index, striped pillars represents cost, and gray pillars show FRVS. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

### Sensitivity analysis

In the current study, all 76 parameters and variables used in the stochastic simulation optimization process were used to be assessed in the sensitivity analysis. A critical point that should be considered is that instead of analyzing the FRVS, authors used standard violation as an objective function to assess its sensitivity; the reason is related to the nature of these two different functions, despite their similarity. The second phase of sensitivity analysis leads to the generation of a series of N values that could be used to make a CDF or PDF, hence comparing the different CDFs or PDFs of outputs estimated from different input leads to carrying out the sensitivity analysis.

In uncertainty analysis, a considerable number of inputs will be entered into the model to have a large number of output values. These values in turn will be used for generating CDFs or PDFs of model outputs. The definition of FRVS (Equation (15)) indicates that the final answer of this function is a result of assessing the occurrence probability of different concentrations that could be measured at the checkpoint and is a scholar value instead of a series of numbers. One has to use the model output's PDF to calculate the fuzzy risk. In contrast, in the second phase of sensitivity analysis one has to reach a PDF or CDF, and this is not possible, hence instead of using FRVS, which is the final outcome of uncertainty analysis and is a scalar value, one has to use a function whose outcome of uncertainty analysis will be a series of scalar values instead of a single scalar value. Due to this point, a standard violation function, which has a type the same as p(c) in the FRVS (Equation (7)), had been used. The results of the uncertainty analysis are depicted in Figures 810. In these figures, Q, T, and BOD represent the flow, temperature, and biochemical oxygen demand of different dischargers flow to the riverine system, which includes industries, municipal wastewater treatment facilities, and the rivers tributaries: the upstream of Sefidrood, Tootakbon, Tarikrood, and Zilakirood. The subscripts of the above-mentioned parameters and variables refer to the number of each discharger; for instance, among them, Q1 and Q12 are the discharge of the first and twelfth dischargers, respectively, which can be industry or municipal wastewater dischargers. In these figures, u, ka, and kd are stream velocity, reareation rate, and deoxygenation rate in river's nine reaches, respectively.

Figure 8

Sensitivity index of 76 inputs for the standard violation function.

Figure 8

Sensitivity index of 76 inputs for the standard violation function.

Figure 9

Sensitivity index of 76 inputs for the cost function.

Figure 9

Sensitivity index of 76 inputs for the cost function.

Figure 10

Sensitivity index of 76 inputs for the inequity index function.

Figure 10

Sensitivity index of 76 inputs for the inequity index function.

The sensitivity analysis result (Figure 8) shows that standard violation has the most sensitivity to upstream flow (Qupstream), BOD concentration of the ninth discharger, BOD concentration of Tarikrood, BOD concentration of Zilakirood, and Tarikrood flow, respectively. Due to the direct effect of flow (Q) on the concentration of BOD in the riverine system, this objective function has the highest sensitivity to the river's flow, which contains upstream, Tarikrood, Zilakirood, and Tootakbon tributaries flow. The relatively high sensitivity of the standard violation to the concentration of the ninth discharger is due to the vicinity of this discharger to the checkpoint and its higher Q in contrast to other dischargers. At a constant concentration of the BOD, the higher the Q of a discharger leads to the entrance of the higher amount of waste load to the river.

In the case of the cost function, the results of sensitivity analysis (Figure 9) show that this function has the highest sensitivity to upstream flow (QUpstream), BOD concentration of Tarikrood tributary, BOD concentration of ninth discharger (BOD9), and its flow (Q9), respectively. The reason for this ordering is the same as the reasons mentioned above for standard violation.

Based on the sensitivity index driven from the PAWN method, which is depicted in Figure 10, the inequity index function has the highest sensitivity to the BOD concentration and flow of the ninth discharger's (BOD9 and Q9), the seventh discharger's flow (Q7), and the ninth discharger's flow (Q9), respectively. The order of the parameters sensitivity index indicates that after the previous parameters, the function has the most sensitivity to the different dischargers' Q and BOD concentrations in the next place. It should be noticed that the upstream discharge (QUpstream) has the least sensitivity index, among other parameters and variables. Since the definition of inequity (Equation (14)) is just dependent on the amount of waste load , and removal percent , such ordering of the sensitivity index can be justified.

According to the fact that the sensitivity index (the average of Ti index) is a dimensionless value, if the average of this index is calculated, it will reveal that the proposed methodology has the most sensitivity to upstream flow (QUpstream), the ninth discharger's BOD concentration and flow (BOD9 and Q9), and the seventh discharger's flow (Q7). It has also revealed that the methodology has less sensitivity to the reareation rate at the second reach of the river. An overall view of the sensitivity analysis for all 76 inputs is displayed in Figure 11.

Figure 11

Sensitivity index of all parameters and variables (76 inputs) for three objectives.

Figure 11

Sensitivity index of all parameters and variables (76 inputs) for three objectives.

To have a comprehensive view about the variability of output's CDF, Figures 12 and 13 show the variability of output's CDF based on the change in velocity of the seventh reach of the river (u7) and upstream flow (QUpstream) as the most and the least effective inputs on the standard violation function, respectively. In these figures, the red line indicates the unconditional CDF of output values while the gray lines depict the conditional CDF. A black and white spectrum exists at the right side of the graphs, which is an indicator of velocity in the seventh reach in Figure 12 and upstream flow (QUpstream) in Figure 13, both at the conditional mode. It can be understood from Figure 13 that as the upstream flow receives a more significant value at the conditional status, Q = 0.551 CMS, the standard violation will lead to the smallest value. This CDF, the first gray CDF of the graph from the left side, indicates the lowest variability of CDF for standard violation. On the other hand, the smallest amount of upstream discharge, Q = 0.1766 CMS, causes the most violations from a specified standard; the CDF for this conditional output is drawn up in black color at the end of the graph.

Figure 12

Variability of standard violation CDF against variability in the velocity of the seventh reach of the river. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

Figure 12

Variability of standard violation CDF against variability in the velocity of the seventh reach of the river. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

Figure 13

Variability of standard violation CDF against variability in upstream discharge. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

Figure 13

Variability of standard violation CDF against variability in upstream discharge. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

### Conflict resolution

Two of the intended objective functions, cost, and FRVS, are in contrast to each other, so that achieving better performance for one of them requires getting away from the best performance for another one. Reducing the FRVS needs higher treatment percentage, and hence increases dischargers' costs, therefore reducing the system's total cost will lead to a higher FRVS amount.

Two different entities rule water management issues in Iran. The Ministry of Energy (MOE) is responsible for supplying water to all consumers as well as providing wastewater facilities. The Department of Environment (DOE) is liable for protecting water bodies' environmental quality. Regulations issued by DOE are not mandatory for MOE because of institutional problems, so there is a conflict among these parties on water quality. The Municipal Waste Water (MWW), as an affiliation of MOE, specifically focuses on wastewater affairs, tries to reduce its costs of construction and operation of wastewater treatment facilities and thus decrease the treating of wastewater's cost. Whereas the Iranian government is responsible for building and operating wastewater treatment facilities, the DOE and its affiliate, MWW, are its representatives.

On the other hand, DOE only considers declining waste load discharged from municipal and private wastewater treatment facilities to the water bodies to increase water quality, specifically in checkpoints.

The third player in this conflict is the private sector (PS), which is representative of a coalition of industries surrounding the river and owned by the non-governmental sector. This player or DM aims to maximize its profit or minimize its cost; thus, it has specific preferences, which would be in contrast to other players' preferences. Table 2 reveals the player's payoffs for various scenarios obtained from the simulation-optimization section and summarized in Table 1. The differences among scenarios' payoffs for each one of the players show the differences in their preferences.

Table 2

Players' payoffs for six scenarios; MWW and PS look at their cost and penalty ($/y) and DOE focus on FRVS ScenariosDMs' payoffs MWW ($/y)PS ($/y)DOE (FRVS)Inequity index S1 1,304,698 286,711 10.80 17.68 S2 2,279,411 77,457 14.01 9.63 S3 2,291,446 182,080 0.04 11.78 S4 1,860,174 89,264 15.90 11.19 S5 1,137,180 160,932 11.29 14.65 S6 1,797,812 110,790 6.72 13.11 ScenariosDMs' payoffs MWW ($/y)PS (\$/y)DOE (FRVS)Inequity index
S1 1,304,698 286,711 10.80 17.68
S2 2,279,411 77,457 14.01 9.63
S3 2,291,446 182,080 0.04 11.78
S4 1,860,174 89,264 15.90 11.19
S5 1,137,180 160,932 11.29 14.65
S6 1,797,812 110,790 6.72 13.11

Players' preferences, which are provided in Table 3, based on their payoffs and interests, introduced to the GMCR+ model. DOE tries to reduce the FRVS, so its choices are sorted based on the minimum amount of this objective up to the maximum ones. MWW and PS both also try to reduce their costs; thus, they will sort their preferences based on the cost that they must pay for scenarios.

Table 3

Players’ preferences

Preference no.Players (DMs)
DOEMWWPS
S3 S5 S2
S6 S1 S4
S1 S4 S6
S5 S6 S5
S4 S2 S3
S2 S3 S1
Preference no.Players (DMs)
DOEMWWPS
S3 S5 S2
S6 S1 S4
S1 S4 S6
S5 S6 S5
S4 S2 S3
S2 S3 S1

In the GMCR+ model, a player should decide about an option clearly and answer explicitly about the selection of an option: Yes or No. The combination of scenarios (options) selected by players creates a state. For instance, one option can be S3, S5, and S2 which means the DOE selects scenario number S3, MWW selects S4, and eventually, PS chooses S2. The presence of three players and six options for each one lead to states in this conflict. Defining impossible states reduced the entire states to 216. One player does not decide about a specific option in a state (null status), one player does not select any options in a state so that it selects No for all options and eventually, selecting more than one choice in a state are impossible states.

Analyzing conflict stability leads to 125 stable equilibrium states, in which the status quo is one of them. In the status quo, which is representative of the worst situation of the river in terms of quality, dischargers treat waste load as little as possible to pay the lowest cost on waste load treatment. It should be noted that GMCR+ finds stable equilibriums or solutions only based on equilibrium definitions and does not consider the status quo. After finding stable equilibriums, the status quo should be selected among states that exist in the model and detect final stable equilibrium or ultimate solution of conflict by tracing a tree graph. In the status quo of this conflict, MWW will pick option S5; PS will select S2, and it is inevitable for DOE to choose S4. The GMCR+ model numbered this state as 64. Analyzing the extracted graph (Figure 14) from the model reveals that, based on the current status quo, the conflict will not obtain a single stable equilibrium.

Figure 14

Graph of players movement from the initial status quo to states which have stable equilibrium (output from GMCR+ model). Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

Figure 14

Graph of players movement from the initial status quo to states which have stable equilibrium (output from GMCR+ model). Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/hydro.2020.181.

In this figure, blue lines show the DOE movements, green lines depict WWC movements, and yellow lines indicate PS movements. Since each player prefers to move to a state with more payoff for itself in comparison to the status quo, DOE can change its selected state from state number 64 to one numbered 61, 62, 63, 65, or 66. WWC will prefer to choose one of the state numbers 40, 46, 52, 58, or 70 to achieve more payoff, while the best payoff for PS will be to move to one of the states numbered 28, 100, 136, 172, or 208. By changing its selection, each player will motivate two other players to improve their choice; thus, the current status quo alongside the players' preferences will not lead to a common option selection by two other players, and a stable equilibrium will not be achieved. For instance, PS comprehends that selecting state number 58 is in its favor; thus, select this state, but two other players try to choose other states, which has more payoff for them in comparison to status number 58.

The intervention of a third party was assessed to comprehend whether it is useful to resolve the conflict. Hence the judiciary entered the conflict as an arbitrator who can force players not to select specific options, including the options with the highest FRVS and the most significant inequity index. The judiciary will ban players from selecting S2 and S4 as the options, with the highest amount of FRVS and S1 as the option with the highest inequity index.

Removing the afore-mentioned triple options causes a change in the status quo, and selected options will move to the nearest ones. In this state, which is representative of the worst situation and numbered 23 by GMCR + , MWT still select S5, DOE selected option will move from S4 to S5, and finally, PS will select S6.

Analyzing the new conflict results from the entrance of arbitrator to previous conflicts shows that all possible stable equilibriums will decline to nine items. By considering the new status quo, GMCR+ outputs prove PS will move from S6 to S5, and consequently, two other players comprehending moving from S5 to another option will not increase their payoffs, and thus they do not change their selection. This state is named number 14 by the model. Therefore, one movement from status number 23 to status number 14 by PS will lead to a stable equilibrium, which is the final resolution (Figure 15).Table 4 shows the twelve dischargers' removal percent for stable equilibrium, which is the final WLA based on the predefined objective functions and triple decision makers' preferences.

Table 4

Removal percent for all dischargers in the selected scenario for waste load allocation

DischargerD1D2D3D4D5D6D7D8D9D10D11D12
Removal percent (0.38 0.28 0.49 0.19 0.24 0.39 0.47 0.57 0.18 0.51 0.59 0.62
DischargerD1D2D3D4D5D6D7D8D9D10D11D12
Removal percent (0.38 0.28 0.49 0.19 0.24 0.39 0.47 0.57 0.18 0.51 0.59 0.62
Figure 15

Graph of final stable equilibrium for conflict resolution; PS move its position from state 23 to 14, and stable equilibrium will be achieved (output from GMCR+ model).

Figure 15

Graph of final stable equilibrium for conflict resolution; PS move its position from state 23 to 14, and stable equilibrium will be achieved (output from GMCR+ model).

S5 will be the option in which the entrance of an arbitrator causes it to be accepted by all parties of the conflict. This option leads to a total cost of 1,298,112 USD per year and MWW should pay 1,137,180 USD, and PS should pay 160,932 USD. Selecting this option also leads to an inequity value of 14.65 and FRVS of 11.29.

## REFERENCES

Afshar
A.
Masoumi
F.
Sandoval Solis
S.
2018
Developing a reliability-based waste load allocation strategy for river-reservoir systems
.
J. Water Resour. Plan. Manage.
144
(
9
),
04018052
.
Allam
A.
Tawfik
A.
Yoshimura
C.
Fleifle
A.
2016
Multi-objective models of waste load allocation toward a sustainable reuse of drainage water in irrigation
.
Environ. Sci. Pollut. Res.
23
(
12
),
11,823
11,834
.
Amakali
M.
2005
Intra-state conflict resolution between local communities and central governments – Namibia Case
. In:
Value of Water–Different Approaches in Transboundary Water Management
.
IHP/HWRP-Sekretariat
,
Koblenz
,
Germany
, pp.
103
108
.
Brams
S. J.
1994
Theory of Moves
.
Cambridge University Press
,
Cambridge
.
Breetz
H. L.
Fisher-Vanden
K.
Garzon
L.
Jacobs
H.
Kroetz
K.
Terry
R.
2004
Water Quality Trading and Offset Initiatives in the US: A Comprehensive Survey
.
Dartmouth College and the Rockefeller Center for the US Environmental Protection Agency, Dartmouth College
,
Hanover, New Hampshire
.
Burn
D. H.
Lence
B. J.
1992
Comparison of optimization formulations for waste-load allocations
.
J. Environ. Eng.
118
(
4
),
597
612
.
Burn
D. H.
Yulianti
J. S.
2001
.
J. Water Resour. Plan. Manage.
127
(
2
),
121
129
.
Buurman
J.
Bui
D. D.
Du
L. T. T.
2020
Drought risk assessment in Vietnamese communities using household survey information
.
Int. J. Water Resour. Dev.
36
(
1
),
88
105
.
R. A.
Miller
A. C.
McDonnell
A. J.
1981
Analysis of waste load allocation procedures
.
JAWRA J. Am. Water Resour. Assoc.
17
(
5
),
760
766
.
Chou
J.
Xian
T.
Zhao
R.
Xu
Y.
Yang
F.
Sun
M.
2019
Drought risk assessment and estimation in vulnerable eco-regions of China: under the background of climate change
.
Sustainability
11
(
16
),
4463
.
Coello
C. A. C. C.
Pulido
G. T.
2001
A micro-genetic algorithm for multiobjective optimization
. In:
Proc., International Conference on Evolutionary Multi-Criterion Optimization
.
Springer
,
Berlin, Heidelberg
, pp.
126
140
.
Coello
C. A. C.
Pulido
G. T.
Lechuga
M. S.
2004
Handling multiple objectives with particle swarm optimization
.
IEEE Trans. Evol. Comput.
8
(
3
),
256
279
.
Corporation
P.
2009
Guide to Using@ RISK.: Risk Analysis and Simulation Add-in for Microsoft Excel
.
,
USA
.
Cukier
R.
Fortuin
C.
Shuler
K. E.
Petschek
A.
Schaibly
J.
1973
Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I Theory
.
J. Chem. Phys.
59
(
8
),
3873
3878
.
Deb
K.
Pratap
A.
Agarwal
S.
Meyarivan
T.
2002
A fast and elitist multiobjective genetic algorithm: NSGA-II
.
IEEE Trans. Evol. Comput.
6
(
2
),
182
197
.
Doherty
J.
2003
Ground water model calibration using pilot points and regularization
.
Groundwater
41
(
2
),
170
177
.
Durillo
J. J.
García-Nieto
J.
Nebro
A. J.
Coello
C. A. C.
Luna
F.
Alba
E.
2009
Multi-objective particle swarm optimizers: An experimental comparison
. In:
Proceedings of the International Conference on Evolutionary Multi-Criterion Optimization
.
Springer
,
Berlin, Heidelberg
, pp.
495
509
.
E.
Mylopoulos
Y.
2008
Game theoretical approach to conflict resolution in transboundary water resources management
.
J. Water Resour. Plan. Manage.
134
(
5
),
466
473
.
Estalaki
S. M.
2010
Codifying Policies for Qualitative Management of Rivers Using Evolutionary Games
.
University of Tehran
,
Tehran
,
Iran
.
Fallah-Mehdipour
E.
O. B.
Mariño
M.
2011
MOPSO algorithm and its application in multipurpose multireservoir operations
.
J. Hydroinform.
13
(
4
),
794
811
.
Feizi Ashtiani
E.
Niksokhan
M.
Ardestani
M.
2015
Multi-objective waste load allocation in river system by MOPSO algorithm
.
Int. J. Environ. Res.
9
(
1
),
69
76
.
Fraser
N. M.
Hipel
K. W.
1984
Conflict Analysis: Models and Resolutions
.
Elsevier-Science Ltd
,
North-Holland
.
Funtowicz
S. O.
Ravetz
J. R.
1990
Uncertainty and Quality in Science for Policy
.
,
The Netherlands
.
Ghosh
S.
Mujumdar
P.
2006
Risk minimization in water quality control problems of a river system
.
29
(
3
),
458
470
.
Ghosh
S.
Mujumdar
P.
2010
Fuzzy waste load allocation model: a multiobjective approach
.
J. Hydroinform.
12
(
1
),
83
96
.
Goudos
S.
Sahalos
J.
2006
Microwave absorber optimal design using multi-objective particle swarm optimization
.
Microw. Opt. Technol. Lett.
48
(
8
),
1553
1558
.
Hathhorn
W. E.
Tung
Y.-K.
1989
Bi-objective analysis of waste load allocation using fuzzy linear programming
.
Water Resour. Manage.
3
(
4
),
243
257
.
He
S.
Kilgour
D.
Hipel
K.
2018
A basic hierarchical graph model for conflict resolution with weighted preference
.
J. Environ. Inform.
31
(
1
),
15
29
.
Hipel
K.
Fang
L.
2005
Multiple participant decision making in societal and technological systems
. In:
Proc., Systems and Human Science, for Safety, Security, and Dependability: Seleced Papers of the 1st International Symposium SSR2003
,
November 2003
,
Osaka, Japan
.
Elsevier
,
The Netherlands
.
Hipel
K. W.
Fang
L.
Kilgour
D. M.
1990
A formal analysis of the Canada-US softwood lumber dispute
.
Eur. J. Oper. Res.
46
(
2
),
235
246
.
Hou
X.
Liu
Y.
Sourina
O.
Chen
C.-H.
Mueller-Wittig
W.
Ang
W. T.
2017
EEG-based human factors evaluation of conflict resolution aid and tactile user interface in future air traffic control systems
. In:
Advances in Human Aspects of Transportation
(N. A. Stanton, S. Landry, G. Di Bucchianico & A. Vallicelli, eds).
Springer
,
Switzerland
, pp.
885
897
.
Howard
N.
1971
Paradoxes of Rationality: Games, Metagames, and Political Behavior
.
MIT Press
,
Cambridge, MA
.
Hwang
I.
Kim
J.
Tomlin
C.
2007
Protocol-based conflict resolution for air traffic control
.
Air Traffic Contr. Q.
15
(
1
),
1
34
.
Jamshidi
S.
G.
A.
2013
An assessment of using water quality trading to improve water quality management
. In:
Proc., 3rd International Conference on Environmental Management and Planning
,
Tehran, Iran (in Persian)
.
Janssen
P.
Petersen
A.
Van der Sluijs
J.
Risbey
J.
Ravetz
J.
2003
RIVM/MNP Guidance for Uncertainty Assessment and Communication: Quickscan Hints & Actions List
.
RIVM/MNP
.
Available from: www.nusap.net
.
Kandel
A.
Zhang
Y.
Borges
P.
1998
Fuzzy prisoner's dilemma
.
Fuzzy Econom. Rev.
3
(
1
),
3
20
.
Karamouz
M.
Akhbari
M.
Moridi
A.
Kerachian
R.
2006
A system dynamics-based conflict resolution model for river water quality management
.
Iran. J. Environ. Health Sci. Eng.
3
(
3
),
147
160
.
Kerachian
R.
Karamouz
M.
2005
Waste-load allocation model for seasonal river water quality management: application of sequential dynamic genetic algorithms
.
Sci. Iran.
12
(
2
),
117
130
.
Kilgour
D. M.
Hipel
K. W.
Fang
L.
1987
The graph model for conflicts
.
Automatica
23
(
1
),
41
55
.
Kinsara
R. A.
Petersons
O.
Hipel
K. W.
Kilgour
D. M.
2015
Advanced decision support for the graph model for conflict resolution
.
J. Decis. Syst.
24
(
2
),
117
145
.
Klauer
B.
Brown
J.
2004
Conceptualising imperfect knowledge in public decision-making: ignorance, uncertainty, error and risk situations
.
Environ. Res. Eng. Manage.
1
,
142
148
.
Knowles
J. D.
Corne
D. W.
2000
M-PAES: A memetic algorithm for multiobjective optimization
. In:
Proceedings of the 2000 Congress on Evolutionary Computation
.
IEEE
,
New York
, pp.
325
332
.
Kolmogorov
A.
1933
On the empirical determination of a distribution law
.
Inst. Ital. Attuari, Giorn.
4
,
83
91
(in Italian)
.
Lei
H.
2011
Application of multi-objective particle swarm optimization (MOPSO) in study of ship's principal parameters
.
J. Ship Mech.
7
,
784
790
.
Li
K. W.
Karray
F.
Hipel
K. W.
Kilgour
D. M.
2001
Fuzzy approaches to the game of chicken
.
IEEE Trans. Fuzzy Syst.
9
(
4
),
608
623
.
Lin
C.
Ma
R.
Wu
Z.
Xiong
J.
Min
M.
2018
Detection of the sensitive inflowing river indicators related to non-point source organic pollution: a case study of Taihu Lake
.
J. Environ. Inf.
32
,
98
111
.
Liu
Q.
Homma
T.
2009
A new computational method of a moment-independent uncertainty importance measure
.
Reliab. Eng. Syst. Sfty.
94
(
7
),
1205
1211
.
Liu
R.
Zhang
K.
Zhang
Z.
Borthwick
A. G.
2018
Watershed-scale environmental risk assessment of accidental water pollution: the case of Laoguan river, China
.
J. Environ. Inform.
31
(
2
),
87
96
.
Loucks
D. P.
Van Beek
E.
Stedinger
J. R.
Dijkman
J. P.
Villars
M. T.
2005
Water Resources Systems Planning and Management: an Introduction to Methods, Models and Applications
.
Unesco
,
Paris
.
K.
Hipel
K. W.
2011
Non-cooperative stability definitions for strategic analysis of generic water resources conflicts
.
Water Resour. Manage.
25
(
8
),
1949
1977
.
Mandel
J.
1984
The Statistical Analysis of Experimental Data
.
Interscience Publications
,
Washington, DC
.
McKay
M. D.
Beckman
R. J.
Conover
W. J.
1979
Comparison of three methods for selecting values of input variables in the analysis of output from a computer code
.
Technometrics
21
(
2
),
239
245
.
H.
Mostert
E.
van de Giesen
N.
2015
Trans-boundary river basin management: factors influencing the success or failure of international agreements
. In:
Conflict Resolution in Water Resources and Environmental Management
(K. W. Hipel, L. Fang, J. Cullmann & M. Bristow, eds).
Springer
,
Switzerland
, pp.
133
143
.
Mirchi
A.
Watkins
D.
Jr
K.
2010
Modeling for Watershed Planning, Management, and Decision Making
.
Nova Science Publishers
,
New York
.
Mohan
S.
Kumar
K. P.
2016
Waste load allocation using machine scheduling: model application
.
Environ. Process.
3
(
1
),
139
151
.
Morris
M. D.
1991
Factorial sampling plans for preliminary computational experiments
.
Technometrics
33
(
2
),
161
174
.
Mujumdar
P.
Sasikumar
K.
2002
A fuzzy risk approach for seasonal water quality management of a river system
.
Water Resour. Res.
38
(
1
),
5
.
Nash
J. F.
1950
Equilibrium points in n-person games
.
36
(
1
),
48
49
.
Nash
J.
1951
Non-cooperative games
.
Ann. Math.
54
(
2
),
286
295
.
Nikoo
M. R.
Kerachian
R.
Niksokhan
M. H.
2012
Equitable waste load allocation in rivers using fuzzy Bi-matrix games
.
Water Resour. Manage.
26
(
15
),
4539
4552
.
Niksokhan
M. H.
Kerachian
R.
Karamouz
M.
2009
A game theoretic approach for trading discharge permits in rivers
.
Water Sci. Technol.
60
(
3
),
793
804
.
Parrachino
I.
2006
Cooperative Game Theory and its Application to Natural, Environmental and Water Resource Issues
.
World Bank Publications
,
Washington, DC
.
Pereira
M.
Athparia
R.
2017
The resilience of tribal conflict resolution systems in Northeast India: a panoramic view
. In:
Legal Pluralism and Indian Democracy
(M. Pereira, B. Dutta & B. Kakati, eds).
Routledge
,
India
, pp.
122
148
.
Pianosi
F.
Sarrazin
F.
Wagener
T.
2015
A Matlab toolbox for global sensitivity analysis
.
Environ. Model. Softw.
70
,
80
85
.
Poeter
E. P.
Hill
M. C.
1998
Documentation of UCODE, A Computer Code for Universal Inverse Modeling
.
DIANE Publishing
,
Denver, CO
.
Pullan
S.
Whelan
M.
Rettino
J.
Filby
K.
Eyre
S.
Holman
I. P.
2016
Development and application of a catchment scale pesticide fate and transport model for use in drinking water risk assessment
.
Sci. Total Environ.
563
,
434
447
.
Qin
X.
Huang
G.
Chen
B.
Zhang
B.
2009
An interval-parameter waste-load-allocation model for river water quality management under uncertainty
.
Environ. Manage.
43
(
6
),
999
1012
.
Rafiee
M.
Lyon
S. W.
Zahraie
B.
Destouni
G.
N.
2017
Optimal wastewater loading under conflicting goals and technology limitations in a riverine system
.
Water Environ. Res.
89
(
3
),
211
220
.
Refsgaard
J. C.
Nilsson
B.
Brown
J.
Klauer
B.
Moore
R.
Bech
T.
Vurro
M.
Blind
M.
Castilla
G.
Tsanis
I.
2005
Harmonised techniques and representative river basin data for assessment and use of uncertainty information in integrated water management (HarmoniRiB)
.
Environ. Sci. Pol.
8
(
3
),
267
277
.
Refsgaard
J. C.
van der Sluijs
J. P.
Højberg
A. L.
Vanrolleghem
P. A.
2007
Uncertainty in the environmental modelling process – a framework and guidance
.
Environ. Model. Softw.
22
(
11
),
1543
1556
.
Rehana
S.
Mujumdar
P.
2009
An imprecise fuzzy risk approach for water quality management of a river system
.
J. Environ. Manage.
90
(
11
),
3653
3664
.
Reveliotis
S. A.
2000
Conflict resolution in AGV systems
.
IIE Trans.
32
(
7
),
647
659
.
Ronco
P.
Zennaro
F.
Torresan
S.
Critto
A.
Santini
M.
Trabucco
A.
Zollo
A.
Galluccio
G.
Marcomini
A.
2017
A risk assessment framework for irrigated agriculture under climate change
.
110
,
562
578
.
M.
Afshar
A.
2007
Waste load allocation modeling with fuzzy goals; simulation-optimization approach
.
Water Resour. Manage.
21
(
7
),
1207
1224
.
Saberi
L.
Niksokhan
M. H.
2017
Optimal waste load allocation using graph model for conflict resolution
.
Water Sci. Technol.
75
(
6
),
1512
1522
.
Sakamoto
M.
Hagihara
Y.
Hipel
K. W.
2005
Coordination process by a third party in the conflict between Bangladesh and India over regulation of the Ganges River
. In:
International Conference on Systems, Man and Cybernetics
.
IEEE
,
Waikoloa, HI
,
USA
, pp.
1119
1125
.
Saltelli
A.
Tarantola
S.
Campolongo
F.
Ratto
M.
2004
Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models
.
John Wiley & Sons
,
Chichester, UK
.
Saltelli
A.
Ratto
M.
Andres
T.
Campolongo
F.
Cariboni
J.
Gatelli
D.
Saisana
M.
Tarantola
S.
2008
Global Sensitivity Analysis: the Primer
.
John Wiley & Sons
,
Chichester, UK
.
Sarang
A.
Lence
B. J.
Shamsai
A.
2008
Multiple interactive pollutants in water quality trading
.
Environ. Manage.
42
(
4
),
620
646
.
Schaefer
T.
Udenio
M.
Quinn
S.
Fransoo
J. C.
2019
Water risk assessment in supply chains
.
J. Clean. Prod.
208
,
636
648
.
Schlager
E.
Heikkila
T.
2014
Water scarcity, conflict resolution, and adaptive governance in federal transboundary river basins
. In:
Federal Rivers: Managing Water in Multi-Layered Political Systems
(D. Garrick, G. R. M. Anderson, D. Connell & J. Pittock, eds).
Edward Elgar Publications
,
Cheltenham
,
UK
, p.
57
.
Shahid
S.
Behrawan
H.
2008
Drought risk assessment in the western part of Bangladesh
.
Nat. Hazards
46
(
3
),
391
413
.
Smirnov
N. V.
1939
On the estimation of the discrepancy between empirical curves of distribution for two independent samples
.
Bull. Math. Univ. Moscou
2
(
2
).
Spetzler
C. S.
Stael von Holstein
C.-A. S.
1975
Exceptional paper – probability encoding in decision analysis
.
Manage. Sci.
22
(
3
),
340
358
.
Van der Sluijs
J.
Janssen
P.
Petersen
A.
Kloprogge
P.
Risbey
J.
Tuinstra
W.
Ravetz
J.
2004
RIVM/MNP Guidance for Uncertainty Assessment and Communication: Tool Catalogue for Uncertainty Assessment
.
Utrecht University
,
The Netherlands
. .
Von Neumann
J.
Morgenstern
O.
1947
Theory of Games and Economic Behavior
, 2nd rev.
Princeton University Press
,
Princeton, NJ, USA
.
Wagener
T.
McIntyre
N.
Lees
M.
Wheater
H.
Gupta
H.
2003
Towards reduced uncertainty in conceptual rainfall-runoff modelling: dynamic identifiability analysis
.
Hydrol. Process.
17
(
2
),
455
476
.
Walker
W. E.
Harremoës
P.
Rotmans
J.
van der Sluijs
J. P.
van Asselt
M. B.
Janssen
P.
Krayer von Krauss
M. P.
2003
Defining uncertainty: a conceptual basis for uncertainty management in model-based decision support
.
Integr. Assess.
4
(
1
),
5
17
.
Yandamuri
S.
Srinivasan
K.
Murty Bhallamudi
S.
2006
Multiobjective optimal waste load allocation models for rivers using nondominated sorting genetic algorithm-II
.
J. Water Resour. Plan. Manage.
132
(
3
),
133
143
.
Zanjanian
H.
H.
Niksokhan
M. H.
Sarang
A.
2018
Influential third party on water right conflict: a Game Theory approach to achieve the desired equilibrium (case study: Ilam dam, Iran)
.
J. Environ. Manage.
214
,
283
294
.