Coupling Two-Stage Stochastic Robust Programming with Improved Export Coefﬁcient for Water Allocation among Industrial Sectors

: Water scarcity and water pollution are essential factors limiting coordinated regional development, especially in water-deprived regions. Industrial restructuring is an effective water management solution to alleviate water scarcity and mitigate water pollution. However, due to widely existing inexact parameter information in the water resource management system, it is challenging to allocate water resources among industrial sectors. To address these problems, an export coefﬁcient coupled with a two-stage stochastic robust programming method (EC-TSRP) was developed through integrating an export coefﬁcient model (ECM), two-stage stochastic programming (TSP) and robust optimization. The proposed EC-TSRP model could effectively deal with the multiple uncertainties expressed as stochastic and the intervals with ﬂuctuation ranges, and enhance the robustness of optimal plans for supporting water resource allocation among industrial sectors under complex uncertainties. It was then applied to Bayan Nur City, in arid north-west China. The optimization alternatives indicate that wheat, sheep and services would be the most sensitive sectors among all industrial sectors, when non-point source (NPS) pollution exports are restricted. In addition, comparing the EC-TSRP results with the deterministic model, the reliability of the system could be improved signiﬁcantly, while the value of the objective function would be decreased slightly. The simulation results were also compared with the historical data from 2012 to 2016. Although the total revenue of Bayan Nur City would decrease by 1.52%, the pollutant loads of total nitrogen, total phosphorus and chemical oxygen demand (TN, TP and COD) would decrease by 14.5%, 7.75% and 2.07%, respectively, and total water allocation also would decrease from 4.6 billion m 3 to 4.23 billion m 3 .


Introduction
Water is one of the most important resources for human survival and development [1,2]. As an essential input of productive activities, water undoubtably has a significant impact on economic development [3]. In recent years, owing to the increasing water demand of industrial sectors activities and human life [4], as well as the excessive water resource exploitation and increasing wastewater discharge, water scarcity is getting increasingly intense worldwide [5]. In many regions, water scarcity exacerbates conflicts between water supply and demand for the advance of urbanization and industrialization, which has become a critical constraint for watershed socioeconomic development [6,7]. problem. However, robust optimization methods were hardly coupled with TSP to address the optimization water resource allocation among industrial sector problems for tackling such uncertain parameters without detailed distributions and unknown upper and lower bounds.
In addition, optimization of water use among industrial sectors involves many aspects, such as society, the economy and the environment, and is a complex system [36,37]. Particularly, the environmental effects brought about by industrial restructuring cannot be ignored [38]. Due to the influence of industry characteristics and water consumption differences, there are obvious differences in the environmental effects of water resources in different industrial sectors [39]. It is necessary to efficiently determine the discharge amounts of pollution loads in the management process of water pollution from industrial sectors. An export coefficient model (ECM) is a promising means to effectively determine the discharge of pollution loads [40]. It can be easily embedded into an optimization framework due to its high computational efficiency [41,42]. In this study, ECM was adopted to quantify non-point source pollution (NPS), but due to the lack of pollutant data, inaccuracy and the uncertain characteristics of NPS exports, the export coefficients were usually expressed as uncertain parameters. The export coefficients were expressed as uncertain intervals with fluctuation range parameters.
Therefore, this paper aims at developing an export coefficient coupled with a twostage stochastic robust programming method (EC-TSRP). It could not only deal with the multiple uncertainties expressed as stochastic and interval within the fluctuation range, but could also enhance the robustness of optimal plans for supporting water resource allocation among industrial sectors under complex uncertainties. The applicability of the proposed method has been demonstrated through a case study in arid and semi-arid areas of China, Bayan Nur City. The optimization results could provide decision-makers with an optimal scheme for water allocation among industrial sectors under different scenarios. This paper is organized as follows: the development of the TSRP method is explained in Section 2; the developed EC-TSRP method is applied to a real-world case in north-west China in Bayan Nur City in Section 3; the results and discussions are interpreted and compared in Section 4; and Section 5 provides the concluding remarks of this research.

Two-Stage Stochastic Robust Programming Model
TSP is effective at addressing problems where an analysis of policy scenarios is desired and the uncertain parameters are described by probability distribution functions [21,23]. However, in the real world, there usually exist multiple uncertainties in the management system. Serval TSP-based studies have been developed to tackle interval numbers and/or fuzzy sets with unknown possibilities and/or probabilistic distributions. However, there are often many uncertain parameters without detailed distributions and unknown upper and lower bounds, which greatly affect management system performances. These studies could not deal with such uncertain parameters. In order to overcome this problem, a TSRP method is proposed in this study. In TSRP, the objective function is based on a TSP model, and the robust optimization method is introduced to the constraints to deal with such uncertain parameters without detailed distribution information.
A general TSRP model can be formulated as follows: The TSRP model solution is mainly divided into two steps. Firstly, a protection function [34] is introduced into the model constraints to realize the degree of protection of uncertain parameters in the constraints and the regulation of the system risk level. Secondly, a nonlinear protection function is transformed into linearity. Through the above two steps, the complex TSRP model is transformed into a simple linear optimization model.
Firstly, a protection function is introduced into the model. Protection function, β i ((X − Y) * , τ i ), is introduced into the model constraints, which is equivalent to adding a "buffer" to the process of the model to reduce the impact of uncertainty of the parameters involved in the constraint and to make the results more reliable. Where Y is the decision variable, and (X − Y) * is a specific set of solutions. τ i is the protection level of constraint i, which refers to the number of intervals with a fluctuation range in constraint i and that are allowed to fluctuate at their bounds under the premise of ensuring the feasibility of constraint i. From Formula (1), constraint ↔ A(X − Y) ≤ ↔ B is the sum of n values, and at most n uncertain parameters may vary within the intervals in the fluctuation range, and then at the protection level, τ ∈ [0, n]. Protection function can be protection τ parameters at most: the obtained model results are feasible when most τ parameters take arbitrary values within the intervals within the fluctuation range and other parameters are equal to the nominal values A; if there are more than τ parameters at the intervals with the fluctuation range of arbitrary values, the model still has a high probability. Model (1) can then be converted into the following problem: Secondly, the protection function is transformed into linearity. The protection function is formulated as follows: where τ i is the protection level of constraint i,τ i ∈ [0, card(J i )]; J i is the set of intervals with the fluctuation range in constraint i; card represents the number of elements in the set; S i is a subset of J i ; τ i denotes the integral part of τ i . When constraint i is tantamount a deterministic one without any protection against infeasibility. The setting of protection level can help regulate the robustness of the model. Robustness means that the model results still have a high reliability when the model is affected by uncertain factors [29,43]. By changing the protection level, comprehensive benefits and reliability of the model can be regulated and weighed.
Water 2022, 14,1947 5 of 26 The model transformation is as follows: where c j , d j , x j , y j , a ij , b i , a i j , b i , are the elements of their corresponding matrices (e.g., a ij denotes the elements of the ith row and jth column of the A matrix); and z i , and p i,j are new variables.

An Improved Export Coefficient Model
In the real-world water allocation problems, it is critical to determine the total discharge amounts of TN and TP from each zone, which vary significantly with the precipitation amount, geomorphology, land use and soil types of the watershed. In this study, the improved ECM was introduced to the optimization framework in order to estimate the total discharge amounts of nutrients. The TN and TP can be calculated as follows: where ↔ L is the annual loss of nutrients (kg); ↔ E i is the export coefficient for nutrient source i (kg/ca/year or kg/km 2 /year); A i is the area of the catchment occupied by land-use type i (km 2 ), or the number of livestock and poultry type i (ca), or the number of population (ca); I i is the nutrient input to source i (kg); and P is the nutrient input from precipitation (kg). Because the probability and possibility of the distribution information of export coefficients (i.e., ↔ E i ) were hardly obtained, intervals within the fluctuation range were selected to deal with the uncertainties in these parameters. Figure 1 illustrates the framework of the EC-TSRP method.

Study Problem
The study area shown in Figure 2 is located in arid and semi-arid areas (40°13′-42°28′ N, 105°12′-109°53′ E), Bayan Nur City, Inner Mongolia, northwest China. It consists of Linhe District (LH), Wuyuan County (WY), Dengkou County (DK), Urad Front Banner (UF), Urad Middle Banner (UM), Urad Rear Banner (UR) and Hanghou Banner (HH), and covers an area of approximately 6.4 × 10 4 km 2 . The cultivated land area is 7.16 × 10 5 hm 2 , and the per capita is 0.35 hm 2 , which is 3.5 times the national average level. The average annual precipitation of the region is only 130~285 mm and is concentrated from July to August. However, the average annual evaporation of this region reaches up to 2030~3180 mm. The average annual volume of water diverted from the Yellow River is about 4 billion m 3 , accounting for 1/7 of the transit water.
As the largest consumer of water in Bayan Nur City, the amount of water consumed by agriculture does not have a positive correlation with the profits generated, and water resources are greatly wasted. According to the statistical yearbook of Bayan Nur City in 2016, for example, the annual water consumption proportion of agriculture was up to 89.55%, but its economic benefit only accounted for 17.28% of the total benefits, while industry contributed to 42.4% of the economic benefits, but only 0.19% water was consumed. Irrational water use among industrial sectors leads to low efficiency of water utilization and seriously hinders local economic development. At the same time, the water quality in this region has deteriorated dramatically due to intensive industrial activities. The return water from farmland containing a large number of pollutants and untreated wastewater discharged by other industries enter the Wuliangsuhai lake, which aggravate the eutrophication of the lake [44].

Study Problem
The study area shown in Figure 2 is located in arid and semi-arid areas ( , and covers an area of approximately 6.4 × 10 4 km 2 . The cultivated land area is 7.16 × 10 5 hm 2 , and the per capita is 0.35 hm 2 , which is 3.5 times the national average level. The average annual precipitation of the region is only 130~285 mm and is concentrated from July to August. However, the average annual evaporation of this region reaches up to 2030~3180 mm. The average annual volume of water diverted from the Yellow River is about 4 billion m 3 , accounting for 1/7 of the transit water.
As the largest consumer of water in Bayan Nur City, the amount of water consumed by agriculture does not have a positive correlation with the profits generated, and water resources are greatly wasted. According to the statistical yearbook of Bayan Nur City in 2016, for example, the annual water consumption proportion of agriculture was up to 89.55%, but its economic benefit only accounted for 17.28% of the total benefits, while industry contributed to 42.4% of the economic benefits, but only 0.19% water was consumed. Irrational water use among industrial sectors leads to low efficiency of water utilization and seriously hinders local economic development. At the same time, the water quality in this region has deteriorated dramatically due to intensive industrial activities. The return water from farmland containing a large number of pollutants and untreated wastewater discharged by other industries enter the Wuliangsuhai lake, which aggravate the eutrophication of the lake [44].

Data Collection and Treatment
Data collection involved in the objective function is as follows. Based on the Bayan Nur statistical yearbook from 1990 to 2017, China Agricultural Information Network, industrial water use quota in Inner Mongolia, relevant references [45], crop yield, crop price, crop irrigation quota, utilization coefficient of water for crop irrigation and the scale of different livestock species were obtained. Based on the Bayan Nur Water Resources Bulletin and the government planning report-the quota of added value of ¥10 4 in the secondary industry and the tertiary industry-the total available water resources of each industry and the total allowable discharge of pollutants were obtained. According to the first National survey of pollution sources, pollutant discharge coefficients of each industry were obtained. Table 1 shows the basic economic parameters of each sector. According to the available water resources of Bayan Nur City from 1992 to 2018, the available water resources were classified as high, medium and low levels by the frequency analysis method, which were 48.1, 49.5 and 5.08 billion m³, respectively. The available water for various industrial sectors under different probability distributions are shown in Table 2. By collecting the export coefficients of relevant and similar areas in China, and by referring to the various researchers, the export coefficients of different crop types and livestock species were further determined according to the field situation and the relevant experimental studies [46][47][48][49][50][51][52][53]; the export coefficients are shown in Table 3. The minimum water requirement of each industry was set by referring to the annual field planning and industrial output values in the study area.
In order to facilitate the calculation of this study, the following assumptions were set up: (1) The probability of low, medium and high flow levels in the forecast year are 20%, 60% and 20% respectively. (2) During the planning period, it is assumed that water consumption, unit output value, cost and other parameters of each industry rise steadily over time, while the target water allocation decreases gradually due to the improvement of the water utilization coefficient. (3) Assuming that all industries are short of water, the target water allocation is the existing water allocation multiplied by a coefficient, which is taken as 1.3 in this study. (4) It is assumed that the saved water due to the reduction of

Data Collection and Treatment
Data collection involved in the objective function is as follows. Based on the Bayan Nur statistical yearbook from 1990 to 2017, China Agricultural Information Network, industrial water use quota in Inner Mongolia, relevant references [45], crop yield, crop price, crop irrigation quota, utilization coefficient of water for crop irrigation and the scale of different livestock species were obtained. Based on the Bayan Nur Water Resources Bulletin and the government planning report-the quota of added value of ¥10 4 in the secondary industry and the tertiary industry-the total available water resources of each industry and the total allowable discharge of pollutants were obtained. According to the first National survey of pollution sources, pollutant discharge coefficients of each industry were obtained. Table 1 shows the basic economic parameters of each sector. According to the available water resources of Bayan Nur City from 1992 to 2018, the available water resources were classified as high, medium and low levels by the frequency analysis method, which were 48.1, 49.5 and 5.08 billion m 3 , respectively. The available water for various industrial sectors under different probability distributions are shown in Table 2. By collecting the export coefficients of relevant and similar areas in China, and by referring to the various researchers, the export coefficients of different crop types and livestock species were further determined according to the field situation and the relevant experimental studies [46][47][48][49][50][51][52][53]; the export coefficients are shown in Table 3. The minimum water requirement of each industry was set by referring to the annual field planning and industrial output values in the study area.
In order to facilitate the calculation of this study, the following assumptions were set up: (1) The probability of low, medium and high flow levels in the forecast year are 20%, 60% and 20% respectively. (2) During the planning period, it is assumed that water consumption, unit output value, cost and other parameters of each industry rise steadily over time, while the target water allocation decreases gradually due to the improvement of the water utilization coefficient. (3) Assuming that all industries are short of water, the target water allocation is the existing water allocation multiplied by a coefficient, which is taken as 1.3 in this study. (4) It is assumed that the saved water due to the reduction of agricultural target water allocation will be distributed to other industries. (5) During the planning period, it is assumed that each export coefficient decreases with time.

Model Formulation
Aiming at the problems of unbalanced water allocation among industrial sectors and prominent contradictions between water supply and demand in the study area, an EC-TSRP model suitable for this area was developed. The objective of this model is to maximize the total economic output value of the whole region of Bayan Nur City, which include those from primary, secondary and tertiary industries. The primary industry considers agriculture and livestock husbandry (the economic output value of agriculture and livestock husbandry accounts for about 95% of the primary industry). The secondary industry considers industry and construction. The tertiary industry is considered as the services sector. The decision variables are the water allocation in the first stage and the water deficit in the second stage for each sector in each subarea under each water flow level. According to the principle of the EC-TSRP model, the first stage of the model is to allocate water (i.e., the water target) to different industrial sectors, directly leading to various water resource availabilities; the second stage will be composed of the sectors involved in the calculation of the loss of benefits resulting from water shortage. The planning period is divided into three periods from 2012 to 2026, and each period is set at five years (2012-2016 is th historical year P1, 2017-2021 is the current year P2, and 2022-2026 is the planning year P3). The formulation of the model is presented as follows: Objection function: where F is the objective function value, which represents the total economic output value of Bayan Nur City (¥); U is the index for subareas (u = 1 for LH District, u = 2 for WY County, u = 3 for DK County, u = 4 for UF Banner, u = 5 for UM Banner, u = 6 for UR Banner, u = 7 for HH Banner); I is the index for planning periods (i = 1 for 2012-2016, i = 2 for 2017-2021, i = 3 for 2022-2026); J is the index for sectors of all industries (primary industry: agriculture, j = 1 for wheat, j = 2 for corn, j = 3 for sunflower, j = 4 for melon, j = 5 for tomato, livestock husbandry, j = 6 for cattle, j = 7 for sheep, j = 8 for pig, secondary industry: j = 9 for industry, j = 10 for construction, tertiary industry: j = 11 for services); H is index for flow level (h = 1 for low flow level, h = 2 for medium flow level, h = 3 for high flow level); NY i is the length of period i (year). β i is the discount factor of period , where δ and θ are the annual interest and inflation rates, respectively, and t represents the t year of the planning period. Therefore, the discount factor for a period is the average of the years involved.) NB uij , NT uij and NC uij are per unit water supply benefits (¥/m 3 ), target water allocation (m 3 ) and unit water shortage penalties (¥/m 3 ) of sector j for subarea u in period i, respectively NS uijh is the water shortage amount of sector j for subarea u in period i in flow level h (m 3 ). p h is the occurrence probability of water flow level. Subject to: (1) Production constraints: Primary industry production constraints: Water 2022, 14,1947 10 of 26 where NZW ui is the annual average economic output value of the primary industry in subarea u in period i (¥10 4 /year); FGL ui is the minimum increase rate of primary industry output value in subarea u from period i to i + 1(%); FGH ui is the maximum increase rate of primary industry output value in subarea u from period i to i + 1(%). Secondary industry production constraints: Industrial production constraints: where SZW ui is the annual average industrial economic output value in subarea u in period i (¥10 4 /year); IGL ui is the minimum increase rate of industrial output value in subarea u from period i to i + 1(%); IGH ui is the maximum increase rate of industrial output value in subarea u from period i to i + 1(%). Construction production constraints: where JZW ui is the annual average construction economic output value in subarea u in period i (¥10 4 /year); JGL ui is the minimum increase rate of construction output value in subarea u in period i (¥10 4 /year) (%); JGH ui is the maximum increase rate of construction output value in subarea u in period i (¥10 4 /year) (%). Tertiary industry production constraints: where TZW ui is the annual average economic output value of the tertiary industry in subarea u in period i (¥10 4 /year); TGL ui is the minimum increase rate of tertiary industry output value in subarea u in period i (¥10 4 /year) (%); TGH ui is the maximum increase rate of tertiary industry output value in subarea u in period i (¥10 4 /year) (%).
(2) Resources constraints: Water supply-demand balance constraints: Water 2022, 14,1947 11 of 26 where Q ih is the total available water resources of period i in flow level h (m 3 ); MAXWATN uih , MAXWATX uih , MAXWATI uih and MAXWATN uih are the total available water resources of agriculture, livestock husbandry, secondary industry and tertiary industry under the inflow level h in subarea u in period i (m 3 ). Sectoral minimum water demand constraints: where IW uij is the minimum water requirement of sector j in subarea u in period i (m 3 ).
(3) Environmental constraints: Dissolved nitrogen loss constraints: where ← − → ECN ij is the export coefficient of TN for sector j in period i (kg/hm 2 /year or kg/ca/year); ← −− → DISN ih is the maximum allowable dissolved nitrogen loss of primary industry under the inflow level h in period i (kg).
Dissolved phosphorus loss constraints: where ← − → ECP ij is the export coefficient of TP for sector j in period i (kg/hm 2 /year or kg/ca/year); ← −− → DISP ih is the maximum allowable dissolved phosphorus loss of primary industry under the inflow level h in period i (kg).
COD discharge constraints: where ← → LC ij is the COD discharge from secondary and tertiary industries of sector j in period i (kg/¥10 4 /year); ←−−−−−→ MAXCOD ih is the maximum allowable COD discharges under the inflow level h in period i (kg).

Results and Discussion
The optimization of water use among industrial sectors at Bayan Nur City were obtained from solving the above EC-TSRP model. For export coefficients and ← → LC ij , it was difficult to obtain deterministic values or precise distributions, which were expressed as intervals with fluctuation ranges. Take ← − → ECN ij , for example: it could be expressed as a fluctuation range between ECN ij − EĈN ij and ECN ij + EĈN ij , defined by a nominal value ECN ij and its fluctuation radius EĈN ij (Figure 3). Meanwhile, in order to quantify the influence of the uncertainty degree of the export coefficients on the results for data deficiency, three fluctuation levels-B30, B20 and B10-were selected; namely, the fluctuation radius was 30%, 20% and 10% of the export coefficient value. Considering the acceptance degree of local decision-makers to the target risk, three protection levels-high, medium and low-were selected and defined as L10, L5 and L1, corresponding to the protection of no more than 10%, 5% and 1% of the constraint weight, respectively.
The scenario settings are shown in Table 4. By comparing different scenarios, the effects of the uncertainty degree of export coefficients (fluctuation level) and the protection degree of environmental constraints (6h), (6i) and (6j) (protection level) on the results were observed and explored respectively. S1 was the scenario with the lowest parameter fluctuation level and constraint protection level; S9 was the scenario with the highest parameter fluctuation level and constraint protection level, and the results of the remaining scenarios were between these two scenarios. Therefore, for the convenience of analysis and discussion of some results, four scenarios-S1, S7, S8 and S9-were selected to illustrate this.

Results and Discussion
The optimization of water use among industrial sectors at Bayan Nur City were ob-  Figure 3). Meanwhile, in order to quantify the influence of the uncertainty degree of the export coefficients on the results for data deficiency, three fluctuation levels-B30, B20 and B10-were selected; namely, the fluctuation radius was 30%, 20% and 10% of the export coefficient value. Considering the acceptance degree of local decision-makers to the target risk, three protection levelshigh, medium and low-were selected and defined as L10, L5 and L1, corresponding to the protection of no more than 10%, 5% and 1% of the constraint weight, respectively. The scenario settings are shown in Table 4. By comparing different scenarios, the effects of the uncertainty degree of export coefficients (fluctuation level) and the protection degree of environmental constraints (6h), (6i) and (6j) (protection level) on the results were observed and explored respectively. S1 was the scenario with the lowest parameter fluctuation level and constraint protection level; S9 was the scenario with the highest parameter fluctuation level and constraint protection level, and the results of the remaining scenarios were between these two scenarios. Therefore, for the convenience of analysis and discussion of some results, four scenarios-S1, S7, S8 and S9-were selected to illustrate this.    Figure 4a shows the industry scale in three planning periods under 4 scenarios of S1-S7, S8 and S9. All of the industries would have an increase in the economic benefits over the planning period in the study area. Among these, the growth of the secondary industry would be the most obvious, especially in the S1 scenario, where the economic benefit would increase from 225.5 billion to 483.7 billion over periods 1-3. This is mainly because, when setting the target water allocation, the model assumes that the target water allocation for the secondary industry increased as the time lapses (reducing the agricultural target water allocation to other sectors), and that the water quota for the secondary industry of the ¥10 4 value added decreased every year. Figure 4b exhibits the contribution of each industry to the total revenue, the share of the secondary industry would gradually increase, the primary industry would continuously decrease and the tertiary industry would tend to be stable. This is also in line with the current situation in the study area, where the secondary industry dominates the share of economic output. Comparing scenarios S1 and S7, the total returns of the primary and tertiary industries would show a decreasing trend from periods 1-3, while the secondary industry would remain largely unchanged, indicating that the increase in data fluctuation level would cause the decline of the economic output value of each industry. Comparing scenarios S7, S8 and S9, the primary industry would show the different decreasing trends in the three periods, and the secondary and tertiary industries would show decreasing trends in period 2 and period 3, and remain basically unchanged in period 1. This indicates that the increase in the constraint protection level would also lead to the decline in economic output value of each industry, and there were differences in the sensitivity of each industry to protection for different planning periods.

Results and Sensitivity Analysis
industry would be the most obvious, especially in the S1 scenario, where the economic benefit would increase from 225.5 billion to 483.7 billion over periods 1-3. This is mainly because, when setting the target water allocation, the model assumes that the target water allocation for the secondary industry increased as the time lapses (reducing the agricultural target water allocation to other sectors), and that the water quota for the secondary industry of the ¥10 4 value added decreased every year. Figure 4b exhibits the contribution of each industry to the total revenue, the share of the secondary industry would gradually increase, the primary industry would continuously decrease and the tertiary industry would tend to be stable. This is also in line with the current situation in the study area, where the secondary industry dominates the share of economic output. Comparing scenarios S1 and S7, the total returns of the primary and tertiary industries would show a decreasing trend from periods 1-3, while the secondary industry would remain largely unchanged, indicating that the increase in data fluctuation level would cause the decline of the economic output value of each industry. Comparing scenarios S7, S8 and S9, the primary industry would show the different decreasing trends in the three periods, and the secondary and tertiary industries would show decreasing trends in period 2 and period 3, and remain basically unchanged in period 1. This indicates that the increase in the constraint protection level would also lead to the decline in economic output value of each industry, and there were differences in the sensitivity of each industry to protection for different planning periods.  The above has described the optimization results and the sensitivity analysis of the scale changes of the primary, secondary and tertiary industries of the EC-TSRP model under different scenarios. In the following section, we analyze each sector, specifically those relating to agriculture, livestock husbandry, industry, construction and services and The above has described the optimization results and the sensitivity analysis of the scale changes of the primary, secondary and tertiary industries of the EC-TSRP model under different scenarios. In the following section, we analyze each sector, specifically those relating to agriculture, livestock husbandry, industry, construction and services and explore the extent to which changes in the level of protection and fluctuation affect each sector, by observing the results of different scenarios for each sector in detail. Figure 5 illustrates the overall agricultural water allocation in Bayan Nur City under the four scenarios of S1, S7, S8 and S9. The target water allocation was set according to policy guidelines, and it decreased gradually over the three planning periods. Based on the set target water allocation, the following model results were obtained. Taking the S1 scenario as an example, the results show that the water allocation for agriculture would decrease over the planning horizon, gradually reducing the water consumption and controlling it at 4 billion m 3 . Among these, the largest proportion of water deficiency in wheat would be 31.89%, followed by sunflower, corn, tomato and melon, with a water deficiency of 16.13%, 9.57%, 0.49% and 0.12%, respectively. Agricultural water allocation first meets cash crops and then crops, because crops contribute to lower economic benefits but consume higher amounts of water and produce more pollutants. For example, wheat consumes 13.81% of the water, but only brings 5.58% of the economic benefits and produces 13.87% of the total nitrogen (TN) and total phosphorus (TP) pollutants, as shown in Figure 6a. of 16.13%, 9.57%, 0.49% and 0.12%, respectively. Agricultural water allocation first meets cash crops and then crops, because crops contribute to lower economic benefits but consume higher amounts of water and produce more pollutants. For example, wheat consumes 13.81% of the water, but only brings 5.58% of the economic benefits and produces 13.87% of the total nitrogen (TN) and total phosphorus (TP) pollutants, as shown in Figure  6a. Comparing scenarios S1 and S7, the water allocation of wheat and corn in the three planning periods would decrease gradually, and sunflowers would show an increasing trend, while melons and tomatoes would remain unchanged. The overall agricultural water allocation in Bayan Nur City would decrease from 4.06 to 3.99 billion m³. This indicates that when the data fluctuation level increases, the overall agricultural water allocation would decrease within the agricultural subsystem, and it would be shifted towards those crop types with low water consumption and pollution exports. Because melons and Comparing scenarios S1 and S7, the water allocation of wheat and corn in the three planning periods would decrease gradually, and sunflowers would show an increasing trend, while melons and tomatoes would remain unchanged. The overall agricultural water allocation in Bayan Nur City would decrease from 4.06 to 3.99 billion m 3 . This indicates that when the data fluctuation level increases, the overall agricultural water allocation would decrease within the agricultural subsystem, and it would be shifted towards those crop types with low water consumption and pollution exports. Because melons and tomatoes are not short of water, some of the water allotment for wheat and corn is deployed to sunflower. Comparing scenarios S7, S8 and S9, the water allocation of wheat would decrease in the three planning periods, and corn and sunflower would show a downward trend in periods 1 and 2, while melons and tomatoes would remain unchanged. The overall agricultural water allocation in Bayan Nur City would decrease from 3.99 to 3.74 billion m 3 . When the protection level increases, the overall water allocation would decrease. The agricultural subsystem would first allocate the water in periods 1 and 2, and would then allocate the water in period 3. This is because the model assumes increasing water benefits per crop as time progresses. Figure 7a shows the export of agricultural TN and TP pollution in Bayan Nur City under nine scenarios. As the level of protection and fluctuation increases, pollution load of TN and TP would decrease. According to the analysis in Figure 7b, in which the pollutants produced by wheat would keep decreasing, the pollutants produced by corn and sunflower would show an overall decreasing trend with a local increase, while melons and tomatoes would remain almost unchanged. This is consistent with the conclusions drawn from the changes in crop water allocation for the different scenarios shown in Figure 5. under nine scenarios. As the level of protection and fluctuation increases, pollution load of TN and TP would decrease. According to the analysis in Figure 7b, in which the pollutants produced by wheat would keep decreasing, the pollutants produced by corn and sunflower would show an overall decreasing trend with a local increase, while melons and tomatoes would remain almost unchanged. This is consistent with the conclusions drawn from the changes in crop water allocation for the different scenarios shown in  Figure 8 illustrates the overall livestock husbandry water allocation in Bayan Nur City under the four scenarios of S1, S7, S8 and S9. In the S1 scenario, for example, the proportion of water deficiency would be 19.6%, 36.7% and 0 for sheep, pigs and cattle, respectively. Pigs would be the most water-scarce sector, due to the fact that they have the under nine scenarios. As the level of protection and fluctuation increases, pollution load of TN and TP would decrease. According to the analysis in Figure 7b, in which the pollutants produced by wheat would keep decreasing, the pollutants produced by corn and sunflower would show an overall decreasing trend with a local increase, while melons and tomatoes would remain almost unchanged. This is consistent with the conclusions drawn from the changes in crop water allocation for the different scenarios shown in  Figure 8 illustrates the overall livestock husbandry water allocation in Bayan Nur City under the four scenarios of S1, S7, S8 and S9. In the S1 scenario, for example, the proportion of water deficiency would be 19.6%, 36.7% and 0 for sheep, pigs and cattle, respectively. Pigs would be the most water-scarce sector, due to the fact that they have the  Figure 8 illustrates the overall livestock husbandry water allocation in Bayan Nur City under the four scenarios of S1, S7, S8 and S9. In the S1 scenario, for example, the proportion of water deficiency would be 19.6%, 36.7% and 0 for sheep, pigs and cattle, respectively. Pigs would be the most water-scarce sector, due to the fact that they have the lowest economic efficiency, consuming 14.83% of the water but contributing only 9.79% of the economy, as shown in Figure 7b. Comparing scenarios S1 and S7, the water allocation for sheep would decrease somewhat in all three planning periods, while the water allocation for pigs would show an opposite upward trend, and the overall water allocation for livestock husbandry would decrease from 0.34 to 0.31 billion m 3 . This indicates that when the data fluctuation level increases, the overall water allocation of livestock husbandry would decrease, and the water allocation within the subsystem would shift towards livestock species with low pollution exports. For example, pigs produce only 5.12% TN pollutants and 9.89% TP pollutants, while consuming 14.83% water. Comparing scenarios S7, S8 and S9, the water allocation of cattle, sheep and pigs would all decrease in these three planning periods. Figure 9 shows the export of TN and TP pollution loads from livestock husbandry in Bayan Nur City under nine scenarios, which is the same as the pollution load export performance of agriculture shown in Figure 6. The pollution exports of TN and TP would decrease as the level of protection and fluctuation increases. Among them, the pollution produced by sheep would decrease continuously, while the cattle and pigs would have an overall decreasing trend, with a local increase. these three planning periods. Figure 9 shows the export of TN and TP pollution loads from livestock husbandry in Bayan Nur City under nine scenarios, which is the same as the pollution load export performance of agriculture shown in Figure 6. The pollution exports of TN and TP would decrease as the level of protection and fluctuation increases. Among them, the pollution produced by sheep would decrease continuously, while the cattle and pigs would have an overall decreasing trend, with a local increase.  these three planning periods. Figure 9 shows the export of TN and TP pollution loads from livestock husbandry in Bayan Nur City under nine scenarios, which is the same as the pollution load export performance of agriculture shown in Figure 6. The pollution exports of TN and TP would decrease as the level of protection and fluctuation increases. Among them, the pollution produced by sheep would decrease continuously, while the cattle and pigs would have an overall decreasing trend, with a local increase.   Figure 10 illustrates the water allocation of industry, construction and services in Bayan Nur City under four scenarios: S1, S7, S8 and S9. Similarly, taking S1 as an example, the water allocation of industry, construction and service sectors would continue to grow in three planning periods. Especially at the high flow level, the industrial water allocation would increase by 24.5% from period 2 to period 3. Comparing scenarios S1 and S7, the water allocation of industry and construction would remain unchanged, and the services sector would decrease. This indicates that industry and construction are not sensitive to changes in data fluctuation levels, which is consistent with the conclusion in Figure 5. Comparing scenarios S7, S8 and S9, the water allocation for construction and services would decrease and that for industry would increase; from S8 to S9, the water allocation for construction and services would remain unchanged and that for industry would decrease. This indicates that as the constraint protection level increases, the total water allocation would decrease, and the system would prioritize the allocation of sectors with high COD pollution exports per unit of output value, and would then allocate other sectors.
As shown in Figure 11, the COD pollution export would decrease continuously as the level of protection and fluctuation increases. The services sector would have the most significant decrease in COD pollution export. This is mainly because the COD export per unit of output value is higher in the service sector compared with industry and construction, and the system would prioritize limiting the water distribution in the service sector to reduce the pollutant export. Figure 11, the COD pollution export would decrease continuously as the level of protection and fluctuation increases. The services sector would have the most significant decrease in COD pollution export. This is mainly because the COD export per unit of output value is higher in the service sector compared with industry and construction, and the system would prioritize limiting the water distribution in the service sector to reduce the pollutant export.

As shown in
(a) S1 and S7 (b) S8 and S9 Figure 10. The overall water allocation of industry, construction and services in Bayan Nur City under four scenarios: S1, S7, S8 and S9 (Ind: industry, Con: construction, Ser: services).
As shown in Figure 11, the COD pollution export would decrease continuously as the level of protection and fluctuation increases. The services sector would have the most significant decrease in COD pollution export. This is mainly because the COD export per unit of output value is higher in the service sector compared with industry and construction, and the system would prioritize limiting the water distribution in the service sector to reduce the pollutant export.
(a) S1 and S7 (b) S8 and S9 Figure 10. The overall water allocation of industry, construction and services in Bayan Nur City under four scenarios: S1, S7, S8 and S9 (Ind: industry, Con: construction, Ser: services). The results and sensitivity analyses for Bayan Nur City overall agriculture, livestock husbandry, industry, construction and services under different scenarios are described above. Each district and county can also be analyzed in a similar way. Take S1 as an example for a brief introduction, as shown in Figure 12.
In agricultural water allocation, wheat would be the most water-scarce crop, facing different degrees of water scarcity in almost all districts and counties, especially in WY, where the proportion of water scarcity would reach 59%. It is followed by corn, which would show different degrees of water scarcity in DK, UF, UM and UR; sunflower would exhibit water scarcity in UF and UM; melons and tomatoes would belong to the least water-scarce crops in planting, indicating that from the point of view of economic benefits and water saving requirements, cash crops (melons and tomatoes) would be a better choice. In the livestock husbandry water allocation, the dominant sector in each district and county is sheep, which would consume about 80% of water for the livestock industry, but would also be the sector with the greatest water shortage, especially in LH and HH, where the proportion of water shortage for sheep would reach 34% and 28%. In the industrial water allocation, there would be differences in the growth trends of each district and county in three planning periods, mainly in two forms. LH, WY, DK and HH all would exhibit a steady upward trend in each period, while UF and UR would show a downward, then an upward trend.

Discussion
In order to assess the performance and merits of the model, this paper discusses three aspects: model comparison, status quo comparison, and PBV risk analysis.

Discussion
In order to assess the performance and merits of the model, this paper discusses three aspects: model comparison, status quo comparison, and PBV risk analysis.  [35]. For example, in pursuit of a balance between economic output value and pollution exports, the planning of decision makers to deal with the actual situation changes from year to year.
The impacts of fluctuation and protection level on the results: Figure 13a compares the changes of total economic output value in nine scenarios. The protection level corresponds to different environmental protection policies; higher levels of protection means more severe environmental management constraints. As the level of protection increases, the total economic output would decrease continuously, but there are differences in the degree of decline at distinct levels of fluctuation. For instance, under fluctuation level L1, when the protection level increases from low to high, the total economic output value would gradually descend from 1.75 trillion to 1.71 trillion. Under fluctuation level L10, these values would gradually descend from 1.73 trillion to 1.56 trillion. The fluctuation level corresponds to different degrees of inaccuracy of parameters, and higher levels of fluctuation mean greater levels of inaccuracy. As the level of fluctuation increases, the total economic output would also decrease continuously, and there are differences in the degree of decline at distinct protection levels. Specifically, under protection level B10, when the fluctuation level increases from low to high, the total economic output value would gradually descend from 1.75 trillion to 1.73 trillion. Under protection level B30, these values would gradually descend from 1.71 trillion to 1.56 trillion. This implies that the benefit and robustness of the system requires a trade-off between the fluctuation level of the data and the protection level of the constraints. A higher level of protection and fluctuation would increase the chance of system feasibility, but would lead to a lower system benefit. Figure 13b shows the scale changes of each industry under nine scenarios. It is found that the scale of the primary and tertiary industries would continue to decline from scenario S1 to S9. The secondary industry, on the other hand, just starts to change at scenario S7. This indicates that the primary and tertiary industries are more sensitive to changes in the system compared to the secondary industry.  Figure 13b shows the scale changes of each industry under nine scenarios. It is found that the scale of the primary and tertiary industries would continue to decline from scenario S1 to S9. The secondary industry, on the other hand, just starts to change at scenario S7. This indicates that the primary and tertiary industries are more sensitive to changes in the system compared to the secondary industry. Model comparison: the EC-TSP model with only random uncertainty was established to solve the same problem. Differing from the EC-TSRP model, the uncertainty of the export coefficients was ignored. Assuming that all parameters would take nominal values, only a simplified set of deterministic values were used as input in the EC-TSP model. In fact, the EC-TSP model is only a special case among the many cases considered by the EC-TSRP model. Unlike the EC-TSRP model, which would protect the model from suboptimality or infeasibility related to data perturbations, the EC-TSP solution would be far from reliable when the value of any input parameter is not a determined value as given by the real problem. Figure 14a illustrates the variations of different indicators when comparing the EC-TSP results against the EC-TSRP solutions. As shown in scenario S7, the total economic benefit of the EC-TSRP model would decrease by 3.19% relative to the EC-TSP model. However, the total water consumption, TN load, TP load and COD discharges all would be reduced by 1.74%, 6.18%, 5.70% and 5.62% accordingly, and EC-TSRP could signifi- Model comparison: the EC-TSP model with only random uncertainty was established to solve the same problem. Differing from the EC-TSRP model, the uncertainty of the export coefficients was ignored. Assuming that all parameters would take nominal values, only a simplified set of deterministic values were used as input in the EC-TSP model. In fact, the EC-TSP model is only a special case among the many cases considered by the EC-TSRP model. Unlike the EC-TSRP model, which would protect the model from suboptimality or infeasibility related to data perturbations, the EC-TSP solution would be far from reliable when the value of any input parameter is not a determined value as given by the real problem. Figure 14a illustrates the variations of different indicators when comparing the EC-TSP results against the EC-TSRP solutions. As shown in scenario S7, the total economic benefit of the EC-TSRP model would decrease by 3.19% relative to the EC-TSP model. However, the total water consumption, TN load, TP load and COD discharges all would be reduced by 1.74%, 6.18%, 5.70% and 5.62% accordingly, and EC-TSRP could significantly reduce the system-violation risk by 37.43% due to the consideration of complex uncertainties. Scenarios S8 and S9 also show the same dynamics. The economic benefits would decrease by 7.31% and 12.12% respectively, and the system-violation risk could decrease by 46.17% and 47.26%. The three scenarios would differ in terms of the rate of improvement of system reliability, while the economic efficiency and pollution exports would decrease. Scenario S7, for example, would significantly improve the reliability of the system when the value of the objective function decreases slightly, which is more acceptable to decision makers. While S8 and S9 would improve the higher reliability and the decrease of the system benefit would be made more obvious, the ability to limit the pollutant export would be greatly improved, which is more suitable in the case of strict environmental policies. Therefore, in the actual situation, the degree of protection for the system should be weighed according to the wishes of decision makers. This shows that the EC-TSRP model can weigh the system benefits, pollution export and system reliability according to the realistic situation. Status quo comparison: Figure 15a shows the comparison of economic output value between the actual situation and planning results of each district and county in 2012-2016 under scenario S1, and Figure 15b shows the comparison of actual situation and planning results in different aspects of Bayan Nur City. The total economic output value of Bayan Nur City would decrease from 84.83 billion to 83.69 billion, a decrease of 1.34%. In agriculture and construction, the Bayan Nur City output value would show a decline, down 10.04%, 15.46% respectively; while the economic output value of livestock husbandry would increase by 30.66%; industrial and service economic output value would remain basically unchanged. The main reasons for the decline in economic benefits of the model were as follows: (1) Compared to the actual situation in 2012-2016, the pollutant emission values under S1 scenario would decrease by 14.5%, 7.75% and 2.07%, respectively, due to the model environmental constraints (constraints (6q), (6r) and (6s)), and the total water consumption in the study area would decrease from about 4.6 to 4.234 billion m³. Among them, the decline in water allocation would be most pronounced in agriculture, from 4459 to 4057 million m³. The secondary and tertiary industries meet this environmental policy by limiting water use in the construction and service sectors. If the saved water resources are used for ecological protection or other environmental protection sectors, it would not only compensate for the reduced economic benefits due to the reduced water allocation, but it also would promote the sustainable development of the regional ecology. (2) When agriculture selects crops, the five crop species with the largest share are selected, so the economic output is bound to be lower than the actual economic output of agriculture. In the actual water allocation, the EC-TSRP model could make a reasonable water allocation spe- Status quo comparison: Figure 15a shows the comparison of economic output value between the actual situation and planning results of each district and county in 2012-2016 under scenario S1, and Figure 15b shows the comparison of actual situation and planning results in different aspects of Bayan Nur City. The total economic output value of Bayan Nur City would decrease from 84.83 billion to 83.69 billion, a decrease of 1.34%. In agriculture and construction, the Bayan Nur City output value would show a decline, down 10.04%, 15.46% respectively; while the economic output value of livestock husbandry would increase by 30.66%; industrial and service economic output value would remain basically unchanged. The main reasons for the decline in economic benefits of the model were as follows: (1) Compared to the actual situation in 2012-2016, the pollutant emission values under S1 scenario would decrease by 14.5%, 7.75% and 2.07%, respectively, due to the model environmental constraints (constraints (6q), (6r) and (6s)), and the total water consumption in the study area would decrease from about 4.6 to 4.234 billion m 3 . Among them, the decline in water allocation would be most pronounced in agriculture, from 4459 to 4057 million m 3 . The secondary and tertiary industries meet this environmental policy by limiting water use in the construction and service sectors. If the saved water resources are used for ecological protection or other environmental protection sectors, it would not only compensate for the reduced economic benefits due to the reduced water allocation, but it also would promote the sustainable development of the regional ecology. (2) When agriculture selects crops, the five crop species with the largest share are selected, so the economic output is bound to be lower than the actual economic output of agriculture. In the actual water allocation, the EC-TSRP model could make a reasonable water allocation specifically based on the overall planning of environmental policies by decision makers.
PBV risk analysis: The PBV risk reflects the constraint (6h), (6i) and (6j) violation of the degree of protection. Figure 14b shows the PBV quantification of the risk of constraint violation for different degrees of protection. It could be seen that the system-violation risk would gradually decrease as the level of protection increases. This paper was relatively conservative, and three protection levels of 1%, 5% and 10% were chosen for calculation. The purpose is to minimize the decline of economic benefits while increasing the system reliability. PBV risk analysis also has some significant implications for policy guidance.
If improving water quality is the way forward, higher pollutant discharge standards are required and higher levels of protection are necessary to deal with the effects of pollution export; if economic development will remain the focus of future development and ecological protection is of secondary importance, the degree of protection of the system can be adjusted downward and it is reasonable to face a greater risk of violations as a way to obtain higher economic benefits. Therefore, in practice, the EC-TSRP model can help decision makers to develop an appropriate water allocation plan based on the policy document, the decision maker's wishes and environmental risks. The developed EC-TSRP method is a widely applicable programming approach that can simultaneously tackle multiple uncertainties in the process of optimization of water use among industrial sectors. This study represents the development and the first application of EC-TSRP to the field of water resources management in Bayan Nur city, where water scarcity and low water use efficiency has hindered local development. The developed methodology can be also applicable in the address of many other management problems, such as waste management and air pollution control and various watersheds.
In detail, two tasks need to be accomplished when applying the developed model to support other management problems. Firstly, determining the value and uncertainty of parameters. Due to the development and management gaps in various areas, the same parameter values will show different uncertainty characteristics. Thus, determining the uncertainty type and value range of parameters is the basis for the successful application of the developed model. Secondly, tailoring the EC-TSRP's constraints and objectives to the site-specific characteristics of other management cases. This study aims to maximize the total economic output value under the impact of environmental, economic and resource constraints. For different fields or areas, the current objective function (Equation (6a)) and constraint conditions (Equations (6b)-(6k)) can be replaced with ones that are in line with the development requirements of a given region. This means that quantifying the regional problems will be a challenge for the practical application of the developed model. The developed EC-TSRP method is a widely applicable programming approach that can simultaneously tackle multiple uncertainties in the process of optimization of water use among industrial sectors. This study represents the development and the first application of EC-TSRP to the field of water resources management in Bayan Nur city, where water scarcity and low water use efficiency has hindered local development. The developed methodology can be also applicable in the address of many other management problems, such as waste management and air pollution control and various watersheds.
In detail, two tasks need to be accomplished when applying the developed model to support other management problems. Firstly, determining the value and uncertainty of parameters. Due to the development and management gaps in various areas, the same parameter values will show different uncertainty characteristics. Thus, determining the uncertainty type and value range of parameters is the basis for the successful application of the developed model. Secondly, tailoring the EC-TSRP's constraints and objectives to the site-specific characteristics of other management cases. This study aims to maximize the total economic output value under the impact of environmental, economic and resource constraints. For different fields or areas, the current objective function (Equation (6a)) and constraint conditions (Equations (6b)-(6k)) can be replaced with ones that are in line with the development requirements of a given region. This means that quantifying the regional problems will be a challenge for the practical application of the developed model.

Conclusions
According to the uncertainty and complexity of the process of optimization of water use among industrial sectors, an EC-TSRP method by means of an integrated export coefficient model (ECM), two-stage stochastic programming (TSP) and robust optimization within a programming framework has been introduced to this study and was developed. It could not only effectively deal with the multiple uncertainties expressed as stochastic and intervals within a fluctuation range, but could also enhance the robustness of optimal plans for supporting water resource allocation among industrial sectors under complex uncertainties. With adjustable levels of protection, the EC-TSRP approach could provide the optimal solution, as a compromise between the desired objective value and the potential system-violation risks. In general, the EC-TSRP model demonstrated the following advantages. (1) It improved upon two-stage stochastic programming by being able to address the uncertain information without the known distributions. (2) It reflected the impact of different degrees of fluctuation of interval values on the boundary to the target. (3) It could provide decision makers with risk-explicit results. (4) It could simplify the complexity of the pollutant transfer process and provide suggestions for policy makers to solve the contradiction between regional economic development and environmental coordination.
The proposed EC-TSRP model was applied to Bayan Nur city (where there is a water scarcity and where pollution is heavy) as the study area and an empirical study was conducted. The model results show that: (1) As the level of data fluctuation increases, the water distribution within the system would be deployed to sectors with low water consumption and a low pollution export; therefore, the export amount of TN, TP and COD pollution would decrease, and the economic benefit would reduce. (2) As the level of constraint protection increases, the water allocation for each sector would decrease, and agricultural systems would prioritize the reduction of water allocation for planning periods 1 and 2. Livestock, industry, construction and services would prioritize lower water allocation in planning period 3; the export of TN, TP and COD pollution would decrease and the economic benefit also would decrease. However, the robustness of the system would gradually increase. (3) Wheat, sheep and service sectors would be the most sensitive sectors for pollutant emissions, with priority given to limiting water allocation when NPS exports are restricted.
The developed EC-TSRP model simulation results were also compared with the historical data from 2012 to 2016. Although the total revenue of Bayan Nur city would decrease by 1.52% compared with the actual situation, the pollutant loads of TN, TP and COD would decrease by 14.5%, 7.75% and 2.07%, respectively, and the total water allocation would also decrease from 4.6 billion m 3 to 4.234 billion m 3 . For Bayan Nur City's actual situation, agriculture should reduce certain crops, increase water consumption of cash crops to improve economic benefit; livestock husbandry should be consistent with the status quo, a sheep should remain the main object of development for livestock husbandry; secondary and tertiary industries should gradually increase water consumption in the future, first in a certain planning period to develop the secondary industry, and should then continue in order to vigorously develop the tertiary industry.
In this study, the empirical model (ECM) was embedded in the regional water-use structure optimization model of Bayan Nur city, but there was little research on the export coefficients of the region, and the accuracy of the export coefficients had a great influence on the sensitivity and rationality of the model. Therefore, the export coefficients need to be studied in greater depth, and this study can be followed up by experiments to seek the appropriate export coefficients for this study area and to increase the accuracy of the model.