Interactive Allocation of Water Pollutant Initial Emission Rights in a Basin under Total Amount Control: A Leader-Follower Hierarchical Decision Model

The initial emission rights allocation is the key measure to achieve the goal of total amount control and deepen the emission trading system. Although many studies have focused on the modeling of initial emission rights allocation, such as using game theory and multi-objective optimization methods, few studies have observed the hierarchical relationship of mutual interference and restriction between watershed management agency and local governments in each subarea during allocation. This relationship directly affects the rationality of the results of regional emission rights allocation. In this study, a leader-follower hierarchical decision model (LFHDM) for allocating initial emission rights in a basin is developed. Based on the bilevel programming approach, the model simulates the interactive decision-making process between the watershed management agency of the upper-level model (LFHDM-U) and the local government of the lower-level model (LFHDM-L) in the allocation under total amount control. A case study of China’s Yellow River Basin is conducted to demonstrate the feasibility and practicality of the model. Findings reveal that, compared with the single-level model, the developed LFHDM has higher satisfaction with the allocation scheme. Under different scenarios, the overall satisfaction of the configuration schemes of COD and NH3-N in each province and autonomous region remains above 0.9. In addition, the allocation volumes of COD and NH3-N in each province of the Yellow River Basin in planning year increase with the enhancement of allowable assimilative capacity of water bodies, but the interval gap of satisfaction with allocation schemes gradually narrows. It shows that when the allowable assimilation capacity of a water body is low, the decision-making of the allocation scheme needs to be more cautious. Moreover, for the Yellow River Basin, apart from Qinghai and Sichuan, the task of reducing water pollutants in other provinces in the next few years is very arduous. The average reduction of total COD and NH3-N in the basin is about 48% and 46%, respectively.


Introduction
Environmental problems are one of the main types of problems faced by human society in the 21st century. When the pollution caused by human social and economic activities increasingly breaks through the carrying capacity of the ecological environment, environmental problems will become an obstacle to the sustainable development of each country [1]. In particular, as a large category of environmental pollution, the pollution of the water environment is affected by the fluidity of water resources and cross-regional characteristics, which have strong complexity and also restrict the economic development of the region. The 2018 United Nations World water development report pointed out that since the 1990s, the water quality of rivers in Africa, Asia and Latin America has Int. J. Environ. Res. Public Health 2023, 20, 1511 2 of 25 been deteriorating, and this deterioration trend will continue for decades in the future, threatening human life and health and the harmonious development of social economy and ecological environment [2]. Among them, developing countries are facing the most serious threat of water quality deterioration due to the rapid development of population and social economy.
Facing the severe environmental situation, countries around the world have taken different measures to regulate the discharge of pollutants such as sewage, exhaust gas and solid waste, and emission rights are widely used as an effective measure of environmental governance [3,4]. Emission right, also known as pollution right, was first proposed by the American economist John Dales in 1968. It is defined as the right of the polluter to discharge certain pollutants into the environment according to the obtained emission indicators due to the production or domestic emission needs. This right is based on the property right of environmental capacity resources and stipulates the time, place, quantity category and mode, which is a kind of possession and using right [5]. Montgomery further demonstrated the importance and effectiveness of the establishment of emission rights under the total amount control policy on this basis [6].
Initial emission rights allocation refers to the formulation of total amount control objectives under the condition of limited environmental capacity resources, and on this basis, the quantity of pollutants allowed to be discharged into the environment by different emission subjects is determined in the form of emission permits [7]. Emission subjects can only discharge pollutants after obtaining emission permits. Although, Coase Theorem [8] points out that when there is no transaction cost, the pareto optimal state of the resource market has nothing to do with the initial allocation of property rights; zero transaction cost is difficult to achieve in reality, so the initial emission rights allocation has always been a technical and political problem [9][10][11].
There are three main allocation modes of emission rights: free allocation, paid allocation and mixed allocation. Among the three distribution modes, the free allocation mode is that the government adopts a certain distribution method to distribute the emission rights to different subjects free of charge, such as the grandfather system [12,13], the benchmarking method [14], etc. Paid allocation usually refers to the purchase of emission rights from the government by the emission subject, mainly through public auction [15][16][17] and government pricing [18]. In the mixed allocation mode, some emission rights are distributed free of charge, and the rest are paid [19]. Compared with the latter two allocation modes, under China's current national and water conditions, the free allocation mode is more operable in the initial emission rights allocation in provinces, which can reduce the obstacles to the implementation of the initial allocation in the basin [20,21]. Although there are differences in the forms of the three allocation modes in which the emission subject obtains the quota, in the process of the initial emission rights allocation of the basin, there is a vertical relationship between the upper basin management agency or government exercising the allocation power and the lower local government or department obtaining the emission right quota [22,23]. How to balance such a vertical relationship in the initial allocation and form a master-slave hierarchical interactive decision-making is still an issue to be explored.
As such, under the total amount control, this paper designs an interactive allocation model of initial emission rights based on a bilevel programming method [24,25], which can reflect the interactive decision-making between the upper watershed management agency and the lower local government. In the upper model, the multi-objective optimization method is embedded to maximize the economic, social and ecological benefits of the basin from the perspective of equity. In the lower model, the zero-sum gains data envelopment analysis (ZSG-DEA) method [26] is used to maximize the allocation efficiency of initial emission rights under total amount control. Through the interactive solution of the upper model and the lower model, the initial emission rights allocation scheme with the greatest satisfaction is obtained. Finally, a case study of the Yellow River Basin is carried out. The three main contributions of this study are as follows: (1) Based on fairness and efficiency principles, a leader-follower hierarchical decision model of upper-and lower-level governments is proposed to allocate the initial emission rights of the basin, which can provide a reasonable allocation scheme. (2) The distribution scheme is expressed in the form of interval number, and the distribution model is solved by using the satisfaction function of comprehensive benefits in group decision making, which is conducive to solving the increasingly complex problem of initial emission rights allocation in the basin. (3) In combination with the policy background of China, the total amount control policy of water pollutants is included in the constraint conditions of the model to help decision-makers make more scientific plans.
The rest of this study is organized as follows: Section 2 provides a literature review; Section 3 introduces the leader-follower hierarchical decision model for initial emission rights allocation and describes the study area and data; the detailed results and discussion are illustrated in Section 4; Section 5 summarizes the main conclusions.

Literature Review
The initial allocation of emission rights is the premise for the smooth implementation of the emission trading policy [27]. At present, many scholars have conducted a lot of research on the allocation principle of initial emission rights, mainly focusing on the fairness principle [28][29][30], efficiency principle [31][32][33] or multi-dimensional principle [34][35][36][37], to find the optimal allocation standard for establishing the optimal allocation of initial emission rights. The fairness principle can not only make the polluters willing to accept the distribution scheme, but also maximize the overall social benefits through emission trading, taxation and other measures to achieve the goal of win-win [38]. The efficiency principle is based on economic benefits in water pollution load distribution, which can achieve the best combination of economic benefits and environmental benefits to a certain extent. Therefore, it is often favored by environmental management departments and favored by researchers. However, the simple pursuit of fairness or the principle of maximizing efficiency can no longer meet the increasingly scarce initial emission right allocation. More scholars have actively explored and effectively studied it by combining fairness and efficiency principles [39][40][41].
Based on considering the allocation principle, the previous research on the allocation model of initial emission rights can be divided into four branches. In the first branch, the single objective decision-making model is widely used [42]. For example, Kang, Li [43] analyzed the impact of water quality on water resource allocation and studied the resource allocation scheme under the water quality protection goal of the water function zone based on the water resource supply model. Li, Zhang [44] determined the ammonia nitrogen and total phosphorus load distribution under environmental capacity constraints according to the equal proportion distribution model. The second branch includes studies pertaining to a multi-objective optimization model [45,46]. For example, Yandamuri, Srinivasan [47] proposed a multi-objective decision-making model for optimal allocation of water pollutants considering the total treatment cost, fairness among waste emitters and comprehensive performance indicators reflecting the violation characteristics of dissolved oxygen (DO). Liu, Guo [48] developed a waste load optimization allocation model based on a genetic algorithm with the goal of maximizing economic benefits, minimizing water shortage and maximizing waste load. The third branch includes studies about a game theory model [49,50]. Poorsepahy-Samian, Kerachian [51] used a game model to allocate water and discharge permits in the agricultural areas of the Karoon Dez River in Iran. Nikoo, Beiglou [52] developed a new cooperative game model to determine the waste load distribution strategy in rivers. The fourth branch includes studies about another game theory model. Wu, Gao [53] combined the Gini coefficient with the linear interactive general optimization method to build a coupling allocation model of ammonia nitrogen emission permit in the Songhua River Basin from the perspective of basin and region. Xie, Xu [54] used Data Envelopment Analysis (DEA) and Nash's non-cooperative game theory to build a coupled fixed cost allocation model for water pollutant discharge permits.
Overall, the initial allocation model has gradually developed from a single objective decision-making method to a coupling configuration method based on multi-objective, multi-indicator and multi-model, and the allocation principle has evolved from a single perspective of fairness or efficiency to combining fairness and efficiency. However, the following limitations still exist in the existing literature: (1) The Gini coefficient method is mostly used in the research of the emission rights allocation to reflect the principle of equity, and other equity factors between regions, such as the adaptability of economic development volume and emission rights allocation, are not considered. (2) In terms of the efficiency principle, it is easy to ignore the principle of total emission control in the existing literature, that is, to maximize the allocation efficiency under limited allocation. (3) The existing allocation model less reflects the interactive decision-making relationship between the upper and lower governments.

Problems Statement
The allocation of initial emission rights under total amount control is a complex system with multi-objective and multi-criteria, which requires interactive decision-making between the watershed management agency and local governments in each subarea. To characterize this process, the leader-follower hierarchical decision model (LFHDM) for allocating initial emission rights is formulated, which can be described by a bilevel programming approach. The decision-makers in the upper-(LFHDM-U) and lower-level (LFHDM-L) model respectively correspond to the watershed management organization and regional governments in the basin. The watershed management organization has the right to preallocate the initial emission rights of water pollutants to each regional government based on the principle of fairness. However, the regional distribution efficiency needs to be considered. On the premise of obeying the watershed management organization, the regional government will feedback information according to the allocation plan, which can help the watershed management agency better adjust the initial allocation scheme. Through leaderfollower hierarchical decision-making of the upper-and lower-level, a reasonable allocation scheme of the initial emission rights can be obtained. The framework of the LFHDM is depicted in Figure 1. basin and region. Xie, Xu [54] used Data Envelopment Analysis (DEA) and Nash's non cooperative game theory to build a coupled fixed cost allocation model for water pollutan discharge permits.
Overall, the initial allocation model has gradually developed from a single objectiv decision-making method to a coupling configuration method based on multi-objective multi-indicator and multi-model, and the allocation principle has evolved from a singl perspective of fairness or efficiency to combining fairness and efficiency. However, th following limitations still exist in the existing literature: (1) The Gini coefficient method i mostly used in the research of the emission rights allocation to reflect the principle of eq uity, and other equity factors between regions, such as the adaptability of economic de velopment volume and emission rights allocation, are not considered. (2) In terms of th efficiency principle, it is easy to ignore the principle of total emission control in the exist ing literature, that is, to maximize the allocation efficiency under limited allocation. (3 The existing allocation model less reflects the interactive decision-making relationship be tween the upper and lower governments.

Problems Statement
The allocation of initial emission rights under total amount control is a complex sys tem with multi-objective and multi-criteria, which requires interactive decision-makin between the watershed management agency and local governments in each subarea. T characterize this process, the leader-follower hierarchical decision model (LFHDM) fo allocating initial emission rights is formulated, which can be described by a bilevel pro gramming approach. The decision-makers in the upper-(LFHDM-U) and lower-leve (LFHDM-L) model respectively correspond to the watershed management organizatio and regional governments in the basin. The watershed management organization has th right to pre-allocate the initial emission rights of water pollutants to each regional gov ernment based on the principle of fairness. However, the regional distribution efficienc needs to be considered. On the premise of obeying the watershed management organiza tion, the regional government will feedback information according to the allocation plan which can help the watershed management agency better adjust the initial allocatio scheme. Through leader-follower hierarchical decision-making of the upper-and lower level, a reasonable allocation scheme of the initial emission rights can be obtained. Th framework of the LFHDM is depicted in Figure 1.

Objective Functions
In the upper-level model, the watershed management agency at upper-level makes the initial allocation decision, then the local government at lower-level executes its initial allocation plan. Under total amount control, the overall objective function of upper-level planning model (LFHDM-U) is to maximize the comprehensive benefits of emission rights allocation, which can be divided into three sub objectives: basin economic benefits, basin social benefits and basin eco-environmental benefits.
(1) Objective function of basin economic benefits The economic benefits goal is to maximize the overall economic benefit of pollutants discharge through the initial emission rights allocation among different subareas in the basin, calculated by: where superscript ± denotes the upper and lower limits of interval parameters; F U1 represents the economic benefits of pollutant d discharge in the planning year t; f i1 (R ± idt ) denotes the emission economic benefits of pollutant d at subarea i in the planning year t; n is the number of subareas, i = 1, 2, · · · , n; R ± idt is the decision variable of LFHDM-U, which denotes the allocation amount of pollutant emission rights d at subarea i in planning year t; θ i is the decision variable of LFHDM-L, which indicates the allocation efficiency of emission permits in the subarea i; BR ± idt is the unit emission revenue of the pollutant d obtained by the subarea i in the basin. Assume that the economic index of subarea i in planning year t is G it (R ± idt ), which can be represented by GDP indicators. Let the emission performance function of subarea can be solved by an exponential fitting method based on the ratio of the pollutant emission performance function V idt (R ± idt ) to R ± idt . (2) Objective function of basin social benefits The social benefits goal is to maximize the equity of emission rights allocation among subareas, so as to realize the balanced development of each subarea in the basin. The water environmental Gini coefficient (WEGC) is applied to measure the equality degree of water pollutant emission intensity loaded by the unit population index of each subarea to ensure that the allocated water pollutant emission rights match the regional population size [55,56]. The smaller the WEGC, the fairer the distribution of emission rights. The objective function of basin social benefits can be calculated by: where F U2 represents the social benefits of pollutant d discharge in the planning year t; X ± it is the cumulative proportion of population in the ith subarea of watershed in planning year t, which is expressed as represents the population of the ith subarea in the watershed in the planning year t; Y ± idt is the cumulative proportion of water pollutant d of ith subarea after the allocation in planning year t, which can be where R ± idt represents the allocation amount of pollutant d of regional pollutants in the planning year t. When i = 1, both X (i−1)t and Y d(i−1) can be regarded as (0, 0). The eco-environmental objective is to ensure the minimization of the amount of pollutants discharged in the river basin through the initial emission rights allocation of each subarea, so as to realize the optimization of the eco-environmental benefits of the river basin, calculated by: where F U3 is the eco-environmental benefits of pollutant discharge in planning year t; f i3 (R ± idt ) is the eco-environmental benefits of pollutant d at subarea i in the planning year t.
Combining the three sub objective functions, the overall objective function F U of LFHDM-U for basin initial emission rights allocation can be expressed as follows:

Constraints
Constraints for total amount of pollutant discharge, regional coordinated development, social equity, as well as nonnegativity, are considered at the upper level, given by: (1) Constraint for total amount of pollutant discharge In planning year t, the emission rights of pollutant d allocated to subarea i should be no more than its maximum allowable total emission amount.
where TR ± dth is total amount control range of pollutant d in planning year t under scenario h of allowable assimilative capacity of water bodies.
(2) Constraint for regional coordinated development To ensure the coordination of the upper-level distribution, the emission rights allocated to each subarea should match its coordinated allocation coefficient.
where R ± kdt denotes the emission rights of pollutant d obtained by subarea k in planning year t; SES ± idt and SES ± kdt are, respectively, the coordinated allocation coefficient for pollutant d of subarea i and subarea k in planning year t. The process of setting the indicates that the change direction of the ratio of coordination allocation coefficients in two subareas must be consistent with the ratio of their emission rights allocation quantity. η min represents the minimum coordination allocation coefficient, which can be determined by democratic consultation among subareas and administrative arbitration of the watershed management agency according to the characteristics of current emission rights allocation in the basin.
(3) Constraint for social justice In planning year t, the Gini coefficient of water pollutants discharge based on population index should not exceed a certain threshold.
where µ d is the Gini coefficient threshold corresponding to water pollutant d. It can be determined according to the social average level of Gini coefficient of water pollution discharge.
(4) Nonnegative constraint In the process of initial emission rights allocation in the basin, the water pollutant emission rights allocated to each subarea must be greater than zero.
The allocation efficiency of the initial emission rights of water pollutants in LFHDM-L is an important basis for judging whether the allocation plan made in LFHDM-U is reasonable. In LFHDM-L, this study applies the zero-sum gains data envelopment analysis (ZSG-DEA) model to estimate the allocation efficiency of each subarea.
According to the classic DEA model, each subarea can be regarded as a decisionmaking unit (DMU). For any subarea, its input (or output) variables have high degrees of freedom and are completely independent of each other. The change of input (or output) of the subarea will not affect the input (or output) of other subareas. However, in the field of resource allocation, the input (or output) variables of a subarea will be subject to the constraints of its total input and usually cannot obey the assumption of complete independence between variables. Moreover, under the total amount control in the basin planning year, the distribution of initial emission rights among subareas is relevant. In other words, if a subarea reduces the unexpected output to improve the allocation efficiency, it will inevitably affect other subareas.
The ZSG-DEA model can well overcome the problem that classic DEA cannot meet the established constraints of total amount, reallocate resources and change the frontier corresponding to all DMUs, but keep the sum of allocated variables unchanged [57]. Based on this method, this paper calculates the allocation efficiency of initial emission rights obtained by LFHDM-U and redistributes the inefficient units.
We assume that the emission rights of water pollutant d are allocated to each subarea DMU i in planning year t by the LFHDM-U model. As an undesirable output, the allocated emission rights quota R ± idt can be used as the only input variable of the ZSG-DEA model to calculate the allocation efficiency [58]. The population, water consumption and gross regional product of each subarea are taken as output variables, which can be expressed by a, b and c, respectively. If the distribution efficiency of subarea DMU 0 in LFHDM-L is inefficient, in order to improve the efficiency, the allocated emission rights quota R ± 0dt of subarea DMU 0 must be reduced by: where R ± Z0dt is the reduced emission rights quota of subarea DMU 0 ; θ Z0 ± is the distribution efficiency of subarea DMU 0 . When θ Z0 ± = 1, the allocation of emission rights in subarea DMU 0 is effective; otherwise, it is invalid.
Following the proportion increase strategy [59], all the other subareas DMU i (i = 0) share the reduced emission rights quota of subarea DMU 0 according to the initial input proportion; that is, the increase of emission rights quota obtained by subarea DMU i from subarea DMU 0 can be measured as: where R ± idt+Z0 is the increased emission rights quota of subarea DMU i .
After all adjustments are completed, the redistribution quota R ± idt of emission rights for DMU i is calculated by: where θ Zi ± is the distribution efficiency of subarea DMU i ; R ± idt (1 − θ Zi ± ) represents the emission rights quota reduced by subarea DMU i when the θ Zi ± is less than 1. According to the input-oriented model of ZSG-DEA, the objective function F L of LFHDM-L is to maximize the allocation efficiency of emission rights in each subarea, which can be expressed as follows: where θ Z0 ± is the relative allocation efficiency value of subarea DMU 0 . The greater the θ Z0 ± , the closer the DMU 0 is to the frontier of ZSG-DEA. There is a linear relationship between the ZSG-DEA model and classic DEA model when a single input is competitive [59], which can be characterized as follows: where θ basic ± is the relative technical efficiency calculated by the classic DEA model; W is a set of DMUs whose relative technical efficiency calculated by the classic DEA model is not 1; ρ i0 is the technical efficiency ratio between subarea DMU i and subarea DMU 0 .

Constraints
Constraints for total distributable amount, output, convex, as well as nonnegativity, are considered at the lower level.
(1) Constraint for total distributable amount The total amount of the emission rights added by DMU i with effective allocation of emission rights shall not exceed the amount of emission rights allocated by invalid DMU 0 .
where R ± idt is the input index value of DMU i ; λ i is the weight coefficient, which represents the allocation ratio of DMU i in all DMU(excluding DMU 0 ).
(2) Output constraint In planning year t, the output value of effective combination DMU i should be greater than or equal to that of invalid combination DMU 0 .
where a it ± , b it ± and c it ± are the output index values of DMU i ; a 0t ± , b 0t ± and c 0t ± are the output index values of DMU 0 .

Solution of Model
Based on LFHDM-U and LFHDM-L, the problem can be represented as a leaderfollower hierarchical decision model for initial emission rights allocation of the water pollutants in a basin. Combining the objective functions and constraints, the global IPBM is formulated as follows: LFHDM-U model : LFHDM-L model : In the IPBM, the objective function and constraints of LFHDM-U and LFHDM-L are interrelated and restricted. The optimal solution of LFHDM-U for watershed management organization depends not only on the upper decision-making variables, but also on the optimal solution of LFHDM-L for regional governments in each subarea. The upper decision variables also have an impact on the Optimization of LFHDM-L. To this end, a leader-follower hierarchical interactive iterative algorithm can be constructed to find the relative optimal solution based on satisfying the double-layer objectives [23].
The specific steps are as follows: Step 1: Set initial round K = 1.
Step 2: Under total amount control, the watershed management organization determines the first-round allocation scheme R ± idtK of initial emission rights according to the constraints of LFHDM-U.
Step 3: Substitute the first-round allocation scheme R ± idtK into LFHDM-L and calculate the initial allocation efficiency θ Z0 ± . If the allocation efficiency of all subareas reaches the effective frontier (i.e., the θ Z0 ± = 1), proceed to Step 5; if the allocation scheme of a certain subarea is invalid (i.e., the θ Z0 ± = 1), it shall be adjusted by the iterative method to obtain the new allocation plan R ± idtK under the condition of maximizing the technical efficiency of each subarea.
Step 4: Substitute the R ± idtK and θ Z0 ± into the objective function of LFHDM-U and judge whether it meets their constraints. If not, set the second round (i.e., K = K + 1), return to Step 2 and perform a new round of allocation according to the constraints of LFHDM-U. Otherwise, the value of basin economic benefits F U1 , basin social benefits F U2 and basin eco-environmental benefits F U3 can be determined.
Step 5: According to the determined values of F U1 , F U2 and F U3 , the satisfactory-degree function is formulated. It is a population decision-making method, which can combine the opinions of government personnel, experts in the field of water resources and mass representatives to evaluate the comprehensive benefit satisfaction of the initial emission rights allocation plan in the basin and further optimize the allocation plan for regional governments. The function can be expressed as follows: where µ K (F U ) denotes the comprehensive benefit satisfaction of the Kth round of emission rights distribution; µ K (F U1 ) is the economic benefit objective satisfaction function; µ K (F U2 ) is the social benefit objective satisfaction function; µ K (F U3 ) is the eco-environmental benefit objective satisfaction function; W 1 ,W 2 and W 3 represent the weights of the three single objective functions. Since the three single objective benefits jointly affect the allocation of emission rights, the weights are equally divided (i.e., ] * indicate the expected target level values, which can be determined in conjunction with the national development policy report and the opinions of government agencies, experts in the field of water resources and representatives of the masses. Step 6: A multi-agent organization composed of watershed management institutions, third-party research institutions and environmental protection public welfare organizations evaluates the comprehensive benefit satisfaction µ K (F U ). If the µ K (F U ) is greater than 90% and the constraints are met, the scheme adjustment will be stopped and the obtained distribution results are the final satisfactory solutions of the model. Otherwise, set K = K + 1, and go to Step 7.
Step 7: Take the µ K (F U ) as an added constraint (i.e., µ K+1 (F U ) ≥ µ K (F U )), so that the comprehensive benefit of the initial distribution is always on the optimization path. Then, go to Step 2 and carry out a new round of allocation according to the constraints of LFHDM-U. The average annual precipitation in YRB is about 400 mm, while the average annual runoff is only 58 billion cubic meters, accounting for 2% of the total river runoff nationwide.
YRB is one of the most ecologically and economically valuable rivers in China. The cultivated land irrigated by water sources in the basin accounts for 15% of China's farmland and provides water for 12% of the population. Its annual gross regional product (GDP) is nearly 8 trillion yuan, accounting for about 14% of China's GDP. However, as the YRB flows through coal-and oil-producing areas such as Inner Mongolia (NM), Shanxi (SX) and Shaanxi (SaX) and there are many high-energy-consumption and high-polluting enterprises, the amount of sewage flowing into the basin is very large. With the increase of water demand, the discharge of water pollutants also increases significantly, and the watershed water environment management faces important challenges.

Data Collection and Parameter Identification
(1) Data collection Considering the availability of data, as well as the water environment status and water pollution characteristics of the Yellow River Basin, the control indicators of water pollutant emission rights in the basin are defined as COD and NH 3 -N. The reasons are as follows. The pollutant control indicators in the Yellow River Basin mainly include COD, NH 3 -N, TP and suspended solids. Among them, COD and NH 3 -N are the key control indicators to eliminate the black and odorous water in the river [60]. The statistics of these two indicators are mainly carried out in the statistical data of water pollutants in the Yellow River Basin, so the availability of data and data is the highest. Due to the different planning scope and calculation caliber of the pollutant carrying capacity of the Yellow River Basin in different planning schemes, the total amount control value of major water pollutants in 2030 is not unified. In addition, the pollutant carrying capacity of the Yellow River Basin in the planning year is affected by the change trend of the annual water volume and major pollutants emission in the historical statistical interval of the basin. The pollutant carrying capacity of the water area has a positive correlation with the basin water volume and its distribution and a negative correlation with the annual discharge and distribution of major pollutants in the historical statistical interval. Therefore, according to the different water inflow and pollutant discharge in the planning year, this paper divides the total amount control target values under the three scenarios h(h = 1, 2, 3) of water bodies' allowable assimilative capacities. Under the scenario h = 1, the allowable assimilative capacity of water bodies is low and the water pollutant carrying capacity control is high; that is, the total amount control target value is low. Under the scenario h = 2, the allowable assimilative capacity of water bodies is moderate and the pollutant carrying control of water pollutants is moderate; that is, the total amount control target value is moderate. Under the scenario h = 3, the water area has high allowable assimilative capacity and low water pollutant carrying control; that is, total amount control target value is high. Based on the statistics data and the missing value supplement of the statistical caliber of pollution, the total amount control target of pollutant emission limitation in the water function area of the Yellow River Basin in 2030 is expressed by interval numbers according to scenarios. Table 1 shows the total amount control target values of COD and NH 3 -N under different scenarios of allowable assimilative capacity of water bodies in the Yellow River Basin in 2030. (2) Parameter identification The research purpose of this paper is to obtain the initial emission rights allocation scheme of the Yellow River Basin in 2030, so the planning year t can be set as 1. For convenience of expression, it is omitted in the following text. Based on the Comprehensive Planning for Yellow River Basin (2012-2030) and Yellow River Water Resources Bulletin, the BR ± i1 of COD and the BR ± i2 of NH 3 -N are calculated by using the exponential fitting method (see Appendix B). The SES ± i1 for COD and the SES ± i2 for NH 3 -N in each region are estimated according to the characterization indicators (see Appendix C).

Solutions of the LFHDM Model
First, in the scenario of h = 1 where the water pollutant is COD, we analyze the lower-limit initial allocation interval of emission rights of each province in the Yellow River in 2030. In this paper, the population, regional gross domestic product (GDP), regional area and water resource consumption are selected as output indicators, and the quota of COD into rivers and lakes as input indicators [54]. Since the initial configuration scheme in LFHDM-U is interval number, according to interval DEA theory [61], the interval ZASG-DEA model can be transformed into deterministic ZASG-DEA when setting LFHDM-L. When solving the optimal efficiency value, the lower interval bound of the input index is taken as the input value of the ZASG-DEA model, and the upper interval bound of the output index is taken as the output value. In this case, the deterministic programming model of the upper limit for optimal efficiency value can be obtained. Otherwise, the deterministic programming model of the lower limit for optimal efficiency value can be obtained according to the upper interval bound of the input index and the lower interval bound of the output index. With the leader-follower hierarchical decision model for initial emission rights allocation of the water pollutants in the Yellow River Basin, the cyclic coupling distribution of COD emission right between subareas under total amount control can be set. According to the requirements of coordinated matching between emission rights allocation and social economic development, the minimum matching coefficient η min of SES ± idt is taken as 0.8. In addition, since the Gini coefficient value is better to be lower, by the threshold limit of the environmental Gini coefficient, the upper limit of the Gini coefficient µ of water pollutant discharge is set to 0.4.
As it is analyzed in Section 3.4, when K = 1, the initial solution R i1 − of LFHDM-U for COD in 2030 under the scenario h = 1 is shown in Table 2. Combined with the initial solution and the corresponding efficiency value, the F 1 (R − 1 ), F 2 (R − 1 ) and F 3 (R − 1 ) are calculated as 596.740, 0.236 and 58.63, respectively. We take the emission performance when the maximum emission demand is reached in the planning year as the economic expected target value and make [F 1 (R ± dt )] * = 919.625. The Gini coefficient of 0.3 is taken as the social expectation target [F 2 (R ± dt )] * , which is mainly because the setting of the expectation target threshold is too high and may not be applicable to the initial allocation of emission rights in the development stage. Moreover, the ecological target value [F 3 (R ± 1 )] * is the lower emission limit of COD with the lowest and strictest pollutant carrying capacity in the water area, which is equal to 58.63. The setting method of the expectation target level of NH3-N is similar. According to Equation (20), the comprehensive benefit satisfaction µ K=1 F(R − 1 ) is calculated as 0.913.
In the first round, the initial DEA efficiency of the emission allocation efficiency of Ningxia and Shanxi do not meet the technical efficiency frontier, so the allocation plan should be adjusted by the ZSG-DEA model. After two rounds of iteration, the efficiency of emission rights allocation in all provinces (autonomous regions) has reached the frontier. However, when the adjusted allocation scheme is fed back to LFHDM-U, it cannot meet the constraint for regional coordinated development, so it needs to enter the second round of solution. When K = 2, a new round of adjustment and allocation of emission rights in the Yellow River Basin is carried out according to the constraints of LFHDM-U. After the efficiency solution, the allocation of all provinces (autonomous regions) is technically effective, and the iteration is completed. Table 2 displays the iteration results and efficiency values of COD allocation lower limit under 2030 scenario h = 1. After two rounds of iteration, the lower limit allocation scheme of COD emission rights under scenario h = 1 satisfies the constraints of LFHDM-U. According to Table 2, the F 1 (R − 1 ), F 2 (R − 1 ) and F 3 (R − 1 ) in the second iteration are calculated as 590.327, 0.189 and 58.63, respectively. This corresponding comprehensive benefit satisfaction µ K=2 F(R − 1 ) is 0.933, which is higher than 0.913 in the initial solution. Referring to the above calculation process, all initial distribution schemes of COD and NH 3 -N under scenario h = 1, h = 2 and h = 3 are calculated. The calculation results are given in Table 3.

Comparisons with Single-Level Models
Taking the lower limit allocation of COD in the Yellow River Basin under scenario h = 1 as an example, the solution results of LFHDM-U, LFHDM-L and IPBM are illustrated in Table 4. For LFHDM-U, it aims to optimize the overall benefits of the basin. According to Table 4, the comprehensive benefit satisfaction in LFHDM-U is 0.911, which meets the trend of regional coordinated development. However, it cannot guarantee that the emission efficiency of each subarea in the basin is effective, and the emission demand of some regions is sacrificed in the allocation. For example, the emission technical efficiency of Ningxia and Shanxi has not reached the forefront of ZSG-DEA, and the allocation efficiency is insufficient. In the implementation of the initial emission rights allocation policy, some provinces (autonomous regions) would be likely to respond negatively to the policy.
For LFHDM-L, it only considers the emission efficiency of each province. The emission efficiency of each province is effective, and the economic benefit is better than the other two models. However, its emission rights allocation scheme does not meet the evaluation trend of regional coordinated development in Equation (6), which is not conducive to reducing the development gap between provinces. It can easily lead to the well-developed provinces having more emission rights and more serious pollution, while the less-developed provinces have less emission rights and more backward development, which cannot reflect the fairness of the initial distribution of emission rights.
However, LFHDM takes into account the overall benefits of the basin and the emission efficiency of each subarea. Through the continuous adjustment of the initial distribution scheme, its comprehensive benefit satisfaction reached 0.933, which was better than 0.913 of LFHDM-U and 0.918 of LFHDM-L and met the constraint for regional coordinated development.

Allocation and Variation Trend of Initial Emission Rights
According to the results in Table 3, the variation trend charts of allocation interval amount for COD and NH 3 -N under three scenarios are drawn, as shown in Figures 2 and 3.

Int. J. Environ. Res. Public Health 2023, 20, x FOR PEER REVIEW 15 of 26
For LFHDM-U, it aims to optimize the overall benefits of the basin. According to Table 4, the comprehensive benefit satisfaction in LFHDM-U is 0.911, which meets the trend of regional coordinated development. However, it cannot guarantee that the emission efficiency of each subarea in the basin is effective, and the emission demand of some regions is sacrificed in the allocation. For example, the emission technical efficiency of Ningxia and Shanxi has not reached the forefront of ZSG-DEA, and the allocation efficiency is insufficient. In the implementation of the initial emission rights allocation policy, some provinces (autonomous regions) would be likely to respond negatively to the policy.
For LFHDM-L, it only considers the emission efficiency of each province. The emission efficiency of each province is effective, and the economic benefit is better than the other two models. However, its emission rights allocation scheme does not meet the evaluation trend of regional coordinated development in Equation (6), which is not conducive to reducing the development gap between provinces. It can easily lead to the well-developed provinces having more emission rights and more serious pollution, while the lessdeveloped provinces have less emission rights and more backward development, which cannot reflect the fairness of the initial distribution of emission rights.
However, LFHDM takes into account the overall benefits of the basin and the emission efficiency of each subarea. Through the continuous adjustment of the initial distribution scheme, its comprehensive benefit satisfaction reached 0.933, which was better than 0.913 of LFHDM-U and 0.918 of LFHDM-L and met the constraint for regional coordinated development.

Allocation and Variation Trend of Initial Emission Rights
According to the results in Table 3      According to Figures 1 and 2, in 2030, the allocation amount of COD and NH3-N in each province (autonomous region) of the Yellow River Basin will gradually increase with the enhancement of the allowable assimilative capacity of water bodies. Shaanxi has the highest allocation of COD and NH3-N among the nine provinces, while SC has the lowest allocation of emission rights due to its small area and population in the basin. The allocation amount of COD in other provinces is Gansu, Shanxi, Henan, Inner Mongolia, Qinghai, Shandong and Ningxia in descending order. In addition, the allocation amount of NH3-N in other provinces is Gansu, Shanxi, Inner Mongolia, Henan, Qinghai, Shandong and Ningxia in the order from small to large.
In addition, shown in Figure 4, under different scenarios, the comprehensive benefit satisfaction of allocation schemes for COD and NH3-N in various provinces (autonomous regions) remains above 0.9, indicating that the leader-follower hierarchical decision-making allocation scheme is relatively reasonable. When 1 h = , it means that the water area has low allowable assimilative capacity and strict discharge restriction in the basin, but the satisfaction of the distribution scheme for COD and NH3-N is the highest, with an average value of 0.929 and 0.924, respectively. With the increase of allowable assimilative capacity of water bodies, the satisfaction of the optimal allocation scheme decreases when the total amount control is gradually relaxed. This is mainly because the satisfaction of ecological benefits reduced by the increase of water pollutants discharge is higher than that of economic benefits increased. Moreover, with the change of different scenarios, from scenario 1 to scenario 3, the interval gap of the satisfaction for the allocation scheme gradually narrows, indicating that when the emission restriction policy is more stringent, the change of the allocation amount of water pollutant emission rights has a greater impact on the satisfaction of the comprehensive benefits and the uncertainty of the allocation scheme decision increases. Conversely, when the emission restriction policy is more relaxed, the change of the assigned amount interval has less impact on the comprehensive benefit satisfaction and the uncertainty of the allocation scheme decision is relatively small. Therefore, when the allowable assimilative capacity of water bodies in the Yellow River Basin is low and emission limitation is strict, the allocation of initial emission rights among regions needs more careful decision-making.  According to Figures 1 and 2, in 2030, the allocation amount of COD and NH 3 -N in each province (autonomous region) of the Yellow River Basin will gradually increase with the enhancement of the allowable assimilative capacity of water bodies. Shaanxi has the highest allocation of COD and NH 3 -N among the nine provinces, while SC has the lowest allocation of emission rights due to its small area and population in the basin. The allocation amount of COD in other provinces is Gansu, Shanxi, Henan, Inner Mongolia, Qinghai, Shandong and Ningxia in descending order. In addition, the allocation amount of NH 3 -N in other provinces is Gansu, Shanxi, Inner Mongolia, Henan, Qinghai, Shandong and Ningxia in the order from small to large.
In addition, shown in Figure 4, under different scenarios, the comprehensive benefit satisfaction of allocation schemes for COD and NH 3 -N in various provinces (autonomous regions) remains above 0.9, indicating that the leader-follower hierarchical decision-making allocation scheme is relatively reasonable. When h = 1, it means that the water area has low allowable assimilative capacity and strict discharge restriction in the basin, but the satisfaction of the distribution scheme for COD and NH 3 -N is the highest, with an average value of 0.929 and 0.924, respectively. With the increase of allowable assimilative capacity of water bodies, the satisfaction of the optimal allocation scheme decreases when the total amount control is gradually relaxed. This is mainly because the satisfaction of ecological benefits reduced by the increase of water pollutants discharge is higher than that of economic benefits increased. Moreover, with the change of different scenarios, from scenario 1 to scenario 3, the interval gap of the satisfaction for the allocation scheme gradually narrows, indicating that when the emission restriction policy is more stringent, the change of the allocation amount of water pollutant emission rights has a greater impact on the satisfaction of the comprehensive benefits and the uncertainty of the allocation scheme decision increases. Conversely, when the emission restriction policy is more relaxed, the change of the assigned amount interval has less impact on the comprehensive benefit satisfaction and the uncertainty of the allocation scheme decision is relatively small. Therefore, when the allowable assimilative capacity of water bodies in the Yellow River Basin is low and emission limitation is strict, the allocation of initial emission rights among regions needs more careful decision-making.

Reduction Plan of Emission Rights in Each Province
Considering the current situation of COD and NH3-N emissions in the Yellow Basin in 2018, the pollutant reduction plan in the basin in 2030 is very difficult. The age reduction range of COD and NH3-N under different allowable assimilative cap of water bodies are about 48% and 46%, respectively. However, the Comprehensive ning for Yellow River Basin (2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027)(2028)(2029)(2030) has proposed that the total amount of CO NH3-N in the water function area of the Yellow River Basin in 2030 should be reduc about 70% compared with the current situation. This is mainly because this plan 2007 as the status quo year. In recent years, with the implementation of relevant en mental policies, the discharge of water pollutants in the basin has been controlle certain extent, has a downward trend year by year and the reduction pressure h creased compared with 2007. Tables 5 and 6, respectively, show the COD and NH3 duction plans of nine regions in the Yellow River Basin in 2030 under different scen

Reduction Plan of Emission Rights in Each Province
Considering the current situation of COD and NH 3 -N emissions in the Yellow River Basin in 2018, the pollutant reduction plan in the basin in 2030 is very difficult. The average reduction range of COD and NH 3 -N under different allowable assimilative capacities of water bodies are about 48% and 46%, respectively. However, the Comprehensive Planning for Yellow River Basin (2012-2030) has proposed that the total amount of COD and NH 3 -N in the water function area of the Yellow River Basin in 2030 should be reduced by about 70% compared with the current situation. This is mainly because this plan takes 2007 as the status quo year. In recent years, with the implementation of relevant environmental policies, the discharge of water pollutants in the basin has been controlled to a certain extent, has a downward trend year by year and the reduction pressure has decreased compared with 2007. Tables 5 and 6, respectively, show the COD and NH 3 -N reduction plans of nine regions in the Yellow River Basin in 2030 under different scenarios.  It can be seen from Tables 5 and 6 that the current COD emissions of the nine provinces (autonomous regions) in the Yellow River Basin in 2018 are ranked from high to low as Shaanxi, Shanxi, Inner Mongolia, Henan, Gansu, Ningxia, Shandong, Qinghai and Sichuan. The order of highest to lowest NH 3 -N emission is Shaanxi, Shanxi, Gansu, Henan, Inner Mongolia, Ningxia, Shandong, Qinghai and Sichuan. The three provinces with the largest COD emissions in the whole basin are Shaanxi, Shanxi and Inner Mongolia, accounting for 53.2% of the total. In 2018, the three provinces with the largest current emissions of NH 3 -N were Shaanxi, Shanxi and Gansu, accounting for 56.8% of the total. However, the three provinces with the lowest emissions of COD and NH 3 -N in the whole basin are Sichuan, Qinghai and Shandong. The COD and NH 3 -N emissions of the three provinces account for 9.7% and 9.3% of the whole basin, respectively.
Based on the current annual pollutant discharge, the reduction rate of water pollutants in the Yellow River Basin in 2030 is analyzed. Although the pollutant emissions of Shaanxi and Shanxi are among the top two in the basin, the pollutant reduction rates of the two provinces are not the highest in the basin under the three scenarios of allowable assimilative capacities of water bodies in the planning year. This is mainly because the two provinces also account for a large proportion of the population in the basin, and the economic development level and water pollution control efficiency are higher than other provinces. From the perspective of fairness and efficiency of distribution, the reduction task is more reasonable. Instead, Ningxia's COD and NH 3 -N emissions are only ranked sixth in the basin, but the reduction rate is the highest. This is mainly because the amount of water resources (only about 8% of the water resources in the basin) and land area in Ningxia account for a relatively small proportion in the basin, the economic development is relatively backward and the efficiency of pollution control is not high. As a result, its initial emission rights quota in 2030 is less, and the reduction pressure is at the top of all provinces and regions. However, due to the abundant total amount of water resources and low population density, Qinghai's current annual discharge capacity has been lower than the initial emission quota of water pollutants under the low allowable assimilative capacity of water bodies in the planning year. Therefore, the initial emission quota of the two pollutants in Qinghai in 2030 can basically meet the emission demand under three scenarios in the planning year. Sichuan accounts for a small proportion of the Yellow River Basin, and the current discharge of pollutants is small, which has a small impact on the overall discharge of water pollutants in the Yellow River Basin. Therefore, its emission rights quotas of COD and NH 3 -N in 2030 also meet the regional economic development and emission requirements under the three scenarios.
From the perspective of the three scenarios of allowable assimilative capacity of water bodies, when the scenario is h = 1, the provinces where the lower limit of the COD reduction rate range reaches more than 50% in 2030 include Ningxia, Shanxi, Inner Mongolia and Shaanxi. Except for Sichuan and Qinghai, the lower limit of the COD reduction rate in the other seven province exceeds one-third of the current year. When the scenario is h = 2, the provinces with the lower limit of COD reduction rate range of more than 50% in 2030 include Ningxia, Shanxi and Inner Mongolia. Except for Sichuan and Qinghai, only Gansu Province has the lower limit of the reduction rate range below onethird (30.0%). As for the reduction rate of NH 3 -N, the provinces with the lower limit of the interval above 50% include Ningxia, Shanxi and Shaanxi. Except for Sichuan and Qinghai, only Inner Mongolia has the lower limit of the cut rate below one-third (31.2%). When the scenario is h = 3, the province with the lower limit of the COD emission reduction rate range above 50% in 2030 is Ningxia (Inner Mongolia and Shaanxi are close to 50%, and their reduction rates are 49.2% and 49.1%, respectively). Except Sichuan and Qinghai, only Gansu province has the lower limit of the reduction rate range below one-third. For the NH 3 -N reduction rate, the provinces with the lower limit of the interval above 50% include Ningxia and Shanxi. Except for Sichuan and Qinghai, only the lower limit of the reduction rate range in Inner Mongolia is less than one-third. Therefore, except for Sichuan and Qinghai, most other provinces in the Yellow River Basin are facing more than one-third of the reduction rate for water pollutants emission rights under three scenarios, and the pressure of emission reduction in the next few years is huge.

Conclusions
The initial emission rights allocation in the basin is not only an important prerequisite for the implementation of emission trading, but also the most controversial and difficult issue in the process of the utilization of emission rights. Therefore, under total amount control, this study established the leader-follower hierarchical decision model (LFHDM) for allocating initial emission rights in the basin. Through the interactive decision-making between the watershed management agency and local governments in each subarea, the initial emission rights allocation scheme of the basin that takes into account fairness and efficiency can be obtained. The LFHDM model is demonstrated and tested in a case study of China's Yellow River Basin in different allowable assimilative capacities of water bodies set by allocating initial emission rights among its nine provinces (autonomous regions). Some conclusions and policy implications can be drawn from this study: (1) The comparison between LFHDM and two single-level models (i.e., LFHDM-U and LFHDM-L) shows that the overall satisfaction of LFHDM 's allocation scheme is better than that of LFHDM-U and LFHDM-L. Compared with the distribution scheme of LFHDM, LFHDM-U cannot guarantee that the emission efficiency of all provinces in the basin is effective. The distribution scheme of LFHDM-L is not conducive to narrowing the development gap between provinces and cannot reflect the fairness of the initial distribution of emission rights. When formulating the initial pollution discharge rights allocation plan for water pollutants, the watershed management agency must take into account the principles of fairness and efficiency and pay attention to communication with lower-level local governments. Otherwise, it can easily lead to a negative response to the distribution policy, and it is not conducive to the establishment of the emission rights trading market.
(2) In the planning year, the allocation volumes of COD and NH 3 -N in the provinces of the Yellow River Basin gradually increase with the enhancement of the allowable assimilative capacity of water bodies, and the interval gap in the satisfaction of the allocation scheme gradually narrows. When the Yellow River Basin has low allowable assimilative capacity of water bodies and strict emission restrictions, the allocation of initial emission rights between regions needs more careful decision-making. The initial emission rights allocation scheme of water pollutants under different scenarios shall be formulated to cope with the change of annual runoff of the basin and avoid excessive initial emission rights in some regions, which will lead to over-discharge of water pollutants.
(3) The task of reducing pollutants in the Yellow River Basin in 2030 is very arduous. Except for Sichuan and Qinghai, most other provinces in the Yellow River Basin face more than one-third of the emission reduction rate under the three scenarios of allowable assimilative capacities of water bodies. The average reduction range of the total amount of COD and NH 3 -N in the basin under the three scenarios is about 48% and 46%, respectively. Before entering the emission trading market, all regions should improve the emission reduction technology as much as possible so that regions with high emission reduction pressure can reduce the water pollutant emissions and avoid being unable to meet the local pollution demand due to too little initial emission rights allocation. The regions with low emission reduction pressure can save more initial emission rights and transfer them to regions with scarce emission rights in the trading market.
There are still some limitations in this study. Due to the limited availability of data, this study did not include the impact of uncertain factors such as climate, hydrology and human activities on the model. Future research will further explore the allocation scenarios under different influencing factors. In addition, all industries will be included in the leader-follower hierarchical decision-making model (LFHDM) for multi-level interactive decision-making.
Author Contributions: Q.Y. designed the research and drafted the manuscript; Z.S. and J.S. revised the paper. X.X. and X.C. collected the data and conducted the model simulation. The contributions of X.X. and X.C. are the same, so both are recognized as corresponding authors. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available in article and Appendices A-C.

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

Appendix A
According to [62], the characteristic indexes of SES ± idt are selected and determined by the technique for order preference by a similarity to ideal solution (TOPSIS) decisionmaking method with interval numbers [63]. Table A1 provides the characterization indicators of SES ± idt . Table A1. Characterization indicators of SES ± idt .

Indicators Unit Explanation Attribute
Population size 10,000 persons It reflects the population differences within each region. The larger the value, the more people enjoy equal emission rights.

Income-type
Area size km 2 It reflects the area difference under the jurisdiction of regional water sources. The larger the value, the higher the natural rationality of distribution.

Income-type
Total water resources amount 100 million tons It reflects the difference of regional water resources. The larger the value, the better the regional resource endowment.

Income-type
Water consumption per 10,000 yuan GDP ton It reflects the utilization efficiency of water resources in each region. The larger the value, the more water is used per unit of GDP.

Indicators Unit Explanation Attribute
Pollution emission per 10,000 yuan GDP ton It reflects the pollutant emission efficiency of each region. The greater the value, the greater the pollutant emission per unit of GDP.

Cost-type
Sewage reuse rate % It reflects the sewage reuse efficiency of each region. The larger the value, the higher the sewage treatment capacity of the region.

Income-type
Standard rate of water quality in water function area % It reflects the water quality status of water functional areas in each region. The larger the value, the better the regional water environment.

Income-type
The steps to solve the SES ± idt are as follows: Step 1: Standardizing the index value The attributes of the characteristic indexes for SES ± idt include income-type and cost-type. The cost-type indicators can be converted into income-type indicators through the reciprocal method. According to the algorithm of interval number [64], the index value is further standardized to eliminate the dimensional influence by using the following equation: , i = 1, 2, · · · , n; j = 1, 2, · · · , m; t = 1, 2, · · · , T; j ∈ J 1 (A1) where x ± ijt is the attribute value of index j for region i in planning year t; y ± ijt is the normalized value of x ± ijt ; J 1 is the subscript set of income-type index; n, m and T are the total amount of sample region, sample time and sample indicators, respectively; + refers to the upper limit of index interval value; − refers to the lower limit of index interval value.
Step 2: Determining positive and negative ideal points According to the normalized index value, the positive ideal point and negative ideal point are determined as follows: where Y jtPositive denotes the matrix of positive ideal points; Y jtNegative represents the matrix of negative ideal points.
Step 3: Determining the index weight Based on the interval number dispersion method [65], the objective weight w j of index j is calculated.
Step 4: Determining the relative proximity of each region to the ideal points According to the Manhattan Distance method [66], the relative proximity of each region to the ideal point can be determined as follows.
where S iNegative is the relative proximity of each solution to the negative ideal points for region i; S iPositive is the relative proximity of each solution to the positive ideal points for region i.
Step 5: Calculating the value of SES ± idt The relative proximity of each region to the ideal point is used to represent the value of SES ± idt , which can be calculated by:

Appendix B
According to the Comprehensive Planning for Yellow River Basin (2012-2030) and Yellow River Water Resources Bulletin, the ratio of GDP to COD emissions in the nine provinces (autonomous regions) in the Yellow River Basin from 2006 to 2018, i.e., G i (R ± i1 )/R ± i1 , and the ratio of GDP to NH 3 -N emissions, i.e., G i (R ± i2 )/R ± i2 , are given in Table A2. Table A2. Ratio of GDP to water pollutants emission amount for subareas in the Yellow River Basin.
Year QH SC GS NX IM Year SaX SX HN SD The exponential fitting method is used to fit the G i (R ± i1 )/R ± i1 and G i (R ± i2 )/R ± i2 of nine provinces (autonomous regions) from 2006 to 2018. By deriving the emission performance functions of COD and NH 3 -N in the nine provinces (autonomous regions), the BR ± id can be obtained, namely:

Appendix C
According to Appendix A, taking COD as an example, the interval parameter values of SES ± i1 after standardized treatment are shown in Table A3. The objective weights w j of the seven indexes determined by the interval number deviation method are: w j = (w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 ) = (0.160, 0.195, 0.112, 0.165, 0.141, 0.123, 0.105) Based on Equations (A4) and (A5), the Manhattan Distance between each province (autonomous region) and the ideal point in the Yellow River Basin is shown in Table A4. According to Table A4 and Equation (A6), results of calculating the SES ± i1 for COD and the SES ± i2 for NH 3 -N in each region are shown in Table A5.