Study area
Central Plains Urban Agglomeration(CPUA) (33°08′-36° 02 ′N, 111°08′-115° 15 ′E) is situated in the central and eastern regions of China, encompassing the middle and lower reaches of the Yellow River, including 18 prefectures and cities in Henan Province, Changzhi, Jincheng and Yuncheng in Shanxi Province, Xingtai and Handan in Hebei Province, Liaocheng and Heze in Shandong Province, Huaibei, Bengbu, Suzhou, Fuyang and Bozhou in Anhui Province, totaling 30 cities (Fig. 1). The CPUA is not only related to the high-quality economic growth of the central region, but also plays a linking role in the industrial upgrading of the eastern region and the economic development of the western region. However, the CPUA within the development of regional differentiation is obvious, the conditions of water resources endowment is poor, the amount of water resources is in serious shortage, the time distribution is extremely uneven, the contradiction between supply and demand is prominent. A series of environmental pollution problems triggered by urbanization and industrialization have made the region face the great test of balancing economic development and ecological environment protection.
Central plains urban agglomeration study area.
In the questionnaire stage, 15 experts and scholars in different fields were selected as respondents and the results were analyzed for reliability. The data on water resources and the ecological environment are get from regional water resources bulletins ( surveys; the economic and social data are obtained from regional statistical bulletins on national economic and social development (https://www.henan.gov.cn/zwgk/zfxxgk/fdzdgknr/tjxx/tjgb/).
Index system construction
This study combines expert interviews and literature research to select 35 influencing factors (as shown in Table 1) from three dimensions to analyze the key influencing factors of coordinated development of WEE.
Methods
Synergistic development of WEE coupling system is a hot research topic at present, but the internal relationship between its influencing factors has not been fully studied and needs further exploration. Delphi method is also called expert opinion method, which can give full play to experts ‘advantages and gather experts’ wisdom, so that the selected factors have certain credibility22. Delphi method was used to select representative experts to form scientific and reasonable expert advisory group. Experts were then invited to rate the impact of each factor on the WEE coupling system. DEMATEL-ISM-MICMAC method was used to identify the key factors, and the synergistic development relationship among water resources, economic society and ecological environment was analyzed through coupling coordination degree.
DEMATEL-ISM-MICMAC Method
Decision-Making Trial and Evaluation Laboratory (DEMATEL) is a methodology that leverages graph theory and matrix methods to elucidate the relationships between factors in intricate systems and evaluate the relative significance of each factor23. Interpretive Structural Modeling Method (ISM) can sort out the hierarchical structure of complex systems, so that the logical structure relationship of each influencing factor can be clearly displayed in a hierarchical structure diagram24. Matrix Impacts Cross-reference Multiplication Applied to a Classification (MICMAC) is a methodology that uses the driving-dependence matrix to divide factors into self-consistent factors, dependent factors, associated factors and independent factors to elucidate the critical function of each component within the system25. In this study, we will first use the DEMATEL method to analyze the interaction of factors affecting the synergistic development of WEE coupling system. ISM method was used to establish the layered structure of influencing factors. Finally, MICMAC method was used to determine which factors had the greatest driving force and dependence. The flowchart is shown in Fig. 2.

The specific steps are as follows:
Step 1 19 experts, including experts, scholars and relevant researchers in water resources and ecology, were invited to score the mutual influence relationship among indicators according to their understanding of the three-dimensional frameworks and the symbolic meaning of indicators. Specific scores are: 0 represents no impact, 1 represents a low impact, 2 represents a medium impact, 3 represents a high impact, and 4 represents an extremely high impact.
Step 2 Collect questionnaire data, calculate average value, and construct initial direct impact matrix A.
$$A = [a_{ij} ]_{n \times n}$$
(1)
where n is the number of influencing factors, indicating the degree of influence of factor i on factor j, diagonal indicates the influence of a factor on itself, and all values are 0.
Step 3 Compute normalized direct impact matrix N and synthetic impact matrix T23
$$N = [n_{ij} ]_{n \times n} = A/x$$
(2)
$$x = \max [\mathop {\max }\limits_{1 \le i \le n} \sum\limits_{i = 1}^{n} {a_{ij} ,} \mathop {\max }\limits_{1 \le j \le n} \sum\limits_{i = 1}^{n} {a_{ij} ,} ]$$
(3)
$$T = [{\text{t}}_{ij} ]_{n \times n} = N(E – N)^{ – 1}$$
(4)
where E is the identity matrix.
Step 4 Calculate centrality (Zi) and causation (Di)23.
$$r_{i} = \sum\limits_{i = 1}^{n} {t_{ij} }$$
(5)
$$c_{i} = \sum\limits_{j = 1}^{n} {t_{ij} }$$
(6)
$$Z_{i} = r_{i} + c_{i}$$
(7)
$$D_{i} = r_{i} – c_{i}$$
(8)
where ri denotes the degree of influence of factor i, the sum of the degrees of direct or indirect influence, the sum of the degree to which it directly or indirectly influences other factors, and ci denotes the degree to which factor i is influence, the sum of the degree to which it is influenced by other factors.
Step 5 Constructing adjacency matrix \(B = [b_{ij} ]_{(n + 1) \times (n + 1)}\) and reachability matrix M24.
$$b_{ij} = \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & {} \\ \end{array} } \\ {\begin{array}{*{20}c} 1 & {} \\ \end{array} } \\ \end{array} } \right.\begin{array}{*{20}c} {b_{ij} < \lambda } \\ {b_{ij} \ge \lambda } \\ \end{array} \begin{array}{*{20}c} {} & {i = 1,2, \cdots ,n + 1;\begin{array}{*{20}c} {} & {j = 1,2, \cdots ,n + 1} \\ \end{array} } \\ \end{array}$$
(9)
$$ M = (B + to) ^ {N + 2} = b + in) ^ for • \ '
(10)
where \(\lambda\) is used to eliminate some factors with small influence and is determined by debugging all values in the comprehensive influence matrix.
Step 6 According to the reachability matrix M, calculate \(R(S_{i} ),\;A(S_{i} ),\;C(S_{i} ),\;B(S_{i} ){\text{and}}\;E(S_{i} )\), \(i = 1,2, \cdots ,n + 1\)24.
$$R(S_{i} ) = \left\{ {s_{i} \in S\left| {m_{ij} = 1} \right.} \right\}$$
(11)
$$A(S_{i} ) = \left\{ {s_{i} \in S\left| {m_{ji} = 1} \right.} \right\}$$
(12)
$$C(S_{i} ) = R(S_{i} ) \cap A(S_{i} )$$
(13)
$$B(S_{i} ) = \left\{ {s_{i} \in S\left| {A(S_{i} ) = C(S_{i} )} \right.} \right\}$$
(14)
$$E(S_{i} ) = \left\{ {s_{i} \in S\left| {R(S_{i} ) = C(S_{i} )} \right.} \right\}$$
(15)
Step 7 According to the hierarchical classification criteria, \(C(S_{i} ) = R(S_{i} )\) it is the first level, then delete the used elements and repeat the classification according to the \(C(S_{i} ) = R(S_{i} )\) until all elements are classified into different hierarchical levels, resulting in a multi-level ISM model.
Step 8 Calculate driving force (Qi) and dependence (Yi). Driving force indicates the degree of influence on other factors, and dependence indicates the degree of influence by other factors25.
$$Q_{i} = \sum\limits_{i = 1}^{n + 1} {m_{ij} }$$
(16)
$$Y_{i} = \sum\limits_{j = 1}^{n + 1} {m_{ij} }$$
(17)
Among them, high driving force and high dependence factors belong to correlation factors; low driving force and high dependence factors belong to dependence factors; low driving force and low dependence factors belong to self-consistent factors; high driving force and low dependence factors belong to independent factors.
Step 9 Plot driving forces and dependencies in quadrant diagrams to visually represent relationships between factors.
Coupled coordinated development model
Different indexes have different influence on the system. This research utilizes the entropy method to assign weights to the index system. Its calculation is publicized as follows26:
$$e_{j} = – \frac{1}{\ln n}\sum\limits_{i = 1}^{n} {P_{ij} } ln\left( {P_{ij} } \right)\left( {j = 1,2, \cdots ,m} \right)\;\;\;\;\;$$
(18)
$$w_{j} = \frac{{1 – e_{j} }}{{\sum\limits_{j = 1}^{m} {\left( {1 – e_{j} } \right)} }}\left( {j = 1,2, \cdots ,m} \right)$$
(19)
where: \ (E_ {j} \) is the information entropy; \(P_{ij}\) is the proportion of the j index in the i year; is the index weight.
According to the index system of WEE, the comprehensive evaluation index of WEE is constructed, namely:
$$W\left( x \right) = \sum\limits_{k = 1}^{m} {A_{k} \cdot X_{ik} }$$
(20)
where: \(W\left( x \right)\) is the water resources system development index; \(A_{k}\) is the weight value of each index of the water resources system; \ (X_ {I} \) is the standardized value of each index of the water resources system in the i year; k is the number of indicators of the water resources system.
$$S\left( y \right) = \sum\limits_{g = 1}^{m} {A_{g} \cdot Y_{ig} }$$
(21)
where: \(S\left( y \right)\) is the development index of economic and social system; \(A_{g}\) is the weight value of each indicator of economic and social system; \ (Y_ {ig} \) is the standardized value of each indicator of economic and social system in the i year; g is the number of indicators of economic and social system.
$$E\left( z \right) = \sum\limits_{q = 1}^{m} {A_{q} \cdot Z_{iq} }$$
(22)
where: \(E\left( z \right)\) is the development index of ecological environment system; \(A_{q}\) is the weight value of each index of ecological environment system; \(Z_{iq}\) is the standardized value of each index of ecological environment in the i year; q is the number of indicators of ecological environment system.
According to the coupling concept in physics, WEE is regarded as the coupling systems, and the extent of coupling and the level of coordination between these three systems are assessed, which can objectively reflect the coordinated development of the systems27. The coupling model is:
$$C = \frac{{3\sqrt[3]{W\left( x \right)S\left( y \right)E\left( z \right)}}}{W\left( x \right) + S\left( y \right) + E\left( z \right)}$$
(23)
where C is the degree of coupling,\(C \in [0,1]\). When C = 0, it indicates that the system develops into disorder state, and when C = 1, it indicates that the three systems are in the optimal coupling state.
Coupling degree can only reflect the interaction degree of the three systems, but cannot explain their coordinated development level. Thus, a coupling coordination degree model is applied to provide an objective measure of the integrated advancement among the three systems.
$$D = \sqrt {C \times T}$$
(24)
where D is the degree of coordination,\(D \in [0,1]\).